dada2物种分类鉴定

生信喵

物种分类鉴定

1. #方法一:dada2自带注释函数,官网下载
   
2. #https://benjjneb.github.io/dada2/training.html
   
3. # taxa <- assignTaxonomy(seqtab.nochim, "MiSeq_SOP/tax/silva_nr99_v138.1_train_set.fa.gz", multithread=12)
   
4. # taxa <- addSpecies(taxa, "MiSeq_SOP/tax/silva_species_assignment_v138.1.fa.gz")
   
5. #
   
6. # taxa.print <- taxa
   
7. # rownames(taxa.print) <- NULL
   
8. # head(taxa.print)
   
9. 
   
10. #方法二： DECIPHER
    
11. download.file(url="http://www2.decipher.codes/Classification/TrainingSets/SILVA_SSU_r138_2019.RData", destfile="SILVA_SSU_r138_2019.RData")
    
12. 
    
13. #加载CECIPHER SILVA库
    
14. load("SILVA_SSU_r138_2019.RData")
    
15. 
    
16. #加载DECIPHER
    
17. library(DECIPHER)
    
18. packageVersion("DECIPHER")
    
19. 
    
20. #创建dna序列
    
21. dna <- DNAStringSet(getSequences(seqtab.nochim))
    
22. 
    
23. ids <- IdTaxa(dna, trainingSet, strand="top", processors=12, verbose=FALSE)
    
24. ranks <- c("domain", "phylum", "class", "order", "family", "genus", "species") # ranks of interest
    
25. # Convert the output object of class "Taxa" to a matrix analogous to the output from assignTaxonomy
    
26. taxid <- t(sapply(ids, function(x) {
    
27.   m <- match(ranks, x$rank)
    
28.   taxa <- x$taxon[m]
    
29.   taxa[startsWith(taxa, "unclassified_")] <- NA
    
30.   taxa
    
31. }))
    
32. colnames(taxid) <- ranks
    
33. rownames(taxid) <- getSequences(seqtab.nochim)

复制代码

评估结果

1. unqs.mock <- seqtab.nochim["Mock",]
   
2. unqs.mock <- sort(unqs.mock[unqs.mock>0], decreasing=TRUE)
   
3. cat("DADA2 inferred", length(unqs.mock), "sample sequences present in the Mock community.\n")
   
4. 
   
5. mock.ref <- getSequences(file.path(path, "HMP_MOCK.v35.fasta"))
   
6. match.ref <- sum(sapply(names(unqs.mock), function(x) any(grepl(x, mock.ref))))
   
7. cat("Of those,", sum(match.ref), "were exact matches to the expected reference sequences.\n")

复制代码

画树

1. library(msa)
   
2. 
   
3. seqs <- getSequences(seqtab.nochim)
   
4. names(seqs) <- seqs
   
5. mult <- msa(seqs, method="ClustalW", type="dna", order="input")
   
6. 
   
7. library("phangorn")
   
8. phang.align <- as.phyDat(mult, type="DNA", names=getSequence(seqtab))
   
9. dm <- dist.ml(phang.align)
   
10. treeNJ <- NJ(dm) # Note, tip order != sequence order
    
11. fit = pml(treeNJ, data=phang.align)
    
12. ## negative edges length changed to 0!
    
13. fitGTR <- update(fit, k=4, inv=0.2)
    
14. fitGTR <- optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE,
    
15.                     rearrangement = "stochastic", control = pml.control(trace = 0))
    
16. fitGTR$tree

复制代码