1. 安装虚拟机及QIIME1
2. 虚拟机增强功能及共享文件夹
简单，百度都有，目的在于实现win系统和linux系统可以自由复制粘贴、全屏和文件共享
3. 准备数据文件
图1. 下载文件
params.txt
map.txt
diversityParams.txt
4. OTU分析
4.1. map文件检查
validate_mapping_file.py -o vmf-map/ -m map.txt
检查结果
结果：No errors or warnings were found in mapping file.
同时生成一个vmf-map文件夹
图2. 生成vmf-map文件夹
4.2. 序列拼接
overlap大于6bp且错配率小于8%才有可能完成匹配
join_paired_ends.py -f SRR1370913_1.fastq.gz -r SRR1370913_2.fastq.gz -o SRR1370913
join_paired_ends.py -f SRR1370914_1.fastq.gz -r SRR1370914_2.fastq.gz -o SRR1370914
join_paired_ends.py -f SRR1370915_1.fastq.gz -r SRR1370915_2.fastq.gz -o SRR1370915
join_paired_ends.py -f SRR1370920_1.fastq.gz -r SRR1370920_2.fastq.gz -o SRR1370920
拼接命令必选参数
生成文件保存在SRR1230913、SRR1230914、SRR1230915、SRR1230920文件夹内
生成拼接后文件夹
jion文件内容查看
4.3. 过滤
(1)去除过短序列；（2）N charater序列；（3）转化为fasta文件
split_libraries_fastq.py -i SRR1370913/fastqjoin.join.fastq,SRR1370914/fastqjoin.join.fastq,SRR1370915/fastqjoin.join.fastq,SRR1370920/fastqjoin.join.fastq --sample_id SRR1370913,SRR1370914,SRR1370915,SRR1370920 -o slout/ -m map.txt -q 19 --barcode_type 'not-barcoded' --phred_offset=33
过滤命令必选项
过滤命令可选项
所有的文件信息保存在seqs.fna文件内
运行结果
4.3.1. log文件
log文件内容
log文件记录了输入的文件信息及有效序列数目等信息（可以整理为一个质控表格）
4.4. 去除嵌合体
嵌合体的去除与否目前仍存在争议，如果不需要去除，可直接跳过该步骤。
4.5. OTU分析和注释
pick_open_reference_otus.py -o otus/ -i slout/seqs.fna -p params.txt
包括以下四步：
OTU分析
4.5.1. 结果解读
outs文件
index.html
4.5.2. OTU数目
# 生成简要内容
biom summarize-table -i otus/otu_table_mc2_w_tax_no_pynast_failures.biom
# 生成简要内容并保存
biom summarize-table \
-i otus/otu_table_mc2_w_tax_no_pynast_failures.biom \
-o otu_count.txt
OTU数目
4.5.3. 转化.biom文件为.txt
biom convert -i otus/otu_table_mc2_w_tax_no_pynast_failures.biom -o otus/table.from_biom_w_taxonomy.txt --to-tsv --header-key taxonomy
生成txt文件
4.6 Venn图及热图
4.6.1 绘制Venn图
#install.packages("VennDiagram")
sampleCol=c("SRR1370913","SRR1370914","SRR1370915","SRR1370920")
library(VennDiagram)
rt=read.table("otus/table.from_biom_w_taxonomy.txt",sep="\t",header=T,skip=1,row.names=1,comment.char = "")
list1=list()
for(i in sampleCol){
rt1=rt[rt[,i]!=0,]
print(head(rownames(rt1)))
list1[[i]]=rownames(rt1)
}
venn.diagram(list1,filename="venn.tiff",fill=rainbow(length(sampleCol)))
4.6.2 绘制热图
#install.packages("gplots")
lineNum=100
library('gplots')
rt=read.table("otus/table.from_biom_w_taxonomy.txt",sep="\t",header=T,skip=1,row.names=1,comment.char = "")
sampleNum=ncol(rt)-1
rt=rt[,1:sampleNum]
rt1=rt[order(rowSums(rt),decreasing=TRUE),]
rt2=rt1[1:lineNum,]
y=as.matrix(log10(rt2+1))
pdf(file="heatmap.pdf",height=12)
par(oma=c(3,3,3,5))
heatmap.2(y,col='greenred',trace="none",cexCol=1)
dev.off()
5. α多样性和β多样性
core_diversity_analyses.py -o cdout/ -i otus/otu_table_mc2_w_tax_no_pynast_failures.biom -m map.txt -t otus/rep_set.tre -e 193760 -p diversityParams.txt
包括四个小程序：
1.single rarefaction.py------Rarify the OTU table
2.beta diversity.py-------beta多样性
3.principal coordinates.py------PCoA图形化
4.alpha diversity.py-------Alpha多样性
5.make rarefaction plots.py-------稀释曲线
6.summarize taxa.py------Summarize Taxonomy
7.plot taxa summary.py------Taxa summary plots
其中，193685通过以下命令获得：
biom summarize-table -i otus/otu_table_mc2_w_tax_no_pynast_failures.biom
193760的由来
5.1 结果解读
α及β多样性
生物学分类
α多样性指标
PCoA
6. 系统树构建
系统发生树（英文：phylogenetic tree或evolutionary tree）是表明被认为具有共同祖先的各物种相互间演化关系的树，又被称为系统发育树、系统演化树、系统进化树、种系发生树、演化树、进化树、系统树。它用来表示系统发生研究的结果，用它描述物种之间的进化关系。
6.1. 主要步骤
1.filter otus from otu table.py------挑选丰度>1‰的OTU，获得otu_table_mc2_0.001_fraction.biom文件
filter_otus_from_otu_table.py -i otus/otu_table_mc2_w_tax_no_pynast_failures.biom -o otus/otu_table_mc2_0.001_fraction.biom --min_count_fraction 0.001
2.filter_fasta.py------获得丰度>1‰OTU的fasta文件，获得rep_set_0.001_fraction.fasta文件
filter_fasta.py -f otus/pynast_aligned_seqs/rep_set_aligned_pfiltered.fasta -b otus/otu_table_mc2_0.001_fraction.biom -o otus/rep_set_0.001_fraction.fasta
3.make_phylogeny.py-------构建系统树，获得rep_set_0.001_fraction.tre文件
make_phylogeny.py -i otus/rep_set_0.001_fraction.fasta -o otus/rep_set_0.001_fraction.tre
4.ggtree.R------系统树可视化
如果R版本较低，可以将otus\table.from_0.001_fraction.txt和otus\rep_set_0.001_fraction.tre复制到windows下进行分析。
#if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
#BiocManager::install()
#BiocManager::install('ggtree')
setwd("C:/Users/DELL/Desktop/tree")
library("ggplot2")
library("ggtree")
library("colorspace")
cls=list()
rt=read.table("table.from_0.001_fraction.txt",sep="\t",header=T,skip=1,row.names=1,comment.char = "")
for(i in 1:nrow(rt)){
otu=rownames(rt[i,])
phylum=strsplit(as.character(rt$taxonomy[i]),"\\; |p\\_\\_")[[1]][3]
cls[[phylum]]=c(cls[[phylum]], otu)
}
phylumNames=names(cls)
phylumNum=length(phylumNames)
tree <- read.tree("rep_set_0.001_fraction.tre")
tree <- groupOTU(tree, cls)
# 这部分代码有问题，后续再更改！
pdf(file="circosTree.pdf")
ggtree(tree, layout="fan", ladderize = FALSE, branch.length="none", aes(color=group)) +
scale_color_manual(values=c(rainbow_hcl(phylumNum+1)), breaks=1:phylumNum, labels=phylumNames ) + theme(legend.position="right") +
geom_text(aes(label=paste(" ",gsub("\\d+\\.\\d+|New\\.|Reference|CleanUp\\.","",label),sep=""), angle=angle+90), size=2.2)
dev.off()
pdf(file="heatmapTree.pdf",width=15,height=18)
p <- ggtree(tree, ladderize = FALSE, branch.length="none", aes(color=group)) +
scale_color_manual(values=c(rainbow_hcl(phylumNum+1)), breaks=1:phylumNum, labels=phylumNames ) + theme(legend.position="left") +
geom_text(aes(label=paste(" ",gsub("\\d+\\.\\d+|New\\.|Reference|CleanUp\\.","",label),sep="")), size=2.2)
gplot(p, log(rt[,1:ncol(rt)-1]+1), font.size=4)
dev.off()
5 OTU差异分析
differential_abundance.py -i otus/otu_table_mc2_w_tax_no_pynast_failures.biom -o diff_otus.txt -m map.txt -c SampleType -x normal -y tumor -a DESeq2_nbinom
diff_otus
网友评论