-1.僅作筆記使用

我所用到的所有腳本可以+Q索要

請(qǐng)實(shí)名，名字+學(xué)校+方向

0. 首先需要獲取基因的相對(duì)豐度表，這里可以略參考 宏基因組-(1)基因預(yù)測(cè)與基因相對(duì)豐度

圖1. 大概是這樣 ，個(gè)人比較喜歡csv格式，第一行是樣品品，第一列是基因名。其實(shí)和OTU表差不多

1. R的讀取以及處理

因?yàn)樵撐募?00M，所以用基礎(chǔ)函數(shù)有點(diǎn)吃不消，當(dāng)然也可以用python

以下是R的代碼

library(data.table)

data <- fread("./genemark.GeneAbun.csv",stringsAsFactors = F,encoding = "UTF-8")? ? ? ? #genemark.GeneAbun.csv即基因相對(duì)豐度表

tdata = t(data)? ? ? ? ? ? ? #行列轉(zhuǎn)置

a <- as.data.frame(tdata)? ? ? #轉(zhuǎn)為數(shù)據(jù)框

colnames(a) <- a[1,]? ? ? #把第一行作為列名，令人奇怪的是，第一列已經(jīng)自動(dòng)作為行名了，所以后續(xù)不管了

a = a[-1,]? ? #第一行作為行名，但是第一行還是和行名一模一樣，所以去除第一行

a = as.data.frame(lapply(a,as.numeric))? #把里面的 數(shù)字從字符串轉(zhuǎn)為數(shù)值型(numeric)

到這，文件就處理完畢了，可以用以后續(xù)的計(jì)算

2.分組信息的處理

group=c(“B”,"B","B","B","B","B","C","C","C","C","C","C")? ? #12個(gè)樣品，兩大組，這里的順序要和數(shù)據(jù)里面的順序排得上,可以看到圖1中是前6個(gè)樣品是B大組，后6個(gè)是C大組

sample_id = rownames(a)? #行列已經(jīng)轉(zhuǎn)置了，把第一列列名作為個(gè)體名，這里的順序依然要和圖1中

data_group = data.frame(sample_id, group)? #構(gòu)建一個(gè)分組信息文件

3.α多樣性 - shannon指數(shù)

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']? #利用分組信息

#檢驗(yàn)是否滿(mǎn)足正態(tài)性

shapiro.test(group_C)? ? ? shapiro.test(group_B)

若p值小于0.05，可以認(rèn)為不滿(mǎn)足正態(tài)性。。我做出來(lái)p都很高

#檢驗(yàn)是否方差齊

bartlett.test(sppshannon~group,data=data_ggplot)

若p值小于0.05，可以認(rèn)為方差不齊

#滿(mǎn)足正態(tài)，方差不齊，兩樣本獨(dú)立且兩組樣本數(shù)相同，故使用? 校正t檢驗(yàn)? -- 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()+? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? #箱線(xiàn)圖

? labs(title="Alpha diversity", x="Group", y="Shannon index")+? ? ? ? ? ? ? ? ? # 坐標(biāo)軸標(biāo)題啥的

? scale_fill_manual(values = c("#70C1B3","#397AA0"))+? ? ? ? ? ? ? ? ? ? ? ? ? ? #挑了兩個(gè)喜歡的顏色

? theme(plot.title=element_text(hjust=0.5), legend.title=element_blank())+? #去除背景網(wǎng)格

? theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+? ? ? #去除背景網(wǎng)格

? geom_signif(comparisons = list(c("B","C")), test = t.test,? ? ? ? ? ? ? ? ? ? ? #添加顯著性標(biāo)識(shí)，ggsignif實(shí)現(xiàn)，甚至都不用單獨(dú)做檢驗(yàn)

? ? ? ? ? ? ? test.args = list(var.equal = F))

B組明顯比C組高

4. Simpson指數(shù)

genesimpson = diversity(a, "simpson")? ? #后面的處理和shannon指數(shù)一模一樣，略

5. Richness指數(shù)

這個(gè)很好理解，就是單純的數(shù)量，最后每個(gè)樣品比對(duì)上了多少個(gè)基因(不是多少種基因)

這里用相對(duì)豐度表(長(zhǎng)度標(biāo)準(zhǔn)化)肯定做不了，因?yàn)榍蠛图悠饋?lái)都是1

所以用數(shù)量表就行了，這里略，畢竟相對(duì)豐度表都算出來(lái)了

之前算物種richness用的數(shù)據(jù)，處理后的，大概就是這樣，還是和OTU表基本一樣，下面貼的是做物種richness的代碼

observed_species = rowSums(tspecies_num_raw)? ? ? #計(jì)算Richness，其實(shí)就是簡(jiǎn)單的對(duì)行求和

rich_ggplot=data.frame(observed_species)? ? ? ? ? ? ? ? #計(jì)算Richness，其實(shí)就是簡(jiǎn)單的對(duì)行求和，轉(zhuǎn)為數(shù)據(jù)框格式

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']

也不用照抄，只要算出來(lái)這些指數(shù)以后，怎么畫(huà)圖方便怎么來(lái)

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)? ? ? ? ? #其實(shí)用colnames.rownames更好

#基于距離矩陣，做主坐標(biāo)分析

#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個(gè)PcoA,最多只有11個(gè)，我猜原因是我有12個(gè)樣本

poi = pcoa$points

eigval = pcoa$eig

pcoa_eig = (pcoa$eig)[]/sum(pcoa$eig)? ? #每個(gè)PCoA的權(quán)重值，這里對(duì)應(yīng)下文ggplot中的lab

poi = as.data.frame(poi)

#將poi導(dǎo)出到excel略作處理再導(dǎo)回來(lái)，主要是添加一列Group分組和把第一列加了個(gè)sample

write.csv(poi,"poi1.csv")

poi2 = read.csv("poi1.csv")

將poi導(dǎo)出，還沒(méi)有做處理，長(zhǎng)這樣

加了一列分組的Group，第一列手動(dòng)加了個(gè)sample做列名,v1即主坐標(biāo)1，v2主坐標(biāo)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")

權(quán)重值，從大到小即PCoA1,PCoA2......

畫(huà)出來(lái)就是這樣的

7.

摘自一篇文獻(xiàn)的物種注釋部分

圖中用的降維算法是NMDS，和PCoA差不多，只是原理不一樣

現(xiàn)在還沒(méi)有實(shí)現(xiàn)的是紅線(xiàn)標(biāo)注的圖，即比較組間差異和組內(nèi)差異，圖中的意思是組內(nèi)的差異小于組間的差異，更一步說(shuō)明了這兩組的物種組成不一樣

2021.12.17更新：https://mp.weixin.qq.com/s/sXlN4m7L72rvvKXGp2WAXA? 紅色標(biāo)注部分參考這個(gè)鏈接，照抄就能做出來(lái)，此處略

8.

物種多樣性的計(jì)算同理

OTU/ASV表同理