I'm going to go back to the slides and we'll talk about example 2.
Genes within a chromosomal region.
What I would like to do is for each gene within a particular region -- I would like to output
several things.
I would like to output symbol, chromosomal accession, start and stop coordinates and
the strand.
I have a particular region in mind for a particular chromosome.
Maybe I have two genes in mind.
And I'd like to have all the genes between those two.
Perhaps I have coordinates of the chromosome in mind.
I would like to get all that genes in that segment of the chromosome.
I need to get the genes symbol.
I need to get the accession for the chromosome and the coordinates.
Where will I get this information?
This is a first task you have to give yourself, is look through the E-Utilities output for
Gene to see where you can find information.
What I have put here is a snippet from a docsum from Gene.
Within the genes docsum, it tells me the chromosome, the chromosome accession.version; you will
realize that this is human chromosome 2.
Here is the start and stop coordinates, and the number of exons.
There is no information about strand.
I need to know the strand.
We can think about that.
The strand is implicit in the start and stop.
If the start is less than the stop, the strand is forward; if the start is greater than the
stop, it is on the reverse strand.
We can use some of that logic to answer the question about strands.
The way that we can think about this is we would like to get the information from the
Gene docsum and take the information out and use of little logic.
I would like to be able to compare the start and stop of the chromosome and make a decision.
It is either strand one or strand two based on which comes first.
We can use things like awk.
Jonathan provided this is an example.
We will look at how we can use awk.
I'm gonna go back one slide.
The gene symbol is not here.
The gene symbol is in another part of the docsum XML and it is called Name.
One of the nice things that EDirect can do is temporarily store values in tags with a
variable you make up.
We can store the name in one variable and output it with the information from the start
and stop coordinates.
We will show you how to do that.
Variables are cool.
The way you do variables in xtract, it is -pattern DocumentSummary.
The document summary is that thing around the article itself.
If you put something in all capitals, xtract will create a variable name and the name of
that variable is & and the name of what it is.
It will store the value of that tag in that variable for later use, within a xtract call.
We have talked about block which allows us to parse the subsection of XML.
(back one slide) This will be the genomic info type.
Here is the overall method.
We will have a search for the particular genes of interest.
We will get the genes from esearch.
It will give us a set of gene IDs.
They are unique IDs.
We will go to esummary to get the docsums for those genes.
That is where the gold is; that is the stop and stop coordinates, accessions are.
We will use xtract to take out the pieces we want, which are Name, ChrAccVer, ChrStart
and ChrStop.
Then we will use awk to find the stuff we want.
It sounds like a lot to do; let's see how to proceed.
I will show you the first example.
We are interested in the Y chromosome.
We will esearch in Gene.
Our query is Homo sapiens[orgn] AND Y [chr].
If you want to do another chromosome you would have to edit this query.
If you want to do a generic program that would allow the user to enter a chromosome you could
adjust.
We will not do that today, we'll just hard code in Y.
Filter is a utility that adds another search term to the previous one.
What this does is it looks at current records; alive[prop] in Gene gives you only current
records.
It is filtering the results of the first search by another query term and giving the output
of that; passing that to esummary.
Now we have our docsums for Gene and now we start doing extractions.
We xtract the document summary.
Here is the name variable; we are defining it; -NAME Name.
You have the -block GenomicInfoType.
Here is a trick.
We want to make sure that the chromosome location is actually Y.
So -match "ChrLoc:Y" allows you to text match within the XML for a particular tag.
If you want to make sure you are getting the right thing, you can use -match.
Finally, tab the new line; -tab "\n".
It will put a new line after each record.
After the name, version, ChrStart and ChrStop, new line.
Each gene and it's data will be on a separate line.
Here is how we use the variable NAME.
-element "&NAME".
We use it in the element tag to say I want to have the name positioned in the same output
line with the version, start and stop.
So I say &NAME.
Xtract treats as if it was in the XML in a position.
Let's take a look at how we would do that.
We will run the example.
When I run the example what you see is a bunch of output.
I get the things that I expected.
Here is the gene name.
Here are some interesting ones.
There is the name, here is the accession and the start and stop.
That is good.
It is a start.
But we would like to get the strand.
We need to do a comparison between the start and stop to see if it is on strand one or
strand two.
That is where awk comes into play.
I add a little bit of awk to the end of the script.
I do not know awk that well.
It is pretty readable.
What it is saying is we will use tab characters to separate the input and the output columns.
Each of these different things are separated by tab.
The output will be separated by tab.
Then we do a thing that looks like Perl.
If $3, which is ChrStart, is less than ChrStop it is on the forward strand.
I will print the value of name, $1, ChrAccVer, $2, $3 +1 and $4 +1.
Why is that there?
This is an NCBI geeky thing.
It turns out that those coordinates are zero offset.
Not one offset.
Just be aware of which one you want.
If you want to make it one offset add one to both of them.
If you want to keep them zero offset that is fine.
We are adding the one to the start and stop to make them one offset so one is actually
the first base, as opposed to zero.
And we print the number one.
If it is not the case that three is less than four, then print it the other way around.
We can sort all of this mess by the start and stop coordinates.
This is a standard UNIX sort command.
The tab is our separator.
We will sort by the third column numerically and the fourth column numerically.
Let us see what happens.
Now what I see is I get something that looks different.
I have a number of different genes.
I see the accession in the second column.
I see the start in the third column.
I see the stop and the strand numbers.
We see there are a number of genes on both strands . We have some information.
It has been sorted as it actually occurs on the Y chromosome.
Good, but I do not want the whole Y chromosome.
I am interested in only a partial section of the Y.
How do I do that?
Two possible ways.
First way -- I have two genes to look at.
You can do this with awk.
This is the same script but we have another awk command at the end.
What this does is it tells me I want to look at the region between ASMT, the first gene,
and IL3RA.
I want all the genes between the two.
This is a bit of awkishness, a little different in terms of syntax.
What the syntax is saying is, if you find ASMT, increment a, which is just a variable.
If you find IL3RA, then increment a again.
If a is not zero, print the line.
If a is greater than one, exit.
If it finds one of these genes it will increment a,a will be non-zero and it will print lines.
It will print lines until it gets to the other gene and then a will become 2, which is greater
than one, and it will exit.
You should be able to get the first gene ASMT and print until it finds the other ones.
Does that work?
Let's see.
Now I get a smaller output.
Where I have IL3RA all the way down to ASMT, I have the accession and a start and stop
and strand.
Let's say that what I really wanted to do was not by genes, what I really know is coordinate
positions.
It is a simple modification to the script.
The command at the end becomes another if statement about whether the start and stop
coordinates are greater or less than a particular chromosome coordinate.
If this is a general program you would have to have input by the user or by some other
process.
The point is you can put this in here and logically do the same thing and to prove the
point -- if we ran the script -- this is not with the genes but with the chromosome coordinates.
We get the same thing.
You can do it either way.
What this is trying to show is that by putting the the E-Utilities on the UNIX command line
you have the power of these wonderful utilities.
Sort, uniq, awk.
Or any other UNIX utility that you have.
Last one.
One more example and we'll be done today.
This one is finding proteins produced by a gene.