宏基因组 - (2)基于基因相对丰度计算α/β多样性
0. 首先需要获取基因的相对丰度表，这里可以略参考 宏基因组-(1)基因预测与基因相对丰度
图1. 大概是这样 ，个人比较喜欢csv格式，第一行是样品品，第一列是基因名。其实和OTU表差不多
1. R的读取以及处理
因为该文件有600M，所以用基础函数有点吃不消，当然也可以用python
以下是R的代码
library(data.table)
data <- fread("./genemark.GeneAbun.csv",stringsAsFactors = F,encoding = "UTF-8") #genemark.GeneAbun.csv即基因相对丰度表
tdata = t(data) #行列转置
a <- as.data.frame(tdata) #转为数据框
colnames(a) <- a[1,] #把第一行作为列名，令人奇怪的是，第一列已经自动作为行名了，所以后续不管了
a = a[-1,] #第一行作为行名，但是第一行还是和行名一模一样，所以去除第一行
a = as.data.frame(lapply(a,as.numeric)) #把里面的 数字从字符串转为数值型(numeric)
到这，文件就处理完毕了，可以用以后续的计算
2.分组信息的处理
group=c(“B”,"B","B","B","B","B","C","C","C","C","C","C") #12个样品，两大组，这里的顺序要和数据里面的顺序排得上,可以看到图1中是前6个样品是B大组，后6个是C大组
sample_id = rownames(a) #行列已经转置了，把第一列列名作为个体名，这里的顺序依然要和图1中
data_group = data.frame(sample_id, group) #构建一个分组信息文件
3.α多样性 - shannon指数
library(vegan)
library(reshape)
geneshannon = diversity(a, "shannon")
data_ggplot=data.frame(geneshannon)
data_ggplot=data.frame(data_ggplot, data_group["group"])
group_C=data_ggplot$geneshannon[data_ggplot$group=='C']
group_B=data_ggplot$geneshannon[data_ggplot$group=='B'] #利用分组信息
#检验是否满足正态性
shapiro.test(group_C) shapiro.test(group_B)
若p值小于0.05，可以认为不满足正态性。。我做出来p都很高
#检验是否方差齐
bartlett.test(sppshannon~group,data=data_ggplot)
若p值小于0.05，可以认为方差不齐
#满足正态，方差不齐，两样本独立且两组样本数相同，故使用 校正t检验 -- t.test(x, y, var.equal = F)
with(data_ggplot, t.test(formula=sppshannon~group,var.equal = F))
library(ggplot2) library(ggsignif)
alpha_boxplot=ggplot(data_ggplot,aes(x=group,y=geneshannon,fill=group))+
geom_boxplot()+ #箱线图
labs(title="Alpha diversity", x="Group", y="Shannon index")+ # 坐标轴标题啥的
scale_fill_manual(values = c("#70C1B3","#397AA0"))+ #挑了两个喜欢的颜色
theme(plot.title=element_text(hjust=0.5), legend.title=element_blank())+ #去除背景网格
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ #去除背景网格
geom_signif(comparisons = list(c("B","C")), test = t.test, #添加显著性标识，ggsignif实现，甚至都不用单独做检验
test.args = list(var.equal = F))
B组明显比C组高
4. Simpson指数
genesimpson = diversity(a, "simpson") #后面的处理和shannon指数一模一样，略
5. Richness指数
这个很好理解，就是单纯的数量，最后每个样品比对上了多少个基因(不是多少种基因)
这里用相对丰度表(长度标准化)肯定做不了，因为求和加起来都是1
所以用数量表就行了，这里略，毕竟相对丰度表都算出来了
之前算物种richness用的数据，处理后的，大概就是这样，还是和OTU表基本一样，下面贴的是做物种richness的代码
observed_species = rowSums(tspecies_num_raw) #计算Richness，其实就是简单的对行求和
rich_ggplot=data.frame(observed_species) #计算Richness，其实就是简单的对行求和，转为数据框格式
rich_ggplot=data.frame(rich_ggplot, rich_group["group"]) #加入分组信息
rich_Chalk=rich_ggplot$observed_species[rich_ggplot$group=='Chalk']
rich_Basalt=rich_ggplot$observed_species[rich_ggplot$group=='Basalt']
也不用照抄，只要算出来这些指数以后，怎么画图方便怎么来
6.β多样性 - BrayCurtis距离
#获取距离矩阵
bray = vegdist(a,method = "bray")
bray <- as.matrix(bray)
write.table(bray,"bray-crutis.txt",sep = "\t") #这两行是把第一列第一行作为列名行名
dis_mat = read.table("bray-crutis.txt",header = T,row.names = 1) #其实用colnames.rownames更好
#基于距离矩阵，做主坐标分析
#pcoa ref https://www.bilibili.com/read/cv10516480
dis_mat = read.table("bray-crutis.txt",header = T,row.names = 1)
pcoa = cmdscale(dis_mat,k=3,eig = T) #这里的K=3就是只算3个PcoA,最多只有11个，我猜原因是我有12个样本
poi = pcoa$points
eigval = pcoa$eig
pcoa_eig = (pcoa$eig)[]/sum(pcoa$eig) #每个PCoA的权重值，这里对应下文ggplot中的lab
poi = as.data.frame(poi)
#将poi导出到excel略作处理再导回来，主要是添加一列Group分组和把第一列加了个sample
write.csv(poi,"poi1.csv")
poi2 = read.csv("poi1.csv")
将poi导出，还没有做处理，长这样
加了一列分组的Group，第一列手动加了个sample做列名,v1即主坐标1，v2主坐标2...
pcoa_plot = ggplot(data = poi2,aes(x=V1,y=V2))+
geom_point(aes(color=group,shape=group),size=3)+
scale_colour_manual(values = c("#70C1B3","#397AA0"))+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
geom_hline(aes(yintercept=0), colour="#BEBEBE", linetype="dashed")+
geom_vline(aes(xintercept=0), colour="#BEBEBE", linetype="dashed")+
labs(x="PCoA1 ",
y="PCoA2")
权重值，从大到小即PCoA1,PCoA2......
画出来就是这样的
7.
摘自一篇文献的物种注释部分
图中用的降维算法是NMDS，和PCoA差不多，只是原理不一样
现在还没有实现的是红线标注的图，即比较组间差异和组内差异，图中的意思是组内的差异小于组间的差异，更一步说明了这两组的物种组成不一样
2021.12.17更新：https://mp.weixin.qq.com/s/sXlN4m7L72rvvKXGp2WAXA 红色标注部分参考这个链接，照抄就能做出来，此处略
8.
物种多样性的计算同理
OTU/ASV表同理
