search

gjospin / PhyloSift

Phylogenetic and taxonomic analysis for genomes and metagenomes
82 stars 17 forks source link

"Walking" Probability distributions not reaching 1 in sequence_taxa_summary.txt #211

Closed hollybik closed 12 years ago

hollybik commented 12 years ago

Noticing some discrepancies in the probability distributions that are reported for placed reads in sequence_taxa_summary.txt

After "walking" up the tree and summing probabilities, the top probability listed for a unique sequence read doesn't always reach 1. Some examples of top probability reported:

HWI-ST640_0078:4:21:13990:121822#AACACCTATCTC/1 32523 no rank TETRAPODA 0.833333333333333

HWI-ST640_0078:4:21:15708:198355#AACACCTATCTC/1 1 no rank ROOT 0.962962962962962

HWI-ST640_0078:4:42:2304:134799#AACACCTATCTC/1 1 no rank ROOT 0.916666666666667 mtDNA_ND1_aa_aligned

HWI-ST640_0078:4:66:10938:60499#AACACCTATCTC/1 32523 no rank TETRAPODA 0.857142857142857

Alternatively, some placed reads have a LOT of taxa reported with a probability of 1 - to me, this looks like PhyloSift is continuing to walk up the tree even after the probability sum has readched its maximum. An example:

HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 32523 no rank TETRAPODA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 32525 no rank THERIA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 33316 no rank COELOMATA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 117571 no rank EUTELEOSTOMI 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 1 no rank ROOT 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 33208 kingdom METAZOA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 117570 no rank TELEOSTOMI 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 2759 superkingdom EUKARYOTA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 7742 no rank VERTEBRATA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 33213 no rank BILATERIA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 314145 superorder LAURASIATHERIA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 131567 no rank CELLULAR ORGANISMS 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 32524 no rank AMNIOTA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 33154 no rank OPISTHOKONTA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 89593 subphylum CRANIATA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 33511 no rank DEUTEROSTOMIA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 7711 phylum CHORDATA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 6072 no rank EUMETAZOA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 7776 superclass GNATHOSTOMATA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 8287 no rank SARCOPTERYGII 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 40674 class MAMMALIA 1
HWI-ST640_0078:4:61:1515:58431#AACACCTATCTC/1 9347 no rank EUTHERIA 1

Examples are all taken from outputs on Edhar at: /home/hollybik/phylosift_devel_20120511/PS_temp/pool5_AACACCT_test.fastq

koadman commented 12 years ago

Upon closer inspection there are no bugs here.

The first issue, that of probabilities not reaching 1 in the summary, results from a portion of the probability mass being unclassifiable. In the companion output file sequence_taxa.txt, it can be seen that the remaining portion of the probability mass is listed as unclassifiable. Thus there is no taxonomic level at which it is possible to achieve a summed probability mass of 1.

The second issue is simply a product of the NCBI taxonomy structure. We could consider an alternate output format where levels are no longer reported once probability mass reaches 1, but creating summaries from such a format may be more difficult. For example in such a format the phylum CHORDATA probability of 1 would not be reported for HWI-ST640_0078 so any script/person/spreadsheet wishing to quickly generate a phylum level summary would have to also parse the NCBI taxonomy whereas that is not required in the current format -- one need only grep for phylum in the appropriate column.