微生物组分析一：稀释曲线和α多样性

简介

微生物组分析是一种研究生态系统中微生物种群组成、功能和相互作用的重要手段。通过对环境样本（如土壤、水体、生物组织）中的微生物DNA或RNA进行高通量测序，可以获得微生物的种类和丰度信息。这些数据能够帮助我们揭示微生物在生态系统中的作用及其对环境变化的响应。

从今天开始，生信分析实习生将陆续分析如何使用R语言进行微生物组数据（16S、ITS、宏基因组数据）基础分析。今天分析稀释曲线和α多样性的分析内容，后台私信发送“微生物组分析一”可领取示例数据与源代码。

step1 数据详情、读取以及预处理

本分析需要三个文件，otu_table.tsv是物种丰度文件，第一列是ASV ID，后续每一列为各样本的物种丰度数据。

group2.txt是分组数据，第一列是样本名，第二列是对应的分组。

taxonomy.tsv则是ASV的注释信息，第一列Feature ID就是第一个文件中的ASV ID一致。第二列Taxon则是ASV的注释信息，最后一列是注释结果的置信度，这一列在后续分析中会剔除掉。

大家可以把自己的数据整理成这个样子，然后套用以下代码就能进行分析。需要注意的是，不一定要把数据整理成tsv格式，可以用csv, txt, xls格式都行，只不过在读取数据的时候需要根据数据格式改变读取方法。tsv, txt，xls格式都可以用read.table()函数读取，csv需要使用read.csv()。

其次需要注意的是，taxonomy.tsv文件中的Feature ID列必须要与otu_table.tsv文件中的ASV ID一致。

library(phyloseq)
library(MicrobiotaProcess)
library(tidyverse)

##建立输出目录
dir.create('./alpha_diversity')
dir.create('./rarecurve')

##读取数据文件
otu_txt <– 'phyloseq/otu_table.tsv'
sample <– "phyloseq/group2.txt"
taxon_txt <– 'phyloseq/taxonomy.tsv'

otu_data <– read.table(otu_txt,header = T,row.names = 1,sep = '\\t')
group_info <– read.table(sample,header = T,row.names = 1,sep = '\\t')

taxon_data <– read.table(taxon_txt,header = T,sep = '\\t',row.names = 1,check.names = F)
taxon_data <– taxon_data[rownames(taxon_data) %in% rownames(otu_data),][–2]

tax_table <– taxon_data %>%
mutate(Kingdom = str_extract(Taxon, "d__[^;]+"),
Phylum = str_extract(Taxon, "p__[^;]+"),
Class = str_extract(Taxon, "c__[^;]+"),
Order = str_extract(Taxon, "o__[^;]+"),
Family = str_extract(Taxon, "f__[^;]+"),
Genus = str_extract(Taxon, "g__[^;]+"),
Species = str_extract(Taxon, "s__[^;]+")) %>%
select(–Taxon)

step2 构建Phyloseq对象进行分析

# 创建 OTU 表格
otu_table <– phyloseq::otu_table(as.matrix(otu_data), taxa_are_rows = TRUE)

# 创建样本数据
sample_data <– phyloseq::sample_data(group_info)

# 创建分类学表格
taxonomy_table <– phyloseq::tax_table(as.matrix(tax_table))

# 创建 phyloseq 对象
physeq <– phyloseq(otu_table, sample_data, taxonomy_table)

step3 抽平

为了保证测序序列的均一性，必须对数据进行抽平处理。可以看到，191个物种被剔除

set.seed(1234)

###如果没有抽平就进行以下步骤
rare.data <– rarefy_even_depth(physeq,replace = T)
mpse2 <– rare.data

step4 α多样性组合图

α多样性是描述单个样本中微生物群落多样性的指标，通常包括物种丰富度（Species Richness）、物种均匀度（Species Evenness）等指标。它反映了样本内微生物种群的复杂性。

## Alpha多样性组合图
mouse.time.mpse <– as.MPSE(mpse2)

mouse.time.mpse %<>%
mp_cal_alpha(.abundance=RareAbundance)

mycol <– c("#B54959", "#2066AC",'lightblue','#FCB2AF','#9BDFDF','#FFE2CE','#C4D8E9')

#alpha diversity
f1 <– mouse.time.mpse %>%
mp_plot_alpha(
.group=group,
test = "wilcox.test",
step_increase = 0.15,
.alpha=c(Observe, Chao1, ACE, Shannon, Simpson, Pielou)
) +
scale_fill_manual(values=mycol, guide="none") +
scale_color_manual(values=mycol, guide="none")

f2 <– mouse.time.mpse %>%
mp_plot_alpha(
.alpha=c(Observe, Chao1, ACE, Shannon, Simpson, Pielou)
)

f=f1 / f2

ggsave(f,filename = './alpha_diversity/alpha_diversity.png',dpi = 300,
width = 10,height = 8)
ggsave(f,filename = './alpha_diversity/alpha_diversity.pdf',
width = 10,height = 8)

step5 稀释曲线和单指标α多样性

## 稀释曲线绘制、Alpha多样性单指标可视化
ps_dada2 <- rare.data
ps_index <- get_alphaindex(ps_dada2)

for (i in colnames(ps_index@alpha)) {
  ## 稀释曲线
  p_rare <- ggrarecurve(obj=ps_dada2,
                        indexNames=i,
                        chunks=100) +
    theme_bw()+
    theme(
      panel.grid = element_blank(),
      strip.background=element_blank(),
      strip.text=element_blank(),
      legend.spacing.y=unit(0.02,"cm"),
      axis.title.y= element_text(size=20, color='black',face='bold'),
      axis.title.x= element_text(size=20, color='black',face='bold'),
      legend.text=element_text(size=15),
      legend.title = element_text(size = 15,face = 'bold'))+
    labs(y=paste0(i),color='')
  ggsave(filename = paste0('./rarecurve/',i,'.png'),dpi = 300,
         width = 8,height = 8)
  ggsave(filename = paste0('./rarecurve/',i,'.pdf'),
         width = 8,height = 8)

  ## Alpha多样性单指标
  p_alpha <- ggbox(ps_index, geom = "boxplot",
                   factorNames="group",
                   indexNames=i,
                   p_textsize = 4,
                   step_increase=0.15,
                   #boxwidth=0.1,
                   signifmap = F,
                   testmethod='wilcox.test') +
    scale_fill_manual(values=mycol)+
    geom_jitter(size=1.5,color="black") +
    theme_bw()+
    theme(
      strip.background=element_blank(),
      strip.text=element_blank(),
      plot.margin= unit(c(10, 10, 5, 10), 'mm'),
      legend.position             = "none",
      panel.grid = element_blank(),
      panel.background    = element_rect(fill='transparent'),
      panel.border        = element_rect(fill="transparent",color="black",size=1),
      axis.text.x         = element_text(angle = 0, hjust = 0.5, size=15, color='black'),
      axis.text.y         = element_text(size=15, color='black'),
      axis.title.x                = element_text(size=10, color='black', face='bold'),
      axis.title.y                = element_text(size=20, color='black', vjust=3),
      axis.ticks.x                = element_line(size = 0, color='white'),
      axis.ticks.y                = element_line(size = 1, color='black'),
      axis.ticks.length.x = unit(0.01, "cm"),
      axis.ticks.length.y = unit(0.01, "cm")
    )+
    labs(y=paste0(i))
  ggsave(filename = paste0('./alpha_diversity/',i,'.pdf'),plot = p_alpha,width = 10,height = 7)
  ggsave(filename = paste0('./alpha_diversity/',i,'.png'),dpi = 300,
         plot = p_alpha,
         width = 10,height = 7)
}

step6 α多样性表格整理

## Alpha_diversity指标表格整理
outdir='./alpha_diversity'
write.table(ps_index, file=paste(outdir,"/alpha_diversity.tmp",sep=""),sep=",",quote=FALSE,col.names=FALSE)
alpha_result<-read.table(paste(outdir,"/alpha_diversity.tmp",sep=""),header=F,sep=",")
colnames(alpha_result)<-c("Sample","Observe","Chao1","ACE","Shannon","Simpson","Pielou","Group")
write.table(alpha_result, file=paste(outdir,"/alpha_diversity.xls",sep=""),sep="\\t",quote=FALSE,row.names=FALSE)

alpha_stat <– data.frame()
for (group1 in unique(alpha_result$Group)) {
alpha_stat_group1 <– c(group1,
paste(round(mean(alpha_result$Observe[alpha_result$Group==group1]),4),round(sd(alpha_result$Observe[alpha_result$Group==group1]),4),sep=" +/- "),
paste(round(mean(alpha_result$ACE[alpha_result$Group==group1]),4),round(sd(alpha_result$ACE[alpha_result$Group==group1]),4),sep=" +/- "),
paste(round(mean(alpha_result$Chao1[alpha_result$Group==group1]),4),round(sd(alpha_result$Chao1[alpha_result$Group==group1]),4),sep=" +/- "),
paste(round(mean(alpha_result$Shannon[alpha_result$Group==group1]),4),round(sd(alpha_result$Shannon[alpha_result$Group==group1]),4),sep=" +/- "),
paste(round(mean(alpha_result$Simpson[alpha_result$Group==group1]),4),round(sd(alpha_result$Simpson[alpha_result$Group==group1]),4),sep=" +/- "),
paste(round(mean(alpha_result$Pielou[alpha_result$Group==group1]),4),round(sd(alpha_result$Pielou[alpha_result$Group==group1]),4),sep=" +/- ")
)
alpha_stat<-rbind(alpha_stat,alpha_stat_group1)
}

colnames(alpha_stat)<-c("Group","Observe","ACE","Chao1","Shannon","Simpson",'Pielou')
write.table(alpha_stat, file=paste(outdir,"/alpha_diversity.stat.xls",sep=""),sep="\\t",quote=FALSE,row.names=FALSE)