其他
助力微生物组学分析--实验室新R包MicrobiotaProcess
1. MicrobiotaProcess
MicrobiotaProcess是本实验室最近开发的一款用于微生物组学分析的R包。这个包的功能主要有三部分构成，第一，最新的聚类软件qiime2与dada2软件的输出数据的导入。第二，常见的微生物组学数据分析，如alpha多样性，稀释曲线分析，物种组成分析,PCA,PCoA与层次聚类分析。第三，也是这款软件的重要功能，就是biomarker查找。目前的结果显示MicrobiotaProcess可以控制极低的假阳性。
2. 实例分析
2.1. 数据导入
数据导入MicrobiotaProcess提供import_dada2与import_qiime2两个函数，支持dada2以及qiime2输出结果的导入。
library(MicrobiotaProcess)
library(phyloseq)
seqtabfile <- system.file("extdata", "seqtab.nochim.rds", package="MicrobiotaProcess")
# seqtabfile是dada2的removeBimeraDenovo输出结果
seqtab <- readRDS(seqtabfile)
taxafile <- system.file("extdata", "taxa_tab.rds",package="MicrobiotaProcess")
# taxafile是dada2的assignTaxonomy输出结果
seqtab <- readRDS(seqtabfile)
taxa <- readRDS(taxafile)
# the seqtab and taxa are output of dada2
sampleda <- system.file("extdata", "mouse.time.dada2.txt", package="MicrobiotaProcess")
ps_dada2 <- import_dada2(seqtab=seqtab, taxatab=taxa, sampleda=sampleda)
ps_dada2
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 232 taxa and 19 samples ]
## sample_data() Sample Data: [ 19 samples by 2 sample variables ]
## tax_table() Taxonomy Table: [ 232 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 232 reference sequences ]
# import data from qiime2 pipeline
# table.qza 与 taxa.qza均为qiime2的输出结果
otuqzafile <- system.file("extdata", "table.qza", package="MicrobiotaProcess")
taxaqzafile <- system.file("extdata", "taxa.qza", package="MicrobiotaProcess")
mapfile <- system.file("extdata", "metadata_qza.txt", package="MicrobiotaProcess")
ps_qiime2 <- import_qiime2(otuqza=otuqzafile, taxaqza=taxaqzafile, mapfilename=mapfile)
ps_qiime2
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 138 taxa and 87 samples ]
## sample_data() Sample Data: [ 87 samples by 23 sample variables ]
## tax_table() Taxonomy Table: [ 138 taxa by 8 taxonomic ranks ]
## refseq() DNAStringSet: [ 138 reference sequences ]
导入结果为phyloseq对象。当然，如果你有phyloseq对象或者物种丰度表，都是可以用MicrobiotaProcess进行下游分析的。
2.2. 稀释曲线分析与alpha多样性比较
稀释曲线分析与多样性比较是微生物组学比较常用的分析之一，前者可用于评估测序深度是否充足，后者可以比较不同群落间的物种均匀度以及丰富度。
# for reproducibly random number
library(ggplot2)
set.seed(1024)
p_rare <- ggrarecurve(obj=ps_dada2,
indexNames=c("Observe","Chao1","ACE"),
chunks=300) +
theme(legend.spacing.y=unit(0.02,"cm"),
legend.text=element_text(size=6))
p_rare
# 首先计算各个样品的alpha多样性值。
set.seed(1024)
alphaobj <- get_alphaindex(ps_dada2)
head(as.data.frame(alphaobj))
## Observe Chao1 ACE Shannon Simpson J sample time
## F3D0 105 105.9091 106.56040 3.842245 0.9635779 0.8255860 F3D0 Early
## F3D1 99 100.5000 100.12294 3.982718 0.9701640 0.8667277 F3D1 Early
## F3D141 73 73.0000 73.00000 3.419822 0.9499830 0.7970760 F3D141 Late
## F3D142 48 48.0000 48.00000 3.118059 0.9386697 0.8054500 F3D142 Late
## F3D143 56 56.0000 56.00000 3.292717 0.9464422 0.8179949 F3D143 Late
## F3D144 47 47.0000 47.17407 2.993200 0.9315811 0.7774246 F3D144 Late
# 然后再用ggbox来进行可视化
p_alpha <- ggbox(alphaobj, geom="violin", factorNames="time") +
scale_fill_manual(values=c("#2874C5", "#EABF00"))+
theme(strip.background = element_rect(colour=NA, fill="grey"))
p_alpha
稀释曲线是对每个样品的测序reads随机抽取一定的reads数，然后计算该reads数下的多样性。chunks就是要对各个样品分成多少份数抽取。alpha多样性比较默认使用wilcox.test进行两两差异分析。
2.3 PCA, PCoA以及层次聚类分析
MicrobiotaProcess软件包提供pca,pcoa计算，同时支持多种距离计算，常见的bray-curtis,(un)weighted unifrac均支持。
pcares <- get_pca(obj=ps_dada2, method="hellinger")
# Visulizing the result
pcaplot <- ggordpoint(obj=pcares, biplot=TRUE, speciesannot=TRUE,
factorNames=c("time"), ellipse=TRUE) +
scale_color_manual(values=c("#2874C5", "#EABF00"))
pcaplot
pcoares <- get_pcoa(obj=ps_dada2, distmethod="euclidean", method="hellinger")
# Visualizing the result
pcoaplot <- ggordpoint(obj=pcoares, biplot=TRUE, speciesannot=TRUE,
factorNames=c("time"), ellipse=TRUE) +
scale_color_manual(values=c("#2874C5", "#EABF00"))
pcoaplot
hcsample <- get_clust(obj=ps_dada2, distmethod="euclidean",
method="hellinger", hclustmethod="average")
# rectangular layout
library(ggtree)
clustplot1 <- ggclust(obj=hcsample,
layout = "rectangular",
pointsize=1,
fontsize=0,
factorNames=c("time")) +
scale_color_manual(values=c("#2874C5", "#EABF00")) +
theme_tree2(legend.position="right",
plot.title = element_text(face="bold", lineheight=25,hjust=0.5))
clustplot1
2.3 物种组成分析
MicrobiotaProcess提供物种分类水平提取以及可视化。
# 提取各个分类物种水平的数据然后进行可视化。
phytax <- get_taxadf(obj=ps_dada2, taxlevel=2)
phybar <- ggbartax(obj=phytax) +
xlab(NULL) + ylab("relative abundance (%)")
phybar
#绝对丰度， facetNames控制分组分面。
phybar2 <- ggbartax(obj=phytax, facetNames="time", count=TRUE) +
xlab(NULL) + ylab("abundance")
phybar2
2.4 biomarker鉴定
MicrobiotaProcess提供了一个灵活可控的功能来进行biomarker鉴定。大体原理与LEfSe相似，具体原理将在以后更新，目前该方法的假阳性控制很低。
#这里我们用发表在GenomeResearch上的CRC的研究数据为例。
#
#`kostic2012crc`是一个phyloseq对象
data(kostic2012crc)
kostic2012crc <- phyloseq::rarefy_even_depth(kostic2012crc,rngseed=1024)
set.seed(1024)
diffres <- diff_analysis(obj=kostic2012crc, class="DIAGNOSIS",
mlfun="lda",
filtermod="fdr",
firstcomfun = "kruskal.test",
firstalpha=0.05,
strictmod=TRUE,
secondcomfun = "wilcox.test",
subclmin=3,
subclwilc=TRUE,
secondalpha=0.01,
lda=3)
diffres
## The original data: 693 features and 177 samples
## The sample data: 1 variables and 177 samples
## The taxda contained 1351 by 7 rank
## after first test (kruskal.test) number of feature (fdr<=0.05):103
## after second test (wilcox.test) number of significantly discriminative feature:40
## after lda, Number of discriminative features: 36 (certain taxonomy classification:26; uncertain taxonomy classication: 10)
该结果是一个S4对象。我们提供了ggdiffbox,ggdiffclade,ggeffectsize以及ggdifftaxbar来对这个结果进行可视化。
plotes_ab <- ggdiffbox(obj=diffres, box_notch=FALSE, colorlist=c("#00AED7", "#FD9347"))
plotes_ab
左边为差异物种在不同分组中的丰度盒形图，右边为该差异物种的效应值的置信区间，默认为95%。结果发现了与原文章一致的Fusobacterium属在肿瘤组中显著富集，还发现Campylobacter也可能在肿瘤中显著富集，但该属效应值的置信区间明显波动较大，说明受到样品影响可能较大。
diffcladeplot <- ggdiffclade(obj=diffres,
layout="radial",
alpha=0.3, size=0.2,
skpointsize=0.6,
taxlevel=3,
settheme=FALSE,
setColors=FALSE) +
scale_fill_manual(values=c("#00AED7", "#FD9347"))+
guides(color = guide_legend(keywidth = 0.1,
keyheight = 0.6,
order = 3,
ncol=1)) +
theme(panel.background=element_rect(fill=NA),
legend.position="right",
plot.margin=margin(0,0,0,0),
legend.spacing.y = unit(0.02, "cm"),
legend.title=element_text(size=7),
legend.text=element_text(size=6),
legend.box.spacing=unit(0.02,"cm"))
diffcladeplot
差异物种分类层次树形图，由内到外分别为界门纲目科属种(取决于输入数据分类情况)，点的大小为-log10(p),p值是第一次比较的显著水平。
软件安装地址
目前MicrobiotaProcess正在bioconductor中排队review，想要试用的用户可以通过github安装,
需要注意的是目前本软件得在R-4.0.0下才能安装。不同版本的R软件环境可以参考《为了提交R包，我构建了个多版本的R开发环境》。
if (! requireNamespace("remotes", quietly=TRUE)){
install.packages("remotes")
}
remotes::install_github("YuLab-SMU/MicrobiotaProcess")
更多文档可以安装完成用如下操作进行查阅。
browseVignettes("MicrobiotaProcess")
问题反馈
用户使用过程中如果发现问题，可以到https://github.com/YuLab-SMU/MicrobiotaProcess/issues提交问题，
同时可以关注本公众号，本公众号将会不定时介绍该包的使用。