View source for Utilities for comparative genomics

== State of the art: Comparative Genomics ==

=== Populating a feature stack ===

Populating a feature stack (eg of a coding exon) is the central task in comparative genomics. Feature stacks collect all information available from all available species on a suitable topic. Ideally the topic, for example evolutionary history of a coding exon, has width in base pairs a fraction of characteristic width of available experimental data. That is, a coding gene will have single exons wholly contained within individual trace reads and multiple exons within gapless GenBank wgs contigs. 

On the other hand, a segmental duplication of 10 genes, say the one on chr20q/chr8, is very unsuited to a feature stack today because only a couple mammals have an adequate assembly over the intrinsic span of the feature. 

Thus it is timely to construct feature stacks for coding exons and the like but premature to consider the comparative genomics of longer features. Consequently all research today that can fully exploit the power of comparative genomics is restricted to computable feature stacks. In other words, if someone can precompute all possible contemporary feature stacks, in effect they have written all possible research papers. This makes them the masters of the comparative genomics universe.

Brian Raney recently has computed this for a genome-based multiz alignment of 28 species, the UCSC 28way [http://genome-test.cse.ucsc.edu/cgi-bin/hgTrackUi?hgsid=1506267&c=chr22&g=multiz28way conservation track] for both [http://hgwdev.cse.ucsc.edu/%7Ebraney/28way.exonNuc/ nucleotides] and [http://hgwdev.cse.ucsc.edu/%7Ebraney/28way.exons/ proteins].

How many feature stacks are there, how difficult is to compute, store, and query them, and what associated precomputed products go with them? Suppose there are 190,000 coding exon and 3 non-coding features per 20,000 gene for 240,000 features of width 500 bp and needed depth of 50 species. That would fit handily within an excel spreadsheet. So the number, storage, and query are non-issues. 

Note thought that the 28way alignment is far from perfect with respect to mis-populating exons with pseudogenes and misalignments and infill completeness, not to mention that 50 vertebrate genomes are actually available in some form. Thus considerable manual curation and infill from the trace archives (long and short) and GenBank database divisions such as nr, est_others, and wgs. In May 2008, typically 43-44 orthologs for a given exon can be located

Here's an example of a feature stack and the implicit paper that practically writes itself:

 <font face="Courier">Comparative Genomics of DRY motifs in exon 3 of RGR Opsins:</font>
 <font face="Courier">1 RWPYGSDGCQAHGFQGFVTALASICSSAAIAW<font color="magenta">G</font>RYHHYCT 1 human</font>
 <font face="Courier">1 RWPYGSDGCQAHGFQGFVTALASICSSAAIAW<font  color="magenta">G</font>RYHHYCT 1 macaque</font>
 <font face="Courier">1 RWPYGSGGCQAHGFQGFTTALASICGSAAIAW<font  color="magenta">G</font>RYHHYCT 1 lemur</font>
 <font face="Courier">1 RWPHGSEGCQVHGFQGFATALASICGSAAVAW<font  color="magenta">G</font>RYHHYCT 1 mouse</font>
 <font face="Courier">1 RWPYGSDGCQAHGFQGFATALASICGSAAIAW<font  color="magenta">G</font>RYHHYCT 1 rabbit</font>
 <font face="Courier">1 RWPYGSDGCQAHGFQGFVTALASICSSAAIAW<font  color="magenta">G</font>RYHHYCT 1 horse</font>
 <font face="Courier">1 RWPYGPDGCQAHGFQGFATALASICSSAALAW<font  color="magenta">G</font>RYHHYCT 1 dog</font>
 <font face="Courier">1 RWPYGSGGCQAHGFQGFAAALASICGSAAVAW<font  color="magenta">G</font>RYHHYCT 1 bat</font>
 <font face="Courier">1 RWPYGSDGCQAHGFQGFVTALASICSCAAIAW<font  color="blue">E</font>RYHHYCT 1 elephant</font>
 <font face="Courier">1 HWPYGSGGCQAHGFQGFTVALASICSCAAIAW<font  color="blue">E</font>RYHHYCT 1 tenrec</font>
 <font face="Courier">1 RWPHGSDSCQAHSFQGFATALASISSSAAIAW<font  color="blue">E</font>RYRHHCT 1 sloth</font>
 <font face="Courier">1 RWPYGSGGCQAHGFQGFVTALASISSSAAIAW<font  color="blue">E</font>RCHRHCI 1 armadillo</font>
 <font face="Courier">1 HWPYGAEGCRLHGFQGFATALASISLSAAIGW<font  color="#00CC66">D</font>RYLRHCS 1 platypus</font>
 <font face="Courier">1 YWPYGSDGCQIHGFHGFLTALTSISSAAAVAW<font  color="#00CC66">D</font>RHHQYCT 1 lizard</font>
 <font face="Courier">1 YWPYGSEGCQIHGFQGFLTALASISSSAAVAW<font  color="#00CC66">D</font>RYHHYCT 1 chicken</font>
 <font face="Courier">1 YWPYGSEGCQIHGFQGFVAALSSIGSCAAIAW<font  color="#00CC66">D</font>RYHQYCT 1 frog</font>
 <font face="Courier">1 YWPYGSDGCQTHGFQGFVTALASIHFIAAIAW<font  color="#00CC66">D</font>RYHQYCT 1 stickleback</font>
 <font face="Courier">1 YWPYGSDGCQTHGFQGFVTALASIHFVAAIAW<font  color="#00CC66">D</font>RYHQYCT 1 fugu</font>
 <font face="Courier">1 YWPYGSEGCQTHGFHGFLTALASIHFIAAIAW<font  color="#00CC66">D</font>RYHQYCT 1 medaka</font>
 <font face="Courier">1 YWPYGSEGCQTHGFHGFLMALASINACAAIAW<font  color="#00CC66">D</font>RYHQNCS 1 elephantshark</font>
 <font face="Courier">1 EWPFGSIGCQLDAFIGMAPTFISIAGAALVAK<font  color="#00CC66">D</font>KYYRICK 1 tunicate </font>


Now half the data needed for a full feature stack is not available at any genome browser but rather at GenBank etc. Whether that data is retrieved by automated pipeline or by one-off queries, there's a universal roadblock in that blast output is terribly formatted for comparative genomics purposes:
 
* the blast algorithm is unaware of splice donors and acceptors and ignores user parsing
* blocks are provided in score order instead of natural query exon order
* blocks have faux match extensions at the ends (exon dribble) complicating clean ortholog extraction
* blocks can contain multiple exons when introns are short leading to nonsense in tblastn
* line width of 60 is way too short, causing string searches to fail because of carriage returns
* gaps necessitate dashes causing string searches to fail
* the match genus and species are inconveniently provided

Fixing Blast output in order to populate a comparative genomics feature stack: No utility web tool exists anywhere on the internet that can repair blast output. Consequently one would make a great addition. The needed algorithm would largely carry over to any pipeline program to populate feature stacks genomewide. The example below shows what needs to happen, as limited to one species. In this case, a single tblastn query to GenBank wgs adds 19 mammals to the feature stack.

<pre>Here is the original query, 3 exons of bovine rhodpsin marked up for reading frame

>RHO1_bosTau Bos taurus (cow)  
2 GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIMGVAFTWVMALACAAPPLVGWSR 2
1 YIPEGMQCSCGIDYYTPHEETNNESFVIYMFVVHFIIPLIVIFFCYGQLVFTVKE 0
0 AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGSDFGPIFMTIPAFFAKTSAVYNPVIYIMMNKQ 0
</pre>
<pre>Here is the raw output match to marsupial

>gb|AAFR03021222.1| Monodelphis domestica cont3.021221, whole genome shotgun sequence Length=94985

 Score =  149 bits (377),  Expect = 9e-34, Method: Compositional matrix adjust.
 Identities = 71/74 (95%), Positives = 74/74 (100%), Gaps = 0/74 (0%)  Frame = +1

Query  119    AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGSDFGPIFMTIPAFFAKTS  178
              AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGS+FGPIFMTIPAFFAK+S
Sbjct  21139  AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGSNFGPIFMTIPAFFAKSS  21318

Query  179    AVYNPVIYIMMNKQ  192
              +VYNPVIYIMMNKQ
Sbjct  21319  SVYNPVIYIMMNKQ  21360 

 Score =  117 bits (294),  Expect = 5e-24, Method: Compositional matrix adjust.
 Identities = 54/59 (91%), Positives = 57/59 (96%), Gaps = 0/59 (0%) Frame = +3

Query  1      GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIMGVAFTWVMALACAAPPLVGWSRYI  59
              GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAI+GVAFTWVMALACA PPL+GWSR +
Sbjct  18648  GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIIGVAFTWVMALACAFPPLIGWSRLV  18824

 Score =  107 bits (268),  Expect = 4e-21, Method: Compositional matrix adjust.
 Identities = 51/59 (86%), Positives = 52/59 (88%), Gaps = 0/59 (0%) Frame = +2

Query  55     WSRYIPEGMQCSCGIDYYTPHEETNNESFVIYMFVVHFIIPLIVIFFCYGQLVFTVKEA  113
              + RYIPEGMQCSCGIDYYT   E NNESFVIYMFVVHF IPLIVIFFCYGQLVFTVKE 
Sbjct  20546  FCRYIPEGMQCSCGIDYYTLKPEVNNESFVIYMFVVHFTIPLIVIFFCYGQLVFTVKEV  20722
</pre>

<pre>Stage 1: order the output blocks according to input exon order and fix line width:

Query  1      GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIMGVAFTWVMALACAAPPLVGWSRYI  59
              GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAI+GVAFTWVMALACA PPL+GWSR +
Sbjct  18648  GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIIGVAFTWVMALACAFPPLIGWSRLV  18824

Query  55     WSRYIPEGMQCSCGIDYYTPHEETNNESFVIYMFVVHFIIPLIVIFFCYGQLVFTVKEA  113
              + RYIPEGMQCSCGIDYYT   E NNESFVIYMFVVHF IPLIVIFFCYGQLVFTVKE 
Sbjct  20546  FCRYIPEGMQCSCGIDYYTLKPEVNNESFVIYMFVVHFTIPLIVIFFCYGQLVFTVKEV  20722

Query  119    AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGSDFGPIFMTIPAFFAKTSAVYNPVIYIMMNKQ  192
              AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGS+FGPIFMTIPAFFAK+S+VYNPVIYIMMNKQ
Sbjct  21139  AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGSNFGPIFMTIPAFFAKSSSVYNPVIYIMMNKQ  21360
</pre>
 <font face="Courier" color="magenta>Stage 2: Find the correct input exons</font>
 <font face="Courier">Query  1      <font  color="magenta">GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIMGVAFTWVMALACAAPPLVGWSR</font>YI  59</font>
 <font face="Courier">              GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAI+GVAFTWVMALACA PPL+GWSR +</font>
 <font face="Courier">Sbjct  18648  GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIIGVAFTWVMALACAFPPLIGWSRLV  18824</font>
 <font face="Courier">Query  55     WSR<font  color="magenta">YIPEGMQCSCGIDYYTPHEETNNESFVIYMFVVHFIIPLIVIFFCYGQLVFTVKE</font>A  113</font>
 <font face="Courier">              + RYIPEGMQCSCGIDYYT   E NNESFVIYMFVVHF IPLIVIFFCYGQLVFTVKE </font>
 <font face="Courier">Sbjct  20546  FCRYIPEGMQCSCGIDYYTLKPEVNNESFVIYMFVVHFTIPLIVIFFCYGQLVFTVKEV  20722</font>
 <font face="Courier">Query  119    <font color="magenta">AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGSDFGPIFMTIPAFFAKTSAVYNPVIYIMMNKQ</font>  192</font>
 <font face="Courier">              AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGS+FGPIFMTIPAFFAK+S+VYNPVIYIMMNKQ</font>
 <font face="Courier">Sbjct  21139  AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGSNFGPIFMTIPAFFAKSSSVYNPVIYIMMNKQ  21360</font>

 <font face="Courier"  color="blue>Stage 3: Transfer the match to marsupial</font>
 <font face="Courier">Query  1      GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIMGVAFTWVMALACAAPPLVGWSRYI  59</font>
 <font face="Courier">              GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAI+GVAFTWVMALACA PPL+GWSR +</font>
 <font face="Courier">Sbjct  18648  <font  color="blue">GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIIGVAFTWVMALACAFPPLIGWSR</font>LV  18824</font>
 <font face="Courier">Query  55     WSRYIPEGMQCSCGIDYYTPHEETNNESFVIYMFVVHFIIPLIVIFFCYGQLVFTVKEA  113</font>
 <font face="Courier">              + RYIPEGMQCSCGIDYYT   E NNESFVIYMFVVHF IPLIVIFFCYGQLVFTVKE </font>
 <font face="Courier">Sbjct  20546  FCR<font  color="blue">YIPEGMQCSCGIDYYTLKPEVNNESFVIYMFVVHFTIPLIVIFFCYGQLVFTVKE</font>V  20722</font>
 <font face="Courier">Query  119    AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGSDFGPIFMTIPAFFAKTSAVYNPVIYIMMNKQ  192</font>
 <font face="Courier">              AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGS+FGPIFMTIPAFFAK+S+VYNPVIYIMMNKQ</font>
 <font face="Courier">Sbjct  21139  <font  color="blue">AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGSNFGPIFMTIPAFFAKSSSVYNPVIYIMMNKQ</font>  21360</font>
 
<pre>Stage 4: Export the desired output with acquired header

>RHO1_monDom Monodelphis domesticus (opossum)
2 GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIIGVAFTWVMALACAFPPLIGWSR 2
1 YIPEGMQCSCGIDYYTLKPEVNNESFVIYMFVVHFTIPLIVIFFCYGQLVFTVKE 0
0 AAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGSNFGPIFMTIPAFFAKSSSVYNPVIYIMMNKQ 0
</pre>

<pre>Stage 5: Fill in the respective feature stacks with marsupial exons (only first is shown):
>RHO1_homSap GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIMGVAFTWVMALACAAPPLAGWSR
>RHO1_bosTau GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIMGVAFTWVMALACAAPPLVGWSR
>RHO1_monDom GEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIIGVAFTWVMALACAFPPLIGWSR
>RHO1_ornAna GEIALWSLVVLAIERYIVVCKPMSNFRFGENHAIMGVAFTWIMALACALPPLVGWSR
>RHO1_galGal GEIALWSLVVLAVERYVVVCKPMSNFRFGENHAIMGVAFSWIMAMACAAPPLFGWSR
>RHO1_anoCar GEMGLWSLVVLAVERYVVICKPMSNFRFGETHALIGVSCTWIMALACAGPPLLGWSR
>RHO1_xenTro GEMALWSLVVLAIERYVVVCKPMANFRFGENHAIMGVVFTWIMALSCAAPPLFGWSR
>RHO1_neoFor GIIALWCLVVLAIERYIVVCKPISNFRFGENHAIMGVVFTWIMALACAGPPLFGWSR
>RHO1_latCha GQVALWALVVLAIERYVVVCKPMSNFRFGENHAIMGVIFTWIMALSCAVPPLFGWSR
>RHO1_takRub GEIALWSLVVLAVERYIVVCKPMTNFRFGEKHAIAGLVFTWIMALTCATPPLLGWSR
>RHO1_leuEri GEVGLWCLVVLAIERYMVVCKPMANFRFGSQHAIIGVVFTWIMALSCAGPPLVGWSR
>RHO1_calMil GEIGLWSLVVLAIERYVVVCKPMSNFRFGTNHAIMGVAFTWVMALACAVPPLMGWSR
>RHO1_petMar DEMSLWSLVVLAIERYIVICKPMGNFRFGSTHAYMGVAFTWFMALSCAAPPLVGWSR
</pre>

An alternative approach -- simply collecting the start and stop numbering in the Sbjct (match) line -- might work better. Another flaw in blast is that no flanking material is provided in the event of the first residue or two in an exon not matching. In the above example, blast actually dropped a perfect matches to the initial hexapeptide AAAQQQ even though the simple sequence filter was off. 

Thus the processing algorithm could gather these spanning numbers. Entrez retrieval allows limitation to these in conjunction with the accesssion numbers. That would allow for orderly recovery of adequately padded exonic regions which could then be concatenated for translation in a consistent frame with validation of intron position and phase. Alignment could then be performed in an altered version of blast more sympathetic to the goals here. This would better address the special problem of split codon reconstruction in the 12 and 21 overhang situations (where the completed codon does not appear in extended translation of either exon.

=== Species commonly available ===

To conveniently collect exons from a given gene, place the list below in one column of a spreadsheet (after search-replace has put in the name of the gene being specifically investigated). The next 3 columns can contain the splice acceptor phase, the exon sequences themselves, and the donor phase. After all the data is collected, it is converted to standard exonic fasta format by replacing tabs with returns. 

The numbers to the left allow the spreadsheet to be sorted by  output order in the UCSC 28way alignment, phylogenetic quasi-ordering with respect to the vertebrate tree (human arbitrarily first), or alphabetically by genus-species. Those three orders are provided as pre-sorts below.


<pre>
>10.10.gene_homSap Homo sapiens (human)
>11.11.gene_panTro Pan troglodytes (chimp)
>99.12.gene_gorGor Gorilla gorilla (gorilla)
>99.13.gene_ponPyg Pongo pygmaeus (orang_sumatran)
>99.14.gene_nomLeu Nomascus leucogenys (gibbon)
>12.15.gene_macMul Macaca mulatta (rhesus)
>99.16.gene_papAnu Papio anubis (baboon)
>99.17.gene_papHam Papio hamadryas (baboon)
>99.18.gene_calJac Callithrix jacchus (marmoset)
>99.19.gene_tarSyr Tarsius syrichta (tarsier)
>13.20.gene_otoGar Otolemur garnettii (bushbaby)
>99.21.gene_micMur Microcebus murinus (mouse_lemur)
>99.22.gene_cynVol Cynocephalus volans (flying_lemur)
>14.23.gene_tupBel Tupaia belangeri (tree_shrew)
>15.24.gene_musMus Mus musculus (mouse)
>16.25.gene_ratNor Rattus norvegicus (rat)
>17.26.gene_cavPor Cavia porcellus (guinea_pig)
>99.27.gene_speTri Spermophilus tridecemlineatus (squirrel)
>99.28.gene_dipOrd Dipodomys ordii (kangaroo_rat)
>18.29.gene_oryCun Oryctolagus cuniculus (rabbit)
>99.30.gene_ochPri Ochotona princeps (pika)
>21.31.gene_canFam Canis familiaris (dog)
>22.32.gene_felCat Felis catus (cat)
>24.33.gene_bosTau Bos taurus (cow)
>99.34.gene_oviAri Ovis aries (sheep)
>99.35.gene_susScr Sus scrofa (pig)
>23.36.gene_equCab Equus caballus (horse)
>99.37.gene_myoLuc Myotis lucifugus (microbat)
>99.38.gene_pteVam Pteropus vampyrus (macrobat)
>99.39.gene_turTru Tursiops truncatus (dolphin)
>99.40.gene_susScr Sus scrofa (pig)
>99.41.gene_vicVic Vicugna vicugna (vicugna)
>19.42.gene_eriEur Erinaceus europaeus (hedgehog)
>20.43.gene_sorAra Sorex araneus (shrew)
>99.44.gene_borAnc Boreoeuthere ancestralis (ancestral)
>25.45.gene_dasNov Dasypus novemcinctus (armadillo)
>99.46.gene_choHof Choloepus hoffmanni (sloth)
>26.47.gene_loxAfr Loxodonta africana (elephant)
>99.48.gene_proCap Procavia capensis (hyrax)
>99.49.gene_echTel Echinops telfairi (tenrec)
>27.50.gene_monDom Monodelphis domestica (opossum)
>99.51.gene_macEug Macropus eugenii (wallaby)
>99.52.gene_triVul Trichosurus vulpecula (possum)
>28.53.gene_ornAna Ornithorhynchus anatinus (platypus)
>99.54.gene_tacAcu Tachyglossus aculeatus (echidna)
>30.55.gene_galGal Gallus gallus (chicken)
>99.56.gene_taeGut Taeniopygia guttata (finch)
>29.57.gene_anoCar Anolis carolinensis (lizard)
>31.58.gene_xenTro Xenopus tropicalis (frog)
>99.59.gene_xenTro Xenopus laevis (frog)
>99.60.gene_neoFor Neoceratodus forsteri (lungfish)
>32.61.gene_danRer Danio rerio (zebrafish)
>33.62.gene_tetNig Tetraodon nigroviridis (pufferfish)
>34.63.gene_takRub Takifugu rubripes (fugu)
>35.64.gene_gasAcu Gasterosteus  aculeatus (stickleback)
>36.65.gene_oryLap Oryzias latipes (medaka)
>99.66.gene_ictPun Ictalurus punctatus (fish)
>99.67.gene_oncMyk Oncorhynchus mykiss (trout)
>99.68.gene_funHet Fundulus heteroclitis (flounder)
>99.69.gene_calMil Callorhinchus milii (elephantfish)
>99.70.gene_squAca Squalus acanthias (spiny dogfish)
>99.71.gene_petMar Petromyzon marinus (lamprey)
>99.72.gene_braFlo Branchiostoma floridae (amphioxus)

>10.10.gene_homSap Homo sapiens (human)
>11.11.gene_panTro Pan troglodytes (chimp)
>12.15.gene_macMul Macaca mulatta (rhesus)
>13.20.gene_otoGar Otolemur garnettii (bushbaby)
>14.23.gene_tupBel Tupaia belangeri (tree_shrew)
>15.24.gene_musMus Mus musculus (mouse)
>16.25.gene_ratNor Rattus norvegicus (rat)
>17.26.gene_cavPor Cavia porcellus (guinea_pig)
>18.29.gene_oryCun Oryctolagus cuniculus (rabbit)
>19.42.gene_eriEur Erinaceus europaeus (hedgehog)
>20.43.gene_sorAra Sorex araneus (shrew)
>21.31.gene_canFam Canis familiaris (dog)
>22.32.gene_felCat Felis catus (cat)
>23.36.gene_equCab Equus caballus (horse)
>24.33.gene_bosTau Bos taurus (cow)
>25.45.gene_dasNov Dasypus novemcinctus (armadillo)
>26.47.gene_loxAfr Loxodonta africana (elephant)
>27.50.gene_monDom Monodelphis domestica (opossum)
>28.53.gene_ornAna Ornithorhynchus anatinus (platypus)
>29.57.gene_anoCar Anolis carolinensis (lizard)
>30.55.gene_galGal Gallus gallus (chicken)
>31.58.gene_xenTro Xenopus tropicalis (frog)
>32.61.gene_danRer Danio rerio (zebrafish)
>33.62.gene_tetNig Tetraodon nigroviridis (pufferfish)
>34.63.gene_takRub Takifugu rubripes (fugu)
>35.64.gene_gasAcu Gasterosteus  aculeatus (stickleback)
>36.65.gene_oryLap Oryzias latipes (medaka)
>99.12.gene_gorGor Gorilla gorilla (gorilla)
>99.13.gene_ponPyg Pongo pygmaeus (orang_sumatran)
>99.14.gene_nomLeu Nomascus leucogenys (gibbon)
>99.16.gene_papAnu Papio anubis (baboon)
>99.17.gene_papHam Papio hamadryas (baboon)
>99.18.gene_calJac Callithrix jacchus (marmoset)
>99.19.gene_tarSyr Tarsius syrichta (tarsier)
>99.21.gene_micMur Microcebus murinus (mouse_lemur)
>99.22.gene_cynVol Cynocephalus volans (flying_lemur)
>99.27.gene_speTri Spermophilus tridecemlineatus (squirrel)
>99.28.gene_dipOrd Dipodomys ordii (kangaroo_rat)
>99.30.gene_ochPri Ochotona princeps (pika)
>99.34.gene_oviAri Ovis aries (sheep)
>99.35.gene_susScr Sus scrofa (pig)
>99.37.gene_myoLuc Myotis lucifugus (microbat)
>99.38.gene_pteVam Pteropus vampyrus (macrobat)
>99.39.gene_turTru Tursiops truncatus (dolphin)
>99.40.gene_susScr Sus scrofa (pig)
>99.41.gene_vicVic Vicugna vicugna (vicugna)
>99.44.gene_borAnc Boreoeuthere ancestralis (ancestral)
>99.46.gene_choHof Choloepus hoffmanni (sloth)
>99.48.gene_proCap Procavia capensis (hyrax)
>99.49.gene_echTel Echinops telfairi (tenrec)
>99.51.gene_macEug Macropus eugenii (wallaby)
>99.52.gene_triVul Trichosurus vulpecula (possum)
>99.54.gene_tacAcu Tachyglossus aculeatus (echidna)
>99.56.gene_taeGut Taeniopygia guttata (finch)
>99.59.gene_xenTro Xenopus laevis (frog)
>99.60.gene_neoFor Neoceratodus forsteri (lungfish)
>99.66.gene_ictPun Ictalurus punctatus (fish)
>99.67.gene_oncMyk Oncorhynchus mykiss (trout)
>99.68.gene_funHet Fundulus heteroclitis (flounder)
>99.69.gene_calMil Callorhinchus milii (elephantfish)
>99.70.gene_squAca Squalus acanthias (spiny dogfish)
>99.71.gene_petMar Petromyzon marinus (lamprey)
>99.72.gene_braFlo Branchiostoma floridae (amphioxus)

>29.57.gene_anoCar Anolis carolinensis (lizard)
>99.44.gene_borAnc Boreoeuthere ancestralis (ancestral)
>24.33.gene_bosTau Bos taurus (cow)
>99.72.gene_braFlo Branchiostoma floridae (amphioxus)
>99.18.gene_calJac Callithrix jacchus (marmoset)
>99.69.gene_calMil Callorhinchus milii (elephantfish)
>21.31.gene_canFam Canis familiaris (dog)
>17.26.gene_cavPor Cavia porcellus (guinea_pig)
>99.46.gene_choHof Choloepus hoffmanni (sloth)
>99.22.gene_cynVol Cynocephalus volans (flying_lemur)
>32.61.gene_danRer Danio rerio (zebrafish)
>25.45.gene_dasNov Dasypus novemcinctus (armadillo)
>99.28.gene_dipOrd Dipodomys ordii (kangaroo_rat)
>99.49.gene_echTel Echinops telfairi (tenrec)
>23.36.gene_equCab Equus caballus (horse)
>19.42.gene_eriEur Erinaceus europaeus (hedgehog)
>22.32.gene_felCat Felis catus (cat)
>99.68.gene_funHet Fundulus heteroclitis (flounder)
>30.55.gene_galGal Gallus gallus (chicken)
>35.64.gene_gasAcu Gasterosteus  aculeatus (stickleback)
>99.12.gene_gorGor Gorilla gorilla (gorilla)
>10.10.gene_homSap Homo sapiens (human)
>99.66.gene_ictPun Ictalurus punctatus (fish)
>26.47.gene_loxAfr Loxodonta africana (elephant)
>99.51.gene_macEug Macropus eugenii (wallaby)
>12.15.gene_macMul Macaca mulatta (rhesus)
>99.21.gene_micMur Microcebus murinus (mouse_lemur)
>27.50.gene_monDom Monodelphis domestica (opossum)
>15.24.gene_musMus Mus musculus (mouse)
>99.37.gene_myoLuc Myotis lucifugus (microbat)
>99.60.gene_neoFor Neoceratodus forsteri (lungfish)
>99.14.gene_nomLeu Nomascus leucogenys (gibbon)
>99.30.gene_ochPri Ochotona princeps (pika)
>99.67.gene_oncMyk Oncorhynchus mykiss (trout)
>28.53.gene_ornAna Ornithorhynchus anatinus (platypus)
>18.29.gene_oryCun Oryctolagus cuniculus (rabbit)
>36.65.gene_oryLap Oryzias latipes (medaka)
>13.20.gene_otoGar Otolemur garnettii (bushbaby)
>99.34.gene_oviAri Ovis aries (sheep)
>11.11.gene_panTro Pan troglodytes (chimp)
>99.16.gene_papAnu Papio anubis (baboon)
>99.17.gene_papHam Papio hamadryas (baboon)
>99.71.gene_petMar Petromyzon marinus (lamprey)
>99.13.gene_ponPyg Pongo pygmaeus (orang_sumatran)
>99.48.gene_proCap Procavia capensis (hyrax)
>99.38.gene_pteVam Pteropus vampyrus (macrobat)
>16.25.gene_ratNor Rattus norvegicus (rat)
>20.43.gene_sorAra Sorex araneus (shrew)
>99.27.gene_speTri Spermophilus tridecemlineatus (squirrel)
>99.70.gene_squAca Squalus acanthias (spiny dogfish)
>99.35.gene_susScr Sus scrofa (pig)
>99.40.gene_susScr Sus scrofa (pig)
>99.54.gene_tacAcu Tachyglossus aculeatus (echidna)
>99.56.gene_taeGut Taeniopygia guttata (finch)
>34.63.gene_takRub Takifugu rubripes (fugu)
>99.19.gene_tarSyr Tarsius syrichta (tarsier)
>33.62.gene_tetNig Tetraodon nigroviridis (pufferfish)
>99.52.gene_triVul Trichosurus vulpecula (possum)
>14.23.gene_tupBel Tupaia belangeri (tree_shrew)
>99.39.gene_turTru Tursiops truncatus (dolphin)
>99.41.gene_vicVic Vicugna vicugna (vicugna)
>99.59.gene_xenTro Xenopus laevis (frog)
>31.58.gene_xenTro Xenopus tropicalis (frog)
</pre>

== Obtaining and Analyzing Sequences from 454 Transcript Runs ==

Sanger genomic sequencing basically shut down at the major sequencing centers in early April 2008. Sequencing of vertebrates has shifted over to faster and cheaper new technologies, predominantly 454. These reads are deposited in a user-impractical form at the NCBI [tp://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi Short Read Archives.]

Notice that gene models are highly mediocre outside species such as human and mouse, yet many mammalian species have essentially no available transcripts. This very much detracts from the values of individual species browsers. Other species such as Tachyglossus have no genome project underway but several hundred thousand transcript reads are available. Thus it would be very advantageous to bring the voluminous 454 data into play. 

On 22 April 08, the following large scale 454 transcript and genome projects were [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?cmd=tracking&f=&m=main&s=tracking available:]
<pre>
homSap SRA000297 Homo sapiens             Transcriptome Study 454 GS FLX       WUGSC
calJac SRA000144 Callithrix jacchus       Transcriptome Study                  WUGSC 
tupBel SRA000295 Tupaia belangeri         Transcriptome Study 454 GS FLX       WUGSC
susScr SRA000267 Sus scrofa               Transcriptome study 454 GS 20        AU-DJ
ornAna SRA000294 Ornithorhynchus anatinus Transcriptome Study 454 GS FLX       WUGSC
galGal SRA000238 Gallus gallus            Transcriptome Study 454 GS FLX       WUGSC
tacAcu SRA000241 Tachyglossus aculeatus   Transcriptome Study 454 GS FLX       WUGSC

turTru SRA000063 Tursiops truncatus       Whole Genome Sequencing Project      BCM
mamPri ......... Mammuthus primigenius    Whole Genome Sequencing Project      PSU
monDom SRA000162 Monodelphis domestica    Whole Genome Sequencing              WUGSC
macEug SRA000041 Macropus eugenii         Whole Genome Sequencing Project      BCM
xenTro SRA000296 Xenopus tropicalis       Genomic Sequencing 454 GS FLX        WUGSC
</pre>

For the average user, there are 4 issues that impede effective access to all this great new data:

* file sizes are enormous, in part because bulky read quality data is mixed in with fasta.
* to extract the fasta sequence, proprietary software called sffinfo must be obtained (or a script written). 
* the reads are fairly short (235 bp or 1-2 exons) and overlapping reads are not assembled.
* no blast server is available to target the reads.

The process works like this (for a biomedical researcher not familiar with linux on a Macintosh):

* download the software sffinfo, place in a folder called SFF
* open a linux command line interface by opening the Terminal program in the Finder's Utilities Folder 
* direct the command line to the folder path using: cd /Users/joeblow/Desktop/SFF
* make ssfinfo recognizable as an executable program: chmod +x sffinfo
* download some data files of interest from NCBI's SRA archive into the SFF foldere
* apply sffinfo to a data file name: ./sffinfo -s EUEMSW405.sff > anyOutputFileName
* split the output file into smaller pieces if needed: split -l 30000 anyOutputFileName

These steps discard the read quality data while retaining the fasta sequences. However that usually results in too large a file for practical purposes such as visualizing entries in ordinary desktop software or submitting to a web tool such as [http://www.proweb.org/proweb/Tools/WU-blast.html uBlast]. Consequently the big fasta file needs to be split up into smaller files, here 30,000 lines each. 

On a typical 454 project, it takes an hour to download one transcript run on a Qwest DSL line. It comes compressed (.gz) but this is automatically decompressed by Stuffit. The resultant .sff file was 282 megabytes. After applying ssfinfo, the file size was reduced to 57 megabytes. After splitting,  35 files result, each 1.6 mbytes and each containing 5800 individual fasta formatted reads. These appear in the SFF folder named as xaa, xab, ... These smaller files can be viewed in any text editor software. However most web tools have a 'browse' option that avoids copy/paste in the Clipboard, simply loading the file invisibly into the web tool for submission.

The SRA is now providing a download button on their [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?cmd=viewer&f=&m=data&s=viewer run browser.] Here the trick is to know the [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?cmd=table&f=run&m=data&s=run run number.]

In the case ofTachyglossus (echidna transcripts), the accession number SRA000241 is a bundle of 12 runs. These must be downloaded individually from the run browser. After unpacking with ssfinfo, they can be concatenated into a single large file since no purpose is served keeping runs separate.

* [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR000647&cmd=viewer&m=data&s=viewer SRR000647] 181241 reads
* [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR000648&cmd=viewer&m=data&s=viewer SRR000648]    378 reads
* [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR000649&cmd=viewer&m=data&s=viewer SRR000649] 204062 reads
* [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR000650&cmd=viewer&m=data&s=viewer SRR000650]   1263 reads
* [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR000651&cmd=viewer&m=data&s=viewer SRR000651]   4498 reads
* [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR000652&cmd=viewer&m=data&s=viewer SRR000652]   1309 reads
* [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR000653&cmd=viewer&m=data&s=viewer SRR000653]   3170 reads
* [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR000654&cmd=viewer&m=data&s=viewer SRR000654]   1663 reads
* [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR000655&cmd=viewer&m=data&s=viewer SRR000655]   2414 reads
* [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR000656&cmd=viewer&m=data&s=viewer SRR000656] 242673 reads
* [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR000657&cmd=viewer&m=data&s=viewer SRR000657]    598 reads
* [http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR000658&cmd=viewer&m=data&s=viewer SRR000658] 242164 reads
* 885433 total reads

Next, it is currently necessary to [http://www.ncbi.nlm.nih.gov/BLAST/download.shtml install a personal Blast server] on your computer so that the database of 454 fasta sequences can be searched with favorite queries without using command lines. That's [http://www.ncbi.nlm.nih.gov/blast/docs/wwwblast.html explained] on the NCBI site. Basically the fasta set has to be reconfigured to a human-unreadable format that the Blast tool reads (formatdb) while still allowing full length sequences involved in Blast matches to be retrieved (fastacmd).

As a practical matter, the vast majority of the 454 reads (98%) do not cover any coding exon. The file can be partitioned into those that do and those that don't. At the level of an individual gene family, procede by concatenating protein sequences for representative members of the family (eg selenoproteins, or opsins). Use that as tBlastn query on the home brew Blast server. Extract all the match lines -- these begin with the word Sbjct in blast output lines. Clean out spaces, gap dashes, and numbers using Jim Kent's  [http://www.soe.ucsc.edu/~kent/protDust/protDust.html Protein Duster] obtaining a concatenate of all translated matches to the gene family. 

Next, to assemble reads for a given protein, use the concatenate as Blastp query against an individual exon of an individual gene of the given gene family at [http://www.proweb.org/proweb/Tools/WU-blast.html uBlast.] This allows ready extraction of exons (though Blast output is by score and so out of sequence order). In comparative genomics, a known already-intronated gene is generally available as paste template. This has the practical effect of allowing visual tiling of reads in the situation where some of the reads have frameshifting indels and no one read covers the exon satisfactorily.




[[Category:Comparative Genomics]]