search

benjjneb / dada2

Accurate sample inference from amplicon data with single nucleotide resolution
http://benjjneb.github.io/dada2/
GNU Lesser General Public License v3.0
464 stars 142 forks source link

Assing taxonomy fails but no error is prompted #126

Closed josemseoane closed 7 years ago

josemseoane commented 7 years ago

Dear Ben,

I have run successfully dada2 and i can see the potential of this tool, so thanks a lot for that. My only issue is that I cannot make assignTaxonomy to work properly on my data.

I am running two samples and all seems to be ok until I get to the taxonomy assignment. I have tried to use only sample and up to 30 and always get the same result, only one archea genus is detected. These are soil communities and I have used the mock community you provided in your example, I don't know if/how that matters or how could I get a more relevant mock community for my kind of samples. I have run the same samples through vamps and it worked ok but I would like to stick to dada2 because integration with other libraries in R is awesome! This is what I get, only reducing the minBoot to crazy low values like 10 will improve the assignment. Could you please give me a suggestion on where to look at? I am totally lost since no error is prompted.

Once more time, thanks a lot in advance for your support!

Jose

seqtab <- makeSequenceTable(mergers[names(mergers) != "Mock"]) The sequences being tabled vary in length.

dim(seqtab) [1] 2 314

Inspect Distribution of sequences lenght

table(nchar(colnames(seqtab)))

230 245 247 248 251 255 271 272 273 298 302 303 304 313 314 316 320 336 339 348 350 360 369 378 379 383 397 1 1 1 1 1 2 13 235 32 1 2 1 1 3 1 1 1 1 3 1 1 2 1 1 2 1 3

Remove bimeras

seqtab.nochim <- removeBimeraDenovo(seqtab, verbose=TRUE) Identified 50 bimeras out of 314 input sequences. dim(seqtab.nochim) [1] 2 264

Calculate contribution of quimeras

sum(seqtab.nochim)/sum(seqtab) [1] 0.9619778

Assign taxonomy:

taxa <- assignTaxonomy (seqtab.nochim, "/home/jmseoane/Documents/Bioinfo/metagenome/Reference/rdp_train_set_14.fa.gz", minBoot=50) colnames(taxa) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus") unname(head(taxa)) [,1] [,2] [,3] [,4] [,5] [,6] [1,] NA NA NA NA NA NA
[2,] NA NA NA NA NA NA
[3,] NA NA NA NA NA NA
[4,] NA NA NA NA NA NA
[5,] "Archaea" NA NA NA NA NA
[6,] NA NA NA NA NA NA
head(seqTab) Error in head(seqTab) : object 'seqTab' not found head(seqtab) GGTATCTAATCCGGTTCGCGACCCCAGCCTTCGTCCCTCACCGTCGGACCCGGTCTGGTCAGACGCCTTCGCTACCGGTGGTCCCCCAAGGATTAATAGATTTCACCCCTACCCCTGGAGTACCTCTGACCTCTCCCGGTCCCAAGCCAAGCAGTATTCCCTGATAGCTTAACCGTTAAGCGGTCAAATTTCCCAAGGAACTTACTCAGCCGGCTACGGACGCTTTAGACCCAATAAGCATGGCTACCACTTGAGCTGCCGGTGTTACCGCGG 152656A-RPCM0138-515fB 11354 P-161691-RPCM0138-162538-8519-A 111600 GGTTTCTAATCCGGTTCGCGACCCCAGCCTTCGTCCCTCACCGTCGGACCCGGTCTGGTCAGACGCCTTCGCTACCGGTGGTCCCCCAAGGATTAATAGATTTCACCCCTACCCCTGGAGTACCTCTGACCTCTCCCGGTCCCAAGCCAAGCAGTATTCCCTGATAGCTTAACCGTTAAGCGGTCAAATTTCCCAAGGAACTTACTCAGCCGGCTACGGACGCTTTAGACCCAATAAGCATGGCTACCACTTGAGCTGCCGGTGTTACCGCGG 152656A-RPCM0138-515fB 4767 P-161691-RPCM0138-162538-8519-A 76423 GGTATCTAATCCCGTTTGCTCCCCTAGCTTTCGCTCCTCAGCGTCAGAGATAGACCAGAAAGCCGCCTTCGCCACTGGTGTTCCTCCCGATATCTACGCATTTCACCGCTACACCGGGAATTCCACTTTCCTCTTCTATCCTCAAGTCTTAAAGTTTCAGATGACCTTATCAGGTTGAGCCTGATAATTTCACACCTGACTTTTAAGACCGCCTACGAGCTCTTTACGCCCAATAATTCCGGATAACGCTTGCTCCTTATGTATTACCGCGG 152656A-RPCM0138-515fB 1409 P-161691-RPCM0138-162538-8519-A 1126 GGTATCTAATCCCGTTTGCTCCCCTAGCTTTCGCTCCTCAGCGTCAGAGATAGACCAGAAAGCCGCCTTCGCCACTGGTGTTCCTCCCGATATCTACGCATTTCACCGCTACACCGGGAATTCCACTTTCCTCTTCTAACCTCAAGTCTTAAAGTTTCAGATGACCTTATCAGGTTAAGCCTGATAATTTCACACCTGACTTTTAAGACCGCCTACGAGCTCTTTACGCCCAATAATTCCGGATAACGCTTGCTCCTTATGTATTACCGCGG 152656A-RPCM0138-515fB 990 P-161691-RPCM0138-162538-8519-A 1233 GGTTTCTAATCCCGTTTGCTCCCCTAGCTTTCGCTCCTCAGCGTCAGAGATAGACCAGAAAGCCGCCTTCGCCACTGGTGTTCCTCCCGATATCTACGCATTTCACCGCTACACCGGGAATTCCACTTTCCTCTTCTATCCTCAAGTCTTAAAGTTTCAGATGACCTTATCAGGTTGAGCCTGATAATTTCACACCTGACTTTTAAGACCGCCTACGAGCTCTTTACGCCCAATA

benjjneb commented 7 years ago

DADA2 does not currently check the orientation of your sequences, and the ones you posted here seem to be in the reverse order (3' to 5'). If you reverse complement your sequences, they should classify correctly.

colnames(seqtab) <- dada2:::rc(colnames(seqtab))
josemseoane commented 7 years ago

Hi Ben,

You were right again, this fixed my problem. Once more time thank you very much for your time and for sharing your work!