Dear Dr. Foley,

I can think in several options without create a program/script:

A) Blat-Option:

0) Concatenate all your genomes with Cat into a single fasta file. Modify each of the ID before the concatenation (for example: "sed -r 's/>/>MySpecies_/' file > new_file").

1) Run the Blat with a tabular format output (psl),

2) Using Cut command, select ID, Start and End of the hit

(If you want to get 500 upstream and downstream you can modify the output with

Perl, Python or simply Awk, something like:

"awk '{ print $1"\t"$2-500"\t"$3}' myfile").

3) Use this script to get the sequences from the original fasta (https://github.com/aubombarely/GenoToolBox/blob/master/SeqTools/FastaExtract).

B) Annotation-Option

0) Concatenate the specific regions of your genomes with the hits.

1) Run Maker... if you don't have a lot of sequences, you can run it online specifying as a protein/cDNA set your gene of interest.

2) The output should contain genes in the GFF, you can use Cut to get ID, Start and End and use FastaExtract (see before) to retrieve all the sequences.

The advantage of the second option is that probably you'll have better models.

Other tools that may be interesting to check are the CoGe Blast (http://genomevolution.org/CoGe/CoGeBlast.pl). You have the option of Fasta-View after you get the results where you can retrieve all the fasta results hits.

I hope this helps,

Best,

Aureliano.

Hi!

I am no expert for eukaryotic sequences, but here is a short overview what I would do:

Take a short part of the beginning your protein sequence (50aa should be enough) and use TBLASTN to search for your sequence in all the nucleotide sequences you want. Then download the results as a Hit Table. This gives you the accession number and the start position. Doing the same for the last 50aa gives you the stop position, start and stop can be associated via the GI. The EMBOSS tool package can then be used to extract the regions.

Use Augustus gene prediction tool, it will give you genes and proteins in fasta format.

You can annotate genes using BLAST2GO after you have fasta based file of genes

Thanks for the answers so far. I am not trying to find the gene in my genomes, I want to retrieve and align the genes in complete genomes or contigs in GenBank, based on annotation of the genomes. As reported in Nature in 2009, there are genome projects underway for more than 1,000 vertebrates:

http://www.nature.com/news/2009/091104/full/462021a.html

In theory, if those genomes are to be very useful, it should be an easy task to compare the genes from a given subset of those genomes (mammals, or carnivores, or primates). It should not require each user to do a massive biocomputing project of uploading all of these genomes as raw data to my laptop (with external terabyte drives) and then re-analyzing all the data to find my gene of interest. Each genome project should be, and is, running tools to identify and annotate the genes.

The TBLASTN option is the closest to working for me, but I get thousands of BLAST hits, most of which are from cDNA/mRNA or pseudogenes, and I then have a large job of culling the mniority of hits that are from genes with introns, and devising a way to extract those genes from the whole chromosome or large contig entries in the database.

I was hoping that with so many billions of dollars being spent on acquiring the genome data, the complete genomes of thousands of species, that there was also good money being spent to make it possible for researchers to easily use the data. Asking questions such as "Do all mammals have the same number of hemoglobin genes?" should be rather easy to answer.

Non scripted options could be just using BiosMART in Ensemble. Their menus are flexible enough to customise your selection.

If you REALLY need programmatic access (so you can't simply download a chunk of sequences to access from your disk) then various sequences databases, including Enseble again, provide PERL or other interfaces to query sequences, which you can then pipe into your particularl analysis.

GenBank provides all the tools needed to get one gene at a time, so there might be something in BioPERL that would automate the process. The gene is on chromosome 19 in human, chromosome 10 in mouse.

After selecting the gene region (plus 1,000 bases on each side of the gene to get some flanking region), and updating the view, I download the GenBank view in order to preserve the annotations. If I download FASTA format, I loose information about where the exons are, who the authors were, etc.

Converting GenBank to FASTA later is trivial, and can be done after concatenating hundreds of the GenBank format files together. Doing this "by hand" using GenBank takes about 2 minutes per species. It is not a huge job really, but if I want thousands of genes it gets a little tedious.

From my experiences not all mammals will have same gene order or genes on the chromosomes, duplications happens before and after separation of species that will lead into new duplicated genes in chromosomes. Example:- human has 13 SerpinA genes where as mouse >25 SerpinA genes.

Leaving aside the technicalities of accessing the data, I think Abhishek comment is very relevant. Different species have a set of common genes (orthologs), but also many specific ones that you won't find in a second species. On the other end, duplications can give you multiple orthologs. Still I would say this is not a real problem for a computer. More serious issue is to find the name of the gene of interest in the other species. Trivial as it sounds, this could be a big problem since naming across species (or even databases, authors, etc) can be very variable. In this respect, I would suggest to start your work pulling as many synonyms of your gene as possible and use them all to identify orthologs across all genomes. (Uniprot can be useful for this).

Great comments. The problem of gene duplications and variations in naming of genes is one of the reasons I like the EF-2 gene for molecular phylogeny. Although all organisms from bacteria to mammals (or at least all that I have looked at) have a second gene encoding a similar "EF-2 like" protein, none so far has been found to have two copies of the EF-2 gene itself (there are some processed pseudogenes though). EF-2 is a "housekeeping" gene, absolutely required for life, so gene loss is not an issue. It is not clear why gene duplication seems impossible, but could have to do with tight regulation of this factor in the control of rates of protein synthesis.

I don't really need the complete gene from hundreds of vertebrates. Just from looking at an alignment of a few of them so far, spread through the mammal radiation, I can see that alignment of the introns is going to be useless for phylogenetic purposes. Intron length seems highly variable. I see almost no detectable homology/similarity between the introns in the mouse gene and the introns in the human gene. This is a bit surprising to me, given the very high degree of similarity in the exons, and given the high degree of similarity of the mitochondrial genomes. Within the primates, most of the changes in intron sequence seem to be insertion/deletion events rather than point mutations. Phylogenetic analysis calculations are not well suited to analyzing in/del events, because there are so many different types. Insertion of an ALU repeat for example is different than other events. Analysis of non-coding DNA (introns or intergenic regions for example) is of interest though, because there should be much less selection pressure confounding the analysis.

As I understand it, you want to pull the full EF-2 gene with introns from a group of mitochondrial genomes in varying stages of annotation and assembly. I don't think there is an easy way to do this.

The easiest way to do this, if enough of the genomes are annotated and have full gene sequences to meet your needs, would be to create blast DB's with the gene sequences, and write a script to search them with a reference EF-2 gene and collect the top hits. I would then take the name of the gene which was the top hit and pull it out. (or collect the ID's of these genes to get via biomart or similar system, if such a thing exists for your genome of interest).

If there aren't enough genomes annotated to this level, I would use an EST alignment program, such as EST2genome (there may be newer, better ones out there), and align a reference EF-2 coding sequence. Then I would take the start/end coordinates for the mapped EST and pull everything in between out. I might implement a check to make sure the start/end sequence of the reference EF-2 CDS mapped, to make sure chunks aren't missing.

You might be able to use an graphical genome annotation tool to do this, but then you would have to do it manually- scripting up something would allow you to do it on over 200 genomes with a few hours worth of work, a graphical tool would let you do a few by hand without knowing how to script.

Mitochondria are derived from bacteria. They are essentially prokayotic, and there are no introns in mitochondrial genes. The elongation factor 2 gene is nuclear, it resides on human chromosome 19, pig chromosome 2, and so on. A major goal of sequencing the complete genomes of thousands of vertebrates, is to be able to do some comparative genomics and understand evolution of vertebrates. So I was guessing that in addition to sequencing these genomes and annotating them as described here:

http://www.nature.com/news/2009/091104/full/462021a.html

they would also be developing tools which would make it simple to do some comparative genomics. The ENSEMBL tools do a lot of good stuff, but don't quite give me what I need. With ENSEMBL, it is easy to find that there are 152 genes annotated as elongation factor 2 gene, but I am not finding a way to download the 152 genes all at once, aligned or not.

ENSEMBL seems to be telling me it could not produce an alignment.

I was unsure what you actually needed, given all the background information, which is why I restated what it looked like you were looking for. The methods I gave will work with or without introns.

If you are specifically looking to download from ensembl, you can retrieve orthologs via biomart, though the interface and what is available is quite limiting. It may be possible to connect to their mysql DB and retrieve sequences annotated as "elongation factor 2".

If you are talking generally... no, there is really a lack of good tools (and simply how data is organized and stored- having the same gene have different ID's for ensembl, entrez, flybase, etc is a headache for programmers) to analyze data on the genome level, the field is really just getting established now that it can be done for thousands of dollars rather than hundreds of thousands or millions.

I would expect that in the years to come, tools will improve, as tools for transcritomic analysis, etc, have. But it is still a bit of putting pieces together ones self for many analyses.

The ENSEMBL site definitely seems to be very close to giving me what I want. I am sure that most of my problem is "operator error", I do not quite know what buttons to push there yet. But it is clear that they have made an alignment of dozens of eukaryotic Elongation Factor 2 genes, and it is a bit easier for me to pull each of them out individually there, than from GenBank. So far it is not giving me the alignment, but if I can get the genes out, I can align them myself.

So, after getting 21 of the EF-2 genes clipped out of the mammalian genomes, it turns out that aligning them is not trivial. The exons are very highly conserved, but the introns have no detectable similarity left. Within the (human chimpanzee gorilla orangutan) clade the introns are of course very similar, but aligning primate to other mammals is going to be essentially useless.

This is not entirely surprising, but it is interesting. The gene order, surrounding the EF-2 gene, was conserved. Almost all mammals had a "Death-Associated Protein Kinase 3 gene" immediately to the 3' (downstream) side of the EF-2 gene.

The protein coding regions of genes are very highly conserved, while the introns and intergenic regions have diversified almost beyond recognition. I have attached a photo showing Bovine (mammal) vs chicken (aves) protein comparison, as an example. Negative selection to conserve the proteins is a bit more extreme than I had thought it might be.

If you're still interested in this a nice approach would be to use the fastq files (ALWAYS use the reads since at least you know the error), align the the reference genes, convert to species -specific genes and use these as a new reference (separately for each species). Repeat alignments until you have all your data. If you need scripts for all this let me know.