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

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

檢驗標準：隨時向別人解釋你的工作，讓別人提出尖銳的問題，看你是不是答不上來。

---

16S概念

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

有點生物學基礎的會知道16S和核糖體有關，但大多數還是搞不清楚它們之間的關系。

先明確一些概念：

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

16S rRNA：這并不是我們的研究對象，因為我們測序的不是它，而是它對應在DNA雙鏈上的基因序列，

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

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

所以，16S測序的大致邏輯就是：

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

分析邏輯就是聚類成OTU，然后注釋（比對已知數據庫），后續分析。

---

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

1、16S

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

因16S rDNA分子大小適中，突變率小，故成為細菌系統發育和分類鑒定最常用的標簽。

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

16S rDNA序列結構

2、OTU

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

通過測序獲得的大量reads，如何才能轉變為我們需要的物種信息呢？首先需要對這些reads進行歸類（cluster），通常在97%的相似水平劃分為不同的OTU，將OTU代表序列與相應的微生物數據庫比對（Silva、RDP、Greengene等），得到每個樣本所含的物種信息，進而進行后續生物信息統計分析。

3、Q值

Q值評估用來測序的堿基質量，Q值與測序錯誤E值之間關系為如果一個堿基的Q值為20，那表示這個堿基的可能測錯的可能性為1%。實際操作中常用Q20/Q30作為標準，Q20大于90%妥妥的。

4、Coverage

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

5、Alpha-diversity

Alpha多樣性用于度量群落生態單樣本的物種多樣性，是反映豐富度和均勻度的綜合指標。

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

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

根據各樣本生成的OTU，對樣本序列進行隨機取樣，以取出的序列數及這些序列所能代表的OTU數構建曲線，計算樣本的Alpha多樣性。

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

6、Beta-diversity

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

7、LEfSe

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

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

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

---

16S分析流程

建庫測序暫且跳過，一下是我們拿到測序數據后的分析步驟：

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

Hiseq 的 16S分析流程已經比較成熟了，網上能搜到一些pipeline。

PacBio的 16S分析是最近才開發出來的。

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

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

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

Mothur manual

還有詳細的數據庫： Alignment database（greengenes 和 SILVA）

---

16S分析細節

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

以下一次詳解：

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

得到OTU后，將原來的序列比對回得到的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

轉化格式

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

待續~

參考：

16Sr DNA 百科

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

附錄：

稀釋性曲線(rarefaction curve)

稀釋性曲線(rarefaction curve)：一般是從樣本中隨機抽取一定數量的個體，統計出這些個體所代表物種數目，并以個體數與物種數來構建曲線。它可以用來比較測序數量不同的樣本物種的豐富度，也可以用來說明樣本的取樣大小是否合理。分析采用對優化序列進行隨機抽樣的方法，以抽到的序列數與它們所能代表OTU的數目構建rarefaction curve。稀釋性曲線圖中，當曲線趨向平坦時，說明取樣的數量合理，更多的取樣只會產生少量新的OTU，反之則表明繼續取樣還可能產生較多新的OTU。因此，通過作稀釋性曲線，可以得出樣品的取樣深度情況。 默認是在 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()

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。