好久没更了，这半年心情很糟糕，也在不断地调整自己。希望借助更新来转移一下注意力，毕竟啥都没有学到的知识重要~~~
刚接触微生物组这个领域的时候，就看了Rob Knight的一系列视频和相关文章，很多分析都离不开他们实验室开发的一些软件和数据库，尤其是用于16S rRNA分析的QIIME和升级后焕然一新的QIIME2。因为之前学过python，而QIIME2本身就是基于python开发的，所以我在刚入门微生物组分析的时候就学习了QIIME2的使用。以下是相关的代码，仅供参考！
- 文件输入与处理
注意，这里的文件是已经质控、去引物以及拼接过的cleandata，双端序列只需稍改以下参数设置即可，具体详见QIIME2的user documentation。
source activate qiime2
work_dir=/Share/user/00.digestive
temp=/Share/user/00.digestive/temp
phyloseq=/Share/user/00.digestive/phyloseq
cd ${work_dir}
# import data
time qiime tools import \
--type 'SampleData[SequencesWithQuality]' \
--input-path manifest.txt \
--output-path ${temp}/demux.qza \
--input-format SingleEndFastqManifestPhred33V2
# summarize data
time qiime demux summarize \
--i-data ${temp}/demux.qza \
--o-visualization ${temp}/demux.qzv
# DADA2
time qiime dada2 denoise-single \
--i-demultiplexed-seqs ${temp}/demux.qza \
--p-trunc-len 72 \
--p-n-threads 32 \
--o-representative-sequences ${temp}/rep-seqs.qza \
--o-table ${temp}/table.qza \
--o-denoising-stats ${temp}/denoising-stats.qza
# visualization feature-table
time qiime feature-table summarize \
--i-table ${temp}/table.qza \
--o-visualization ${temp}/table.qzv \
--m-sample-metadata-file metadata.txt
# visualization rep-seqs
time qiime feature-table tabulate-seqs \
--i-data ${temp}/rep-seqs.qza \
--o-visualization ${temp}/rep-seqs.qzv
# filtering feature-table based on frequency and number of OTU
time qiime feature-table filter-features \
--i-table ${temp}/table.qza \
--p-min-frequency 10 \
--p-min-samples 2 \
--o-filtered-table ${temp}/filtered-table.qza
# filtering rep-seqs based on feature-table
time qiime feature-table filter-seqs \
--i-data ${temp}/rep-seqs.qza \
--i-table ${temp}/filtered-table.qza \
--o-filtered-data ${temp}/filtered-seqs.qza
# export feature-table
time qiime tools export \
--input-path ${temp}/filtered-table.qza \
--output-path ${phyloseq}
# export rep-seqs
time qiime tools export \
--input-path ${temp}/filtered-seqs.qza \
--output-path ${phyloseq}
cd ${phyloseq}
biom convert -i feature-table.biom -o filtered-table.tsv --to-tsv
sed -i '1d' filtered-table.tsv
sed -i 's/#OTU ID/ASV ID/' filtered-table.tsv
less dna-sequences.fasta |paste - -|sed '1i ASVID,seq' > rep.fa
conda deactivate
- 文件格式转换
这里需要切换到R base底下进行相关操作。
# 设置工作路径
setwd("/Share/user/00.digestive/phyloseq")
# 查看路径下的所有文件
list. files()
# 载入包
library(tidyverse)
library(magrittr)
library(stringr)
# 转化文件格式
otu <- read.delim("filtered-table.tsv",check.names = FALSE,header = T,sep="\t")
rown <- paste0("ASV",seq_len(nrow(otu)))
otu[,1] <- rown
colnames(otu)[1] <- paste0("ASV",colnames(data)[1])
write.table (otu,file ="ASV_table.tsv",sep ="\t",row.names = F,quote = F) #更改ASV表名称
rep <- read.delim("rep.fa",check.names = FALSE, row.names = 1)%>%
set_rownames(paste0(">ASV", seq_len(nrow(.))))
write.table (rep,file ="rep.xls", sep ="\t", row.names = T,quote = F) #更改特征序列名称
- 物种注释
因为greengenes数据库好久没更新了（当然，前几天Rob Knight课题组更新了），所以我更常用Silva数据库，而且QIIME2还非常贴心地帮你训练好了相似阈值为99%的V4区域的数据库集，更加方便。
source activate qiime2
work_dir=/Share/user/
temp=/Share/user/00.digestive/temp
phyloseq=/Share/user/00.digestive/phyloseq
cd ${phyloseq}
# changing format
sed -i '1d' rep.xls
sed -i 's/\r//g' rep.xls
less rep.xls|sed 's/\t/\n/g' > filtered-seqs.fasta
mv filtered-seqs.fasta ${temp}/
time qiime tools import \
--type 'FeatureData[Sequence]' \
--input-path ${temp}/filtered-seqs.fasta \
--output-path ${temp}/filtered-seqs.qza
cd ${work_dir}
# classifier
time qiime feature-classifier classify-sklearn \
--i-classifier silva-138-99-nb-classifier.qza \
--i-reads ${temp}/filtered-seqs.qza \
--o-classification ${temp}/taxonomy.qza
# export taxonomy
time qiime tools export \
--input-path ${temp}/taxonomy.qza \
--output-path ${phyloseq}
# remove Eukaryota,Mitochondria,Chloroplast,Unassigned
conda deactivate
- 系统发育树生成以及导出下游分析所需的文件
source activate qiime2
work_dir=/Share/user/00.digestive
temp=/Share/user/00.digestive/temp
phyloseq=/Share/user/00.digestive/phyloseq
cd ${phyloseq}
less rep.xls|sed 's/\t/\n/g' > filtered-seqs.fasta
less taxonomy.tsv|sed '1i ASV,domain,phylum,class,order,family,genus,species'|sed 's/,/\t/g'|sed 's/;/\t/g' |sed 's/[a-z]__//g'> taxa.xls
biom convert -i ASV_table.tsv -o feature-table.biom --to-hdf5 --table-type='OTU table'
mv filtered-seqs.fasta ${temp}/
mv feature-table.biom ${temp}/
cd ${work_dir}
# convert seq.fasta to qza
time qiime tools import \
--type 'FeatureData[Sequence]' \
--input-path temp/filtered-seqs.fasta \
--output-path temp/filtered-seqs.qza
# convert taxonomy.tsv to qza
time qiime tools import \
--type FeatureData[Taxonomy] \
--input-format HeaderlessTSVTaxonomyFormat \
--input-path phyloseq/taxonomy.tsv \
--output-path temp/taxonomy.qza
# convert ASV.txt to qza
time qiime tools import \
--input-path temp/feature-table.biom \
--type 'FeatureTable[Frequency]' \
--input-format BIOMV210Format \
--output-path temp/filtered-table.qza
# construct phylogeny tree
time qiime phylogeny align-to-tree-mafft-fasttree \
--i-sequences temp/filtered-seqs.qza \
--o-alignment temp/aligned-rep-seqs.qza \
--o-masked-alignment temp/masked-aligned-rep-seqs.qza \
--o-tree temp/unrooted-tree.qza \
--o-rooted-tree temp/rooted-tree.qza
# export rooted-tree
time qiime tools export \
--input-path temp/rooted-tree.qza \
--output-path phyloseq/rooted_tree
# visualization ASV
time qiime feature-table summarize \
--i-table temp/filtered-table.qza \
--o-visualization temp/filtered-table.qzv \
--m-sample-metadata-file metadata.txt
# visualizaion taxonomy
time qiime metadata tabulate \
--m-input-file temp/taxonomy.qza \
--o-visualization temp/taxonomy.qzv
# visualization taxa
time qiime taxa barplot \
--i-table temp/filtered-table.qza \
--i-taxonomy temp/taxonomy.qza \
--m-metadata-file metadata.txt \
--o-visualization temp/taxa-bar-plots.qzv
# visualizaion seq
time qiime feature-table tabulate-seqs \
--i-data temp/filtered-seqs.qza \
--o-visualization temp/filtered-seqs.qzv
conda deactivate
QIIME2也可以直接生成alpha和beta多样性指数结果，并进行可视化，但是考虑到个性化分析的需要，以及QIIME2 view的界面实在是不好看，因此我都是利用R和生成的四个文件进行后续的下游分析。
