工作中有個真理：如果你連自己所做的工作的來龍去脈都講不清楚，那你是絕對不可能把這份工作做好的。
這適用於任何行業。如果你支支吾吾，講不清楚，那么說難聽點，你在混日子，沒有靜下心來工作。
檢驗標准：隨時向別人解釋你的工作，讓別人提出尖銳的問題，看你是不是答不上來。
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，這真是個好工具啊，還有詳細的教程。
還有詳細的數據庫： 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 百科
附錄：
稀釋性曲線(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()