工作中有個(gè)真理：如果你連自己所做的工作的來龍去脈都講不清楚，那你是絕對(duì)不可能把這份工作做好的。

這適用于任何行業(yè)。如果你支支吾吾，講不清楚，那么說難聽點(diǎn)，你在混日子，沒有靜下心來工作。

檢驗(yàn)標(biāo)準(zhǔn)：隨時(shí)向別人解釋你的工作，讓別人提出尖銳的問題，看你是不是答不上來。

---

16S概念

什么是16S？S是什么意思？ 16S分析是用來干嘛的？能分析什么？ 16S大致的分析原理是什么？

有點(diǎn)生物學(xué)基礎(chǔ)的會(huì)知道16S和核糖體有關(guān)，但大多數(shù)還是搞不清楚它們之間的關(guān)系。

先明確一些概念：

核糖體：Ribosome，由 RNA(rRNA）和蛋白質(zhì) 組成，配合 tRNA 來翻譯 mRNA。核糖體按沉降系數(shù)來分類，S就是沉降系數(shù)，原核70S，真核80S。我們一般研究微生物，70S，由50S和30S兩個(gè)亞基組成。再細(xì)分為 5S、16S、23S，我們的 16S 就是指核糖體的亞基的一個(gè)組分，16S rRNA。（記住，這是原核生物核糖體的一個(gè)組分）

16S rRNA：這并不是我們的研究對(duì)象，因?yàn)槲覀儨y(cè)序的不是它，而是它對(duì)應(yīng)在DNA雙鏈上的基因序列，

16S rDNA。可以這樣理解，我們所說的16S 就是指 16S rDNA。

分子鐘：即氨基酸在單位時(shí)間以同樣的速度進(jìn)行置換。16S 的進(jìn)化具有良好的時(shí)鐘性質(zhì)，在結(jié)構(gòu)與功能上具有高度的保守性，在大多數(shù)原核生物中rDNA都具有多個(gè)拷貝，5S、16S、23S rDNA的拷貝數(shù)相同。16S rDNA由于大小適中，約1.5Kb左右，既能體現(xiàn)不同菌屬之間的差異，又能利用測(cè)序技術(shù)較容易地得到其序列，故被細(xì)菌學(xué)家和分類學(xué)家接受。（來源百度）

所以，16S測(cè)序的大致邏輯就是：

拿到一個(gè)樣品，我們捕獲其16S區(qū)域（引物PCR），然后測(cè)序，16S既然有極好的保守性，那就可以用于鑒別不同的物種（相當(dāng)于一個(gè)物種的獨(dú)一無二的條形碼）（有很大一部分是鑒定不到物種的）。

分析邏輯就是聚類成OTU，然后注釋（比對(duì)已知數(shù)據(jù)庫），后續(xù)分析。

---

偶然看到一篇好的科普文，轉(zhuǎn)載一下：來自伯豪生物

1、16S

通常所說的16S是指16S rDNA（或16S rRNA），16S rRNA 基因是編碼原核生物核糖體小亞基的基因，長度約1542bp，包括9個(gè)可變區(qū)和10個(gè)保守區(qū)，保守區(qū)序列反映了物種間的親緣關(guān)系，而可變區(qū)序列則能反映物種間的差異。

因16S rDNA分子大小適中，突變率小，故成為細(xì)菌系統(tǒng)發(fā)育和分類鑒定最常用的標(biāo)簽。

16S測(cè)序是指選擇16S rDNA某個(gè)或某幾個(gè)變異區(qū)域，選擇通用引物對(duì)環(huán)境樣本（腸道、土壤、水體等）微生物進(jìn)行PCR擴(kuò)增，然后對(duì)PCR產(chǎn)物進(jìn)行高通量測(cè)序，并將得到的測(cè)序數(shù)據(jù)與已有的16S rDNA數(shù)據(jù)庫進(jìn)行比對(duì)分析，從而對(duì)環(huán)境群落多樣性進(jìn)行研究，核心是物種分析，包括微生物的種類，不同種類間的相對(duì)豐度，不同分組間的物種差異以及系統(tǒng)進(jìn)化等。

16S rDNA序列結(jié)構(gòu)

2、OTU

OTU即Operational Taxonomic Units的縮寫（千萬表手滑寫成OUT，否則就OUT了），在系統(tǒng)發(fā)生學(xué)或群體遺傳學(xué)研究中，為了便于進(jìn)行分析，人為給某一個(gè)分類單元（品系，屬，種、分組等）設(shè)置的同一標(biāo)志。理論上一個(gè)OTU代表一個(gè)微生物物種。

通過測(cè)序獲得的大量reads，如何才能轉(zhuǎn)變?yōu)槲覀冃枰奈锓N信息呢？首先需要對(duì)這些reads進(jìn)行歸類（cluster），通常在97%的相似水平劃分為不同的OTU，將OTU代表序列與相應(yīng)的微生物數(shù)據(jù)庫比對(duì)（Silva、RDP、Greengene等），得到每個(gè)樣本所含的物種信息，進(jìn)而進(jìn)行后續(xù)生物信息統(tǒng)計(jì)分析。

3、Q值

Q值評(píng)估用來測(cè)序的堿基質(zhì)量，Q值與測(cè)序錯(cuò)誤E值之間關(guān)系為如果一個(gè)堿基的Q值為20，那表示這個(gè)堿基的可能測(cè)錯(cuò)的可能性為1%。實(shí)際操作中常用Q20/Q30作為標(biāo)準(zhǔn)，Q20大于90%妥妥的。

4、Coverage

Coverage值是指各樣本文庫的覆蓋率，數(shù)值越高，則樣本中序列被檢測(cè)出來的概率越高，該數(shù)值可反映本次測(cè)序結(jié)果是否代表樣本的真實(shí)情況。數(shù)值越接近于100%，代表本次測(cè)序結(jié)果越符合樣本中微生物的實(shí)際情況。

5、Alpha-diversity

Alpha多樣性用于度量群落生態(tài)單樣本的物種多樣性，是反映豐富度和均勻度的綜合指標(biāo)。

菌群豐富度（Community richness）指數(shù)有：Chao、Ace，Chao或Ace指數(shù)越大，說明菌群豐富度越高。

菌群多樣性（Community diversity）指數(shù)有：Shannon、Simpson，Shannon值越大，說明群落多樣性越高；Simpson指數(shù)值越大，說明菌群多樣性越低。

根據(jù)各樣本生成的OTU，對(duì)樣本序列進(jìn)行隨機(jī)取樣，以取出的序列數(shù)及這些序列所能代表的OTU數(shù)構(gòu)建曲線，計(jì)算樣本的Alpha多樣性。

通常Alpha多樣性指數(shù)均以表格形式呈現(xiàn)，這顯然不利于論文顏值的提升，那么可以將Alpha多樣性指數(shù)結(jié)合起來，出個(gè)箱線圖（Box-plot）或可有所幫助。

6、Beta-diversity

Beta多樣性用于不同生態(tài)系統(tǒng)之間多樣性的比較，利用各樣本序列間的進(jìn)化關(guān)系及豐度信息來計(jì)算樣本間距離，反映樣本（組）間是否具有顯著的微生物群落差異。目前應(yīng)用比較多的是PCA、PCoA、NMDS分析等。由于微生物多樣性研究通常會(huì)涉及到大樣本數(shù)量的樣本，因此通過Beta-diversity分析可以直觀地反映樣本組間的差異情況。距離越遠(yuǎn)，微生物群落差異越大，即相似性越高。

7、LEfSe

LEfSe分析即LDA Effect Size分析，多用于多個(gè)分組（≧2）之間的比較，或者進(jìn)行亞組比較分析，進(jìn)而找到組間在豐度上有顯著差異的物種（biomaker）。

基本過程是首先在多組樣本中采用非參數(shù)因子Kruskal-Wallis秩和檢驗(yàn)檢測(cè)不同分組間豐度差異顯著的物種；然后基于獲得的顯著差異物種，利用成組的Wilcoxon秩和檢驗(yàn)進(jìn)行組間差異分析；最后采用線性判別分析（LDA）對(duì)數(shù)據(jù)進(jìn)行降維并評(píng)估差異顯著的物種的影響力大小（即LDA score）。

LEfSe分析包含了不同分類學(xué)水平（門到屬或者種）的物種信息，并且是顯著差異的物種，因而是一項(xiàng)信息量大、采用率（biger）高的分析，值得擁有~~~

---

16S分析流程

建庫測(cè)序暫且跳過，一下是我們拿到測(cè)序數(shù)據(jù)后的分析步驟：

過濾Filter 連接Connect Tags 聚類OTU Cluster 分類Taxonomy 多樣性Alpha_Beta 。。。

Hiseq 的 16S分析流程已經(jīng)比較成熟了，網(wǎng)上能搜到一些pipeline。

PacBio的 16S分析是最近才開發(fā)出來的。

參考文章：High-resolution phylogenetic microbial community profiling

以下是該文章的分析流程：

可以看到，PacBio 16S中大部分都在使用 mothur，這真是個(gè)好工具啊，還有詳細(xì)的教程。

Mothur manual

還有詳細(xì)的數(shù)據(jù)庫： Alignment database（greengenes 和 SILVA）

---

16S分析細(xì)節(jié)

step01.OTU_Cluster.sh step02.OTU_Taxonomy.sh step03.Alpha_Diversity.sh step04.OTU_Analysis.sh

以下一次詳解：

step01.OTU_Cluster.sh 的核心就是聚類得到 OTU，使用的是mothur，步驟很多。

得到OTU后，將原來的序列比對(duì)回得到的OTU上。

usearch7.0.1090/usearch7.0.1090_i86linux32 -usearch_global ./1.Clean_Tags/All_Sample.fasta -db ./2.OTU_Cluster/OTU.fasta -strand both -id 0.97 -uc ./2.OTU_Cluster/global.map.uc

usearch_global command

轉(zhuǎn)化格式

bin/usearch7.0.1090/uc2otutab.py ./2.OTU_Cluster/global.map.uc > ./2.OTU_Cluster/OTU_shared_all.xls

待續(xù)~

參考：

16Sr DNA 百科

16s rRNA分析流程和工具的介紹

附錄：

稀釋性曲線(rarefaction curve)

稀釋性曲線(rarefaction curve)：一般是從樣本中隨機(jī)抽取一定數(shù)量的個(gè)體，統(tǒng)計(jì)出這些個(gè)體所代表物種數(shù)目，并以個(gè)體數(shù)與物種數(shù)來構(gòu)建曲線。它可以用來比較測(cè)序數(shù)量不同的樣本物種的豐富度，也可以用來說明樣本的取樣大小是否合理。分析采用對(duì)優(yōu)化序列進(jìn)行隨機(jī)抽樣的方法，以抽到的序列數(shù)與它們所能代表OTU的數(shù)目構(gòu)建rarefaction curve。稀釋性曲線圖中，當(dāng)曲線趨向平坦時(shí)，說明取樣的數(shù)量合理，更多的取樣只會(huì)產(chǎn)生少量新的OTU，反之則表明繼續(xù)取樣還可能產(chǎn)生較多新的OTU。因此，通過作稀釋性曲線，可以得出樣品的取樣深度情況。 默認(rèn)是在 97%相似性水平下劃分OUT并制作各樣品的稀疏曲線。

畫OTU Rank curve

data <- read.table("Sample.OTU_rank_percent.xls",header=T,check.names=F)
data <- t(data)
data <- data[rowSums(is.na(data)) != 16,]
library("RColorBrewer")
col=c("#00FF40FF","#FF0000FF","#00FFFFFF")
x_point=as.numeric(rownames(data))
x=max(x_point,na.rm=T)
x=x*1.1
y=max(data[,1:ncol(data)],na.rm=T)

pdf("Sample.OTU_rank.pdf",width=6+1*1,height=6)
par(mai=c(1,1,0.5,1*1))
par("usr")
plot(x,y,xlim=c(0,x),ylim=c(0.0005,100),type="n",xlab="Number of OTUs",ylab="Relative abundance(%)",main="OTU Rank Curve",tcl=-0.2,mgp=c(3,0.3,0),cex.lab=1.2,las=1,log="y",yaxt="n")
axis(2,tcl=-0.2,at=c(0.001,0.01,0.1,1,10,100),labels=c(0.001,0.01,0.1,1,10,100),las=1,hadj=0.6)

points(x_point,data[,1],lwd=2,lty=1,cex=0.6,col=col[1],type="l") #C1
points(x_point,data[,2],lwd=2,lty=1,cex=0.6,col=col[2],type="l") #C3.1
points(x_point,data[,3],lwd=2,lty=1,cex=0.6,col=col[3],type="l") #C3.2

legend(par("usr")[2],10^par("usr")[4],legend=colnames(data),col=col[1:3],pch=15,pt.cex=1.2,cex=0.9,bty="n",ncol=1,xpd=TRUE)
dev.off()

png("Sample.OTU_rank.png",type="cairo",width=6+1*1,height=6,units="in",res=120)
par(mai=c(1,1,0.5,1*1))
plot(x,y,xlim=c(0,x),ylim=c(0.0005,100),type="n",xlab="Number of OTUs",ylab="Relative abundance(%)",main="OTU Rank Curve",tcl=-0.2,mgp=c(3,0.3,0),cex.lab=1.2,las=1,log="y",yaxt="n")
axis(2,tcl=-0.2,at=c(0.001,0.01,0.1,1,10,100),labels=c(0.001,0.01,0.1,1,10,100),las=1,hadj=0.6)

points(x_point,data[,1],lwd=2,lty=1,cex=0.6,col=col[1],type="l") #C1
points(x_point,data[,2],lwd=2,lty=1,cex=0.6,col=col[2],type="l") #C3.1
points(x_point,data[,3],lwd=2,lty=1,cex=0.6,col=col[3],type="l") #C3.2

legend(par("usr")[2],10^par("usr")[4],legend=colnames(data),col=col[1:16],pch=15,pt.cex=1.2,cex=0.9,bty="n",ncol=1,xpd=TRUE)
dev.off()

總結(jié)

以上是生活随笔為你收集整理的16S 基础知识、分析工具和分析流程详解的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。