OTU进化树关系，从qiime2的结果文件入手。更多知识分享请到 https://zouhua.top/。

生成树文件和进化关系文件

  • rep_seq_k5_cln.tree
  • rep_seqs_k5.tax
  1. # install software 
  2. conda install -n qiime -c bioconda qiime -y 
  3. source activate py27
  4. conda install clustalo -y
  6. # step1: count -> relative abundance
  7. biom normalize-table -i feature-table.biom -r -o feature-table_norm.biom
  9. # step2: filter abundance 0.5%
  10. filter_otus_from_otu_table.py --min_count_fraction 0.005 -i feature-table_norm.biom -o feature-table_norm_k5.biom
  12. # step3: get fasta
  13. #qiime tools export --input-path rep-seqs.qza --output-path rep-seqs
  14. filter_fasta.py -f rep_seq.fa -o rep_seq_k5.fa -b feature-table_norm_k5.biom
  16. # step4: multip alignment
  17. clustalo -i rep_seq_k5.fa -o rep_seq_k5_clus.fa --seqtype=DNA --full --force --threads=30
  19. # step5: fastree build tree file
  20. make_phylogeny.py -i rep_seq_k5_clus.fa -o rep_seq_k5.tree
  22. # step6: delete ' signal for R ggtree
  23. sed "s/\'//g" rep_seq_k5.tree > rep_seq_k5_cln.tree
  25. # step7: get annotation information
  26. grep ">" rep_seq_k5_clus.fa | sed 's/>//g' > rep_seq_k5_clus.id
  27. awk 'BEGIN{OFS="\t";FS="\t"} NR==FNR {a[$1]=$0} NR>FNR {print a[$1]}' taxonomy.tsv rep_seq_k5_clus.id | sed 's/; /\t/g' | cut -f 1-5 | sed 's/;/\t/g' | cut -f 1-5 > rep_seqs_k5.tax

可视化

因为qiime2的ASV名称是复杂编码，所以需要对其进行修改，但最好在得到otu table时候就同一修改

  1. library(dplyr)
  2. library(tibble)
  3. library(ggtree)
  4. library(ggplot2)
  6. # load data 
  7. tree <- read.tree("rep_seq_k5_cln.tree")
  8. tax <- read.table("rep_seqs_k5.tax", row.names=1)
  10. # curation
  11. colnames(tax) <- c("kingdom","phylum","class","order")
  12. otuid <- data.frame(OTUID=tree$tip.label, 
  13.                     OTU_ID=paste0("OTU_", seq(1:length(tree$tip.label))))
  14. tax <- inner_join(tax %>% rownames_to_column("OTUID"), otuid, by = "OTUID") %>%
  15.   column_to_rownames("OTU_ID") %>%
  16.   dplyr::select(-OTUID)
  17. tree$tip.label <- paste0("OTU_", seq(1:length(tree$tip.label)))

plot 1

groupInfo <- split(row.names(tax), tax$phylum) 
tree <- groupOTU(tree, groupInfo)
ggtree(tree, aes(color=group))+  
  theme(legend.position = "right")+
  geom_tiplab(size=3)

数据分析：扩增子进化树 - 图1

plot2

ggtree(tree, layout="fan", ladderize = FALSE, size=1.2, 
       branch.length = "none", aes(color=group))+
  geom_tiplab2(size=4)+ 
  theme(legend.position = "right")

数据分析：扩增子进化树 - 图2

参考

  1. 扩增子进化树分析

参考文章如引起任何侵权问题，可以与我联系，谢谢。