0044-【宏基因组】-16S分析qiime1极简教程
1. 数据介绍
FASTQ数据下载
- EBI
- SRA
数据地址:
https://www.ebi.ac.uk/ena/data/view/PRJNA246490
- Study: PRJNA246490
- 16S rRNA sequencing of soil microbial community metagenome 土壤微生物
- SRP041809 副登录号
- Description:The study goal was to investigate relationships between soil microbial community taxonomy, soil geochemistry, and enzyme activities for two biogeographically different soil types, a cabot very stoney silt loam located in VT, USA and a skerry fine sandy loam located in NH, USA.
研究对象为:位于美国佛蒙特州VT的cabot非常粉砂壤土和位于美国新罕布什尔州的斯凯里细砂壤土。
研究目的:看两个不同的微生物群落的分类特征
总样品:110个
Sample:SRS631664-SRS631753
Experiment:SRX585017-SRX585126
Run:SRR1370913-SRR1371022
将TEXT
文本,提取网址部分,在linux端,使用wget下载
wget -i download.txt
2.样品信息文件 (map.txt)
#SampleID BarcodeSequence LinkerPrimerSequence SampleType Description
SRR1370913 AGCTGACTAGTC GTGCCAGCMGCCGCGGTAA tumor SRR1370913
SRR1370914 AGCTGACTAGTG GTGCCAGCMGCCGCGGTAA tumor SRR1370914
SRR1370915 AGCTGACTAGTA GTGCCAGCMGCCGCGGTAA normal SRR1370915
SRR1370920 AGCTGACTAGTT GTGCCAGCMGCCGCGGTAA normal SRR1370920
- 第一列必须是#SampleID,为可区分的数字,字母或 。每一行的这个值
应该是唯一的。 - 第二列必须为BarcodeSequence。每一行的这个值应该是唯一的。
- 第三列必须为LinkerPrimerSequence,为扩增样品的引物。
- 后面的列可以根据样品的特点加以描述,但是每一列必须包含至少两个值,
就是指定样品信息,就可以在这个位置设置了。 - 最后一列必须是Description。每一个样品不一样的地方,必须不一样额。
样品已经经过barocde拆分为单独的fastq,或者后续按barcode、引物的长度进行切除人工序列,这两列可以填8个A、或18个A作为占位。
上面的分组,只是练习需要,非真实。此列可以不要。
3.激活环境
source activate qiime1
4.打印系统版本
print_qiime_config.py
显示
QIIME library version: 1.9.1QIIME script version: 1.9.1
qiime-default-reference version: 0.1.3NumPy version: 1.10.4SciPy version: 0.17.1pandas version: 0.23.1matplotlib version: 1.4.3biom-format version: 2.1.6h5py version: 2.8.0 (HDF5 version: 1.10.1)qcli version: 0.1.1pyqi version: 0.3.2scikit-bio version: 0.2.3PyNAST version: 1.2.2Emperor version: 0.9.51burrito version: 0.9.1burrito-fillings version: 0.1.1sortmerna version: SortMeRNA version 2.0, 29/11/2014sumaclust version: SUMACLUST Version 1.0.00swarm version: Swarm 1.2.19 [Apr 12 2017 17:30:22]gdata: Installed.
5.创建文件夹
分为数据文件夹、参数文件夹
mkdir rawdata
mkdir config
cp /Bio/Rawdata/TGM0001-tutorial-Soil-4-meta16S/*.fastq.gz rawdata/
cp /Bio/Rawdata/TGM0001-tutorial-Soil-4-meta16S/*.txt config/
cp /Bio/Rawdata/TGM0001-tutorial-Soil-4-meta16S/*.R config/
6.检测map文件是否有错
$validate_mapping_file.py -o vmf-map/ -m config/map.txt
No errors or warnings were found in mapping file.
7.拼接
方法默认[default: fastq-join]
MIN_OVERLAP 默认为none,最小为8个bp
PERC_MAX_DIFF 错配默认为none,最大错配为10%$cat work.sh
join_paired_ends.py -m fastq-join -j 8 -p 10 -f rawdata/SRR1370913_1.fastq.gz -r rawdata/SRR1370913_2.fastq.gz -o joined/SRR1370913
join_paired_ends.py -m fastq-join -j 8 -p 10 -f rawdata/SRR1370914_1.fastq.gz -r rawdata/SRR1370914_2.fastq.gz -o joined/SRR1370914
join_paired_ends.py -m fastq-join -j 8 -p 10 -f rawdata/SRR1370915_1.fastq.gz -r rawdata/SRR1370915_2.fastq.gz -o joined/SRR1370915
join_paired_ends.py -m fastq-join -j 8 -p 10 -f rawdata/SRR1370920_1.fastq.gz -r rawdata/SRR1370920_2.fastq.gz -o joined/SRR1370920
8.split
split_libraries_fastq.py
- fastq转fasta
- 去短片段和低质量
- 去N
split_libraries_fastq.py -i joined/SRR1370913/fastqjoin.join.fastq,joined/SRR1370914/fastqjoin.join.fastq,joined/SRR1370915/fastqjoin.join.fastq,joined/SRR1370920/fastqjoin.join.fastq --sample_id SRR1370913,SRR1370914,SRR1370915,SRR1370920 -o slout/ -m config/map.txt -q 19 --barcode_type 'not-barcoded' --phred_offset=33
数据结果:
total 363M
-rw-rw-r-- 1 toucan toucan 551 Jun 14 22:28 histograms.txt
-rw-rw-r-- 1 toucan toucan 363M Jun 14 22:28 seqs.fna
-rw-rw-r-- 1 toucan toucan 2.0K Jun 14 22:28 split_library_log.txt
9.pick OTU
OTU分析与注释–pick_open_reference_otus.py
- pick_otus.py——从fasta文件中提取OTUs,提取方法有cd-hit,blast,Mothur,usearch。
- make_otu_table.py——-将提取出的OTUs制作成out_table,即.biom文件。
- assign_taxonomy.py——使用reference中的reference_seqs和 taxonomy文件对序列分类。
- biom add-metadata——-将taxa数据加入到.biom文件。
$cat params.txt
pick_otus:enable_rev_strand_match True
运行命令:
pick_open_reference_otus.py -o otus/ -i slout/seqs.fna -p config/params.txt
允许正反向聚类分析
相当于不允许的两倍
log文件记录详细的子程序脚本:
otus/log_20180614223759.txt
如100bp,测序读错了3个碱基以上,就不能满足97%相似度。出现一条的视为测序错误。
head otus/step1_otus/seqs_clusters.uc
# uclust --input /tmp/UclustExactMatchFilterM9PRwT.fasta --id 0.97 --tmpdir /tmp --rev --w 8 --stepwords 8 --usersort --maxaccepts 1 --libonly --stable_sort --maxrejects 8 --lib /home/toucan/miniconda3/envs/qiime1/lib/python2.7/site-packages/qiime_default_reference/gg_13_8_otus/rep_set/97_otus.fasta --uc otus//step1_otus/seqs_clusters.uc
qiime1自带数据库为:
/home/toucan/miniconda3/envs/qiime1/lib/python2.7/site-packages/qiime_default_reference/gg_13_8_otus/rep_set:
97_otus.fasta/home/toucan/miniconda3/envs/qiime1/lib/python2.7/site-packages/qiime_default_reference/gg_13_8_otus/rep_set_aligned:
85_otus.pynast.fasta/home/toucan/miniconda3/envs/qiime1/lib/python2.7/site-packages/qiime_default_reference/gg_13_8_otus/taxonomy:
97_otu_taxonomy.txt/home/toucan/miniconda3/envs/qiime1/lib/python2.7/site-packages/qiime_default_reference/gg_13_8_otus/trees:
97_otus.tree
后续使用的biom文件为:
otus/otu_table_mc2_w_tax_no_pynast_failures.biom
非pynast失败的,大于2条count的OTU表
10.画OTU维恩图
将OTU为0的进行舍弃组成的列表
$cat config/venn.R
#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)))
11.画OTU热图
画前100个OTU丰度值
$cat config/heatmap.R
#install.packages("gplots")
lineNum=100library('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()
12.多样性
core_diversity_analyses.py
1. single_rarefaction.py——Rarify the OTU table-测序深度统一校正为最小的样品count数
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
执行命令
core_diversity_analyses.py -o cdout/ -i otus/otu_table_mc2_w_tax_no_pynast_failures.biom -m config/map.txt -t otus/rep_set.tre -e 195965 -p config/diversityParams.txt
参数
-e为OTUcount数最小的样品数,先过来少于-e count数的样品,再统一校正为-e count的测序深度
-p为各子程序参数配置文件
$ cat config/diversityParams.txt
alpha_diversity:metrics observed_species,chao1,shannon,simpson
13.多样性说明
Alpha多样性(样本内多样性)
Alpha多样性是指一个特定区域或者生态系统内的多样性,
常用的度量指标有Chao1 丰富度估计量(Chao1 richness
estimator) 、香农多样性指数(Shannon diversity index)、
辛普森多样性指数(Simpson diversity index)等。
计算菌群丰度:Chao;
计算菌群多样性(丰度和均匀度):Shannon、Simpson。
Chao1 丰富度估计量
Chao1:是用chao1 算法估计群落中含OTU数目的指数,chao1 在生态学中常用来估计物种总数,由Chao (1984)
最早提出。Chao1值越大代表物种总数越多。
计算公式:
Schao1=Sobs+n1(n1-1)/2(n2+1)
香农和辛普森多样性指数
Shannon-Wiener 曲线,是利用shannon指数来进行绘制的,反映样品中微生物多样性的指数,利用各样品的测序量在不同测序深度时的微生物多样性指数构建曲线,以此反映各样本在不同测序数量时的微生物多样性。 当曲线趋向平坦时,说明测序数据量足够大,可以反映样品中绝大多数的微生物物种信息。
Simpson:用来估算样品中微生物的多样性指数之一,在生态学中常用来定量的描述一个区域的生物多样性。Simpson指数值越大,说明群落多样性越高。
Rank-Abundance曲线
同时解释样品所含物种的丰富程度和均匀程度。横坐标代表物种排序的数量;纵坐标代表观测到的
相对丰度。物种的丰富程度由曲线在横轴上的长度来反映,曲线越宽,表示物种的组成越丰富;物种组成的均匀程度由曲线的形状来反映,曲线越平坦,表示物种组成的均匀程度越高。而曲线快速陡然下降表明样本中的优势菌群所占比例很高,多样性较低。
Beta多样性分析(样品间差异分析)
Beta多样性是生态系之间的种多样性,它包含分类单位的比较生态系内的每一独特性。即衡量群落之间的差别。
Beta多样性不仅描述生境内生物种类的数量,同时也考虑到这些种类的相同性及其彼此之间的位置。不同生境间或某一生态梯度上不同地段间生物种类的相似性越差,Beta生物多样性越高。
基于两个样品OTU丰度之间比较;
基于系统发生树进行比较,亲缘关系;如OTU的序列碱基差异比较大,但亲缘关系比较亲,可以抵消差异。
PCoA(主坐标分析)
PCoA分析,即主坐标分析(PCoA,PrincipalCoordinate Analysis)可以从复杂的多维度数据中提取和可视化展示主要的高信息量突变,与PCA的主要思想类似。它将距离矩阵的信息转变成在一系列直交轴中显示,这样某种变量的最大值体现在第一个成分坐标中,而另一种变量的最大值体现在另一个成分坐标上。距离越近的点代表样本间相似度越高,理论上生物学重复的样本应该聚在一起。
如果三个坐标占85%以上,可以表示与数据一致。
unweighted unifrac
对于系统发生树所有枝,考查其指向的叶节点是否只存于同一个群落,那些叶节点只存在于同一群落的枝的枝长和,占整个树的枝长和的比例,就定义为UniFrac距离。
UniFrac计算了仅被一个群落占据的进化历史的相对大小,这个量越大,说明两个群落中独立的进化过程越多,也就说明这两个群落的差别越大。若两个群落完全相同,那么它们没有各自独立的进化过程,UniFrac值为0;若两个群落在进化树中完全分开,即它们是完全独立的两个进化过程,那么UniFrac值为1.
只考虑树枝是否相同
weighted unifrac
UniFrac的定义中,可以看出它只考虑序列是否在群落中出现,而不考虑序列的丰度。若两个群落包含的物种完全相同,那么不管每个物种的丰度是否有差别或者差别的大小,UniFrac值为0。但在某些具体的情况下,研究者感兴趣的恰恰是群落中物种丰度的变化,例如研究人体肠道微生物分布在抗生素治疗下的丰度变化情况,这时UniFrac就不能解决问题了。
weighted unifrac方法,就是在UniFrac的基础上,将序列的丰度纳入考虑,它能够区分物种丰度的差别。
14.系统树
系统发生树(英文:phylogenetic tree或evolutionary tree)是表明被认为具有共同祖先的各物种相互间演化关系的树,又被称为系统发育树、系统演化树、系统进化树、种系发生树、演化树、进化树、系统树。 它用来表示系统发生研究的结果,用它描述物种之间的进化关系。
构建系统树主要步骤:
1. filter_otus_from_otu_table.py——挑选丰度>1‰的OTU
2. filter_fasta.py——获得丰度>1‰OTU的fasta文件
3. make_phylogeny.py——-构建系统树
4. ggtree.R——系统树可视化
执行脚本:
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
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
make_phylogeny.py -i otus/rep_set_0.001_fraction.fasta -o otus/rep_set_0.001_fraction.tre
biom convert -i otus/otu_table_mc2_0.001_fraction.biom -o otus/table.from_0.001_fraction.txt --to-tsv --header-key taxonomy
R脚本系统发生树
$cat config/ggtree.R
#source("https://bioconductor.org/biocLite.R") #source("http://bioconductor.org/biocLite.R")
#biocLite("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 = "")
# 切割门物种,将out传进门的容器列表中
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)# 读取树文件,确定,按门的名字进行分组otu
tree <- read.tree("rep_set_0.001_fraction.tre")
tree <- groupOTU(tree, cls)# 去除NewReference,只显示OTU号
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)# 最后一列为tax,需要-1,取log值标准化,画热图。
gplot(p, log(rt[,1:ncol(rt)-1]+1), font.size=4)
dev.off()
各个门的分布
15.差异OTU分析
根据分组计算log2FoldChange值,计算P值和校正后P值,然后使用校正后P值进行筛选少于0.05的物种。
-a 比较方法
-c 分组类别
-x -y 只限于两组之间比较
#OTU Differential Abundance Testing
differential_abundance.py -i otus/otu_table_mc2_w_tax_no_pynast_failures.biom -o diff_otus.txt -m config/map.txt -c SampleType -x normal -y tumor -a DESeq2_nbinom
注:conda 安装的qiime1缺乏,optparse包,需要额外安装 conda install -n qiime1 optparse
报错显示:To install, open R and run the command: install.packages("optparse")
0044-【宏基因组】-16S分析qiime1极简教程相关推荐
- Nature综述: 宏基因组关联分析-深入研究微生物组
本文由谢忠杰编译,董小橙.江舜尧编辑,本文较长,建议用电脑阅读. "微生太"原创微文,转载已获授权. 导读 问题1:哪些疾病与人体微生物明确相关? 问题2:如何研究人体微生物与健康 ...
- 来来来,一起来pick宏基因组binning分析工具
发表期刊:Computational and Structural Biotechnology Journal(IF=7.271) 发表时间:2021 研究背景 微生物是生物和环境的营养循环和代谢过程 ...
- 易基因|多组学关联研究怎么做? DNA甲基化组+转录组+宏基因组+16S研究思路
大家好,这里是专注表观组学十余年,领跑多组学科研服务的易基因. 本期我们以一篇多组学研究揭示DNA甲基化在植物促生菌-植物-微生物互作关系中作用的文献案例来讲解多组学研究如何做,并对多组学研究分析方法 ...
- 负载分析及问题排查极简教程
作者 | Hollis ,来自 | Hollis 平常的工作中,在衡量服务器的性能时,经常会涉及到几个指标,load.cpu.mem.qps.rt等.每个指标都有其独特的意义,很多时候在线上出现问题时 ...
- 高效sql性能优化极简教程
一,sql性能优化基础方法论 对于功能,我们可能知道必须改进什么:但对于性能问题,有时我们可能无从下手.其实,任何计算机应用系统最终队可以归结为: cpu消耗 内存使用 对磁盘,网络或其他I/O设备的 ...
- CentOS安装使用.netcore极简教程(免费提供学习服务器)
本文目标是指引从未使用过Linux的.Neter,如何在CentOS7上安装.Net Core环境,以及部署.Net Core应用. 仅针对CentOS,其它Linux系统类似,命令环节稍加调整: 需 ...
- 《Kotlin 极简教程 》第5章 集合类
<Kotlin 极简教程 >第5章 集合类 <Kotlin极简教程>正式上架: 点击这里 > 去京东商城购买阅读 点击这里 > 去天猫商城购买阅读 非常感谢您亲爱的 ...
- 《Kotin 极简教程》第13章 使用 Kotlin 和 Anko 的Android 开发
第13章 使用 Kotlin 和 Anko 的Android 开发 最新上架!!!< Kotlin极简教程> 陈光剑 (机械工业出版社) 可直接打开京东,淘宝,当当===> 搜索: ...
- 【“计算机科学与技术”专业小白成长系列】Linux Shell 编程 极简教程
Linux Shell 编程 极简教程 内容摘要 本文是 Linux Shell 编程简单入门.主要内容: Linux 简介 Shell 编程入门 Kotlin 脚本与 Shell 脚本 Linux ...
最新文章
- vc++向txt文件中写入数据,追加数据
- linux autofs ftp,Linux NFS自动挂载autofs配置
- [SDOI2016]生成魔咒
- Python的捕虫笔记
- C#代码创建3D模型
- appium+python 操作APP
- 读书笔记1 : program paradigm
- 爱因斯坦《我的世界观》
- MyBatis 源码解读-mapperElement()
- 用POP动画引擎实现弹簧动画(POPSpringAnimation)
- 利用JDK1.5的工具对远程的Java应用程序进行监测(摘录)
- Codeforces 678E:Another Sith Tournament 状压DP
- 【二叉树】二叉树遍历总结
- java的小于等于符号怎么打_「小于符号」mybatis的一些特殊符号标识(大于,小于,等于,不等于) - seo实验室...
- 技术分享 | show engine innodb status中Pages flushed up to 的含义
- 冒泡排序算法,简单明了哦。
- if lte IE if gte IE 浏览器兼容
- 云栖重磅!阿里云启动视频云V5计划,全面赋能生态合作伙伴
- 数据库 存储过程的建立 调用 加密
- Spring Boot 实现接口的各种参数校验
热门文章
- MATLAB打开后一直在初始化,或者初始化很慢问题
- 海龟交易法则(中译文)
- Pygame实战:家里的小孩数学算数能力很差嘛?别慌—这款“巧算24点小游戏”等你来玩,管用。
- linux系统漏洞升级方法,OpenSSL “Heartbleed”心脏流血漏洞升级方法
- 体脂秤方案开发脂肪秤方案设计
- python随机生成英文字符串_如何用Python语言生成随机字符串 | 学步园
- Neos.Flow UnitTestBootstrap
- linux脚本循环创建用户,shell应用之批量添加用户实例
- LSTM和GRU的对比和分析
- 计算机课情感态度与价值观,浅谈信息技术课中情感态度价值观的培养