写在前面Qiime2是大佬们集中已有的资源和对扩增子领域的思考而成的一个集成的,一个大杂烩的软件。但是,事实上Qiime也是一个大杂烩,不过是命令的大杂烩,将全部不同功能的命令组合到一起,超过200个命令构成了Qiime的全部功能。Qiime2更像是一个在芸芸众生中都进行了一个尝试后,新诞生出来的一个成熟的工具,将不同功能的命令封装起来,从底层解决文件大小问题,解决数据共享问题,解决可重复计算和分析问题等。 从Qiime到Qiime2,遵循事物发展规律,理解从1到2的转变,就可以了解微生物组从开始到成熟的技术演变。从12年到21年,这些年见,最早的一批人已经从摸索阶段到了领域标杆,那么下一个十年,机会在哪里的?不管怎样,先加油!下一个十年属于我们。 环境配置安装conda+Qiime2 # 下载最新版miniconda2 wget -c https://repo./miniconda/Miniconda3-latest-Linux-x86_64.sh # 运行安装程序 bash Miniconda3-latest-Linux-x86_64.sh 数据导入library(tidyverse) #这里我们提取复合要求的指定文件名称
b =dir("./seq/", pattern = c(""), full.names = F, ignore.case = TRUE) b
f = b[b %>% str_detect("_1")] f <- paste("$PWD/seq/",f,sep = "")
r = b[b %>% str_detect("_2")] r <- paste("$PWD/seq/",r,sep = "")
id = b[b %>% str_detect("_2")] %>% strsplit( "_") %>% sapply( `[`, 1)
dat = data_frame( 'sample-id' = id, 'forward-absolute-filepath' = f, 'reverse-absolute-filepath' = r
)
write_tsv(dat,"./manifest.tsv") time qiime tools import \ --type 'SampleData[PairedEndSequencesWithQuality]' \ --input-path smanifest.tsv \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33V2 vsearch+debur# 双端序列合并 time qiime vsearch join-pairs \ --i-demultiplexed-seqs demux.qza \ --o-joined-sequences demux-joined.qza
# 双端去接头 time qiime cutadapt trim-paired \ --i-demultiplexed-sequences demux.qza \ --p-cores 9 \ --p-front-f AACMGGATTAGATACCCKG \ --p-front-r GGAAGGTGGGGATGACGT \ --p-error-rate 0.1 \ --o-trimmed-sequences trimmed-demux.qza \ --verbose
# 序列质控 time qiime quality-filter q-score-joined \ --i-demux demux-joined.qza \ --o-filtered-sequences demux-joined-filtered.qza \ --o-filter-stats demux-joined-filter-stats.qza
# deblur去噪生成特征表 time qiime deblur denoise-16S \ --i-demultiplexed-seqs demux-joined-filtered.qza \ --p-trim-length 250 \ --p-sample-stats \ --o-representative-sequences rep-seqs.qza \ --o-table table.qza \ --o-stats deblur-stats.qza DADA2time qiime dada2 denoise-paired --i-demultiplexed-seqs demux.qza --p-n-threads 0 \ --p-trim-left-f 29 --p-trim-left-r 18 --p-trunc-len-f 0 --p-trunc-len-r 0 \ --o-table table.qza \ --o-representative-sequences rep-seqs.qza \ --o-denoising-stats denoising-stats.qza 外部导入特征表和代表序列导入文件相对麻烦,qiime tools import命令中的参数 —type 多达几十种,所以这里列出来了常用的几种。 # 上传OTU表otutab.txt和代表序列otus.fa # 转换文本为Biom1.0 biom convert -i ./input_otutab/otutab.txt -o ./input_otutab/otutab.biom \ --table-type="OTU table" --to-json # 导入特征表 qiime tools import --input-path ./input_otutab/otutab.biom \ --type 'FeatureTable[Frequency]' --input-format BIOMV100Format \ --output-path ./input_otutab/table.qza # 导入代表序列 qiime tools import --input-path ./input_otutab/otus.fa \ --type 'FeatureData[Sequence]' \ --output-path ./input_otutab/rep-seqs.qza # 注意如果不知道tax格式,那就看下qiime2导出来的是什么格式,同样也可以被导入 qiime tools import --input-path ./input_otutab/taxonomy.txt \ --type 'FeatureData[Taxonomy]'\ --output-path ./input_otutab/tax.qza 特征表格基本统计# 特征表和代表序列统计 qiime feature-table summarize \ --i-table ./input_otutab/table.qza \ --o-visualization ./input_otutab/table.qzv \ --m-sample-metadata-file metadata.tsv
qiime feature-table tabulate-seqs \ --i-data ./input_otutab/rep-seqs.qza \ --o-visualization ./input_otutab/rep-seqs.qzv 导出数据这里十分简单,只需要知道qiime tools export命令接即可,并且将其看作一个解压缩的软件。 #--导出OTU表格 qiime tools export --input-path table.qza --output-path exported-feature-table
biom convert -i exported-feature-table/feature-table.biom \ -o exported-feature-table/table.txt \ --to-tsv qiime tools export --input-path ./rep-seqs.qza --output-path rep set
qiime tools export --input-path ./rooted-tree.qza --output-path tree Alpha和beta多样性分析# 构建进化树用于多样性分析 53s time qiime phylogeny align-to-tree-mafft-fasttree \ --i-sequences rep-seqs.qza \ --o-alignment aligned-rep-seqs.qza \ --o-masked-alignment masked-aligned-rep-seqs.qza \ --o-tree unrooted-tree.qza \ --o-rooted-tree rooted-tree.qza
# 计算核心多样性 8s time qiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table.qza \ --p-sampling-depth 10000 \ --m-metadata-file metadata.tsv \ --output-dir core-metrics-results
# Beta多样性组间显著性分析和可视化 # 可选的beta指数有 unweighted_unifrac、bray_curtis、weighted_unifrac、jaccard # 指定分组是减少计算量,转换检验非常耗时 distance=weighted_unifrac column=Group qiime diversity beta-group-significance \ --i-distance-matrix core-metrics-results/${distance}_distance_matrix.qza \ --m-metadata-file metadata.tsv \ --m-metadata-column ${column} \ --o-visualization core-metrics-results/${distance}-${column}-significance.qzv \ --p-pairwise
# Alpha多样性组间显著性分析和可视化 # 可选的alpha指数有 faith_pd、shannon、observed_otus、evenness index=observed_otus qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/${index}_vector.qza \ --m-metadata-file metadata.tsv \ --o-visualization core-metrics-results/${index}-group-significance.qzv
# Alpha稀疏取线 qiime diversity alpha-rarefaction \ --i-table table.qza \ --i-phylogeny rooted-tree.qza \ --p-max-depth 20000 \ --m-metadata-file metadata.tsv \ --o-visualization alpha-rarefaction.qzv # 结果有observed_otus, shannon, 和faith_pd三种指数可选 训练注释库# 导入参考序列 qiime tools import \ --type 'FeatureData[Sequence]' \ --input-path 99_otus.fasta \ --output-path 99_otus.qza
# 导入物种分类信息 qiime tools import \ --type 'FeatureData[Taxonomy]' \ --input-format HeaderlessTSVTaxonomyFormat \ --input-path 99_otu_taxonomy.txt \ --output-path 99_ref-taxonomy.qza
time qiime feature-classifier fit-classifier-naive-bayes \ --i-reference-reads 99_otus.qza \ --i-reference-taxonomy 99_ref-taxonomy.qza \ --o-classifier gg_99_all_classifier.qza 物种组成分析# 物种注释 1m42s time qiime feature-classifier classify-sklearn \ --i-classifier ./training-feature-classifiers/gg_99_all_classifier.qza \ --i-reads rep-seqs.qza \ --o-classification taxonomy.qza
qiime metadata tabulate \ --m-input-file taxonomy.qza \ --o-visualization taxonomy.qzv
# 堆叠柱状图展示 qiime taxa barplot \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file sample-metadata.tsv \ --o-visualization taxa-bar-plots.qzv 根际互作生物学研究室 简介根际互作生物学研究室是沈其荣教授土壤微生物与有机肥团队下的一个关注于根际互作的研究小组。本小组由袁军副教授带领,主要关注:1.植物和微生物互作在抗病过程中的作用;2 环境微生物大数据整合研究;3 环境代谢组及其与微生物过程研究体系开发和应用。团队在过去三年中在 isme J, Microbiome, PCE,SBB,Horticulture Research等期刊上发表了多篇文章。欢迎关注 微生信生物 公众号对本研究小组进行了解。
|