微生物組學研究-16s測序 從建庫到數(shù)據(jù)分析全流程解析

1）首先想了解在不同組樣本中各有哪些微生物存在和豐富度（對應于菌群鑒定和α多樣性分析）；

2）接著想看不同樣本組間微生物群組成是否存在差異（對應于β多樣性分析）；

3）如果是，那么就有必要找出引起不同組樣本微生物群差異的關鍵菌。如果不是，那說明微生物群比如腸道菌群與疾病或表型可能并不相關（基于已有的研究，這種可能性比較?。?；

4）找到了關鍵菌，在臨床上，很自然會想到，這些（個）關鍵菌是否可以作為Biomarker（對應于疾病診斷模型構建），比如用于區(qū)分糖尿病前期患者與健康組的標志物；

5）以及這些（個）菌是否與臨床指標具有相關性（對應于菌群與臨床指標的相關性分析）；也會進一步想到，既然不同組的微生物群落存在差異，又與疾病具有相關性，

6）那么這些菌群是如何影響宿主的，可能參與了哪些代謝途徑（對應于菌群基因功能預測）；

7）這些預測到的菌群功能是否與疾病有關，通常是肯定的。最后把這些結果整合起來分析，可以初步得出菌群組成的變化是如何與疾病或表型相關的。

順著上述7個生物學問題來看16S測序結果，你會輕松撥開迷霧，直達核心結果。

ITS: 鑒定種以下的水平

18S：鑒定種和種以上的水平

16S：鑒定屬水平，少量可以達到種水平

• 測序物種的差異：16S測序可以對目標區(qū)域（可變區(qū)域）進行擴增，可以鑒定低豐度物種種類，結果高于宏基因組的結果。
• 功能分析的差異：16S 只能預測代謝通路，得不到物種對代謝通路的信息；宏基因組是全基因組的測序，可以得到具體的代謝通路的情況

• 特定時空下，并不是所有的基因都表達
• 宏轉錄組測序：特定時空下， 有活性的微生物群落的組成

• 16S 測序測的是編碼16S rRNA的DNA序列

• 5S: 堿基少 ，不足以反映菌群的差異性
• 23S：堿基數(shù)多，變異多，不容易區(qū)分親緣關系
• 16S：對保守區(qū)域設計探針，然后擴增可變區(qū)域；堿基數(shù)量少，變異少

• 測序結果的reads直接進行97%相似程度歸為劃分小組，一個小組為一個OTU
• 一個OUT可以注釋到一個物種，一個物種可以被多個OUT注釋

#### OUT解釋

OTU(operational taxonomic units) 是在系統(tǒng)發(fā)生學研究或群體遺傳學研究中，為了便于進行分析，人為給某一個分類單元（品系，種，屬，分組等）設置的同一標志。通常按照 97% 的相似性閾值將序列劃分為不同的 OTU，每一個 OTU 通常被視為一個微生物物種。相似性小于97%就可以認為屬于不同的種，相似性小于93%-95%，可以認為屬于不同的屬。樣品中的微生物多樣性和不同微生物的豐度都是基于對OTU的分析。

下表是對每個樣本在分類字水平上的數(shù)量進行統(tǒng)計，并且在表栺中列出了在每個分類字水平上的物種數(shù)目

• 其中SampleName表示樣本名稱；Phylum表示分類到門的OTU數(shù)量；Class表示分類到綱的OTU數(shù)量；Order表示分類到目的OTU數(shù)量；Family表示分類到科的OTU數(shù)量；Genus表示分類到屬的OTU數(shù)量；Species表示分類到種的OTU數(shù)量。

• α多樣性關注【樣本自身】的菌群豐富度和均勻度，而β多樣性關注【樣本間】的菌群組成與分布的差異。只有當樣本（組）間菌群組成存在差異，才有可能進一步探討菌群失調與疾病的關系。
• α多樣性是度量【單個樣本內】有多少種微生物物種，以及每個物種所占比例的指標。
• β多樣性是度量【不同樣本間】菌群組成的相似度大小的指標，即關注各樣本間的菌群組成差異
• a多樣性：樣本內的比較（單個樣本的比較）
• β多樣性：樣本間的比較（樣本間的比較）

Alpha多樣性（樣本內多樣性）

Alpha多樣性是指一個特定區(qū)域或者生態(tài)系統(tǒng)內的多樣性，常用的度量指標有Chao1 豐富度估計量（Chao1 richness estimator） 、香農 - 威納多樣性指數(shù)（Shannon-wiener diversity index）、辛普森多樣性指數(shù)（Simpson diversity index）等。

計算菌群豐度：Chao、ace；??

計算菌群多樣性：Shannon、Simpson。

Simpson指數(shù)值越大，說明群落多樣性越高；Shannon指數(shù)越大，說明群落多樣性越高。

Chao1：是用chao1 算法計算群落中只檢測到1次和2次的OTU數(shù)估計群落中實際存在的物種數(shù)。Chao1 在生態(tài)學中常用來估計物種總數(shù)，由Chao (1984) 最早提出。Chao1值越大代表物種總數(shù)越多。

Schao1=Sobs+n1(n1-1)/2(n2+1)

其中Schao1為估計的OTU數(shù)，Sobs為觀測到的OTU數(shù)，n1為只有一條序列的OTU數(shù)目，n2為只有兩條序列的OTU數(shù)目。

Shannon：用來估算樣品中微生物的多樣性指數(shù)之一。它與 Simpson 多樣性指數(shù)均為常用的反映 alpha 多樣性的指數(shù)。Shannon值越大，說明群落多樣性越高。

Ace：用來估計群落中含有OTU 數(shù)目的指數(shù)，由Chao 提出，是生態(tài)學中估計物種總數(shù)的常用指數(shù)之一，與Chao1 的算法不同。

Simpson：用來估算樣品中微生物的多樣性指數(shù)之一，由Edward Hugh Simpson ( 1949) 提出，在生態(tài)學中常用來定量的描述一個區(qū)域的生物多樣性。Simpson 指數(shù)值越大，說明群落多樣性越高。

Beta多樣性分析（樣品間差異分析）

也許我們有聽說Beta多樣性在最近10年間成為生物多樣性研究的熱點問題之一。具體解釋下：

Beta多樣性度量時空尺度上物種組成的變化, 是生物多樣性的重要組成部分, 與許多生態(tài)學和進化生物學問題密切相關！

PCoA分析

PCoA（principal co-ordinates analysis）是一種研究數(shù)據(jù)相似性或差異性的可視化方法，通過一系列的特征值和特征向量進行排序后，選擇主要排在前幾位的特征值，PCoA 可以找到距離矩陣中最主要的坐標，結果是數(shù)據(jù)矩陣的一個旋轉，它沒有改變樣品點之間的相互位置關系，只是改變了坐標系統(tǒng)。

重要的是，它是可以用來觀察個體或群體間的差異的。

PCA分析

主成分分析（Principal component analysis）PCA 是一種研究數(shù)據(jù)相似性或差異性的可視化方法，通過一系列的特征值和特征向量進行排序后，選擇主要的前幾位特征值，采取降維的思想，PCA 可以找到距離矩陣中最主要的坐標，結果是數(shù)據(jù)矩陣的一個旋轉，它沒有改變樣品點之間的相互位置關系，只是改變了坐標系統(tǒng)。

顯著差異菌群分析

通過β多樣性分析，可以確定不同組間的微生物群落是存在差異的，接著需要進一步找出哪些菌（群）引起了組間的群落差異。只有找出核心菌（群），才能明確下一步的研究方向。在報告中，使用目前在文獻中高頻出現(xiàn)的方法——LEfSe（Linear discriminant analysis Effect Size），來做菌群差異分析，尋找生物標志物（Biomarker）。該方法綜合了統(tǒng)計學上的差異分析和該差異物種對分組結果的影響力得分值，同時強調了統(tǒng)計學意義和生物相關性。

LEFse原理

Step1.?首先在多組樣本中采用?非參數(shù)因子Kruskal-Wallis秩和檢驗?檢測不同分組間豐度差異顯著的物種；也就是圖中按class1 和class2兩個大的分組，每一行都進行檢驗，初步得到差異物種，通過檢驗的打鉤進入step2檢驗；

Step2.?再利用Wilcoxon秩和檢驗，對每一組中的亞組進行兩兩檢驗，具有顯著差異的再進行下一輪檢驗。

Step3.?最后用線性判別分析（LDA）對數(shù)據(jù)進行降維并評估差異顯著的物種的影響力（即LDA score）。

前兩步的Kruskal-Wallis秩和檢驗、Wilcoxon秩和檢驗?比較簡單，類似T檢驗或者方差檢驗等，只不過T檢驗和方差分析為參數(shù)檢驗（要求數(shù)據(jù)符合方差齊性、正態(tài)分布），而在微生物多樣性分析中，樣品物種豐度分布不確定，多采用非參數(shù)檢驗，所以采用非參數(shù)的Kruskal-Wallis秩和檢驗、Wilcoxon秩和檢驗。比較復雜一點的就是最后的LDA分析。

LDA是一種監(jiān)督學習的降維技術，也就是說其數(shù)據(jù)集中的每個樣本是有類別輸出的。是在目前機器學習、數(shù)據(jù)挖掘領域經典且熱門的一個算法這點和PCA不同。PCA是不考慮樣本類別輸出的無監(jiān)督降維技術。LDA是有監(jiān)督的，所以LDA算法可以很好的利用樣本的分組信息，得到的結果更可靠，這就是LDA分析優(yōu)勢。理解了LDA分析的原理，就不難理解LEfSe的分析結果了。

LDA的全稱是Linear Discriminant Analysis （線性判別分析），是一種supervised learning （有監(jiān)督學習）。有些資料上也稱為是Fisher' s Linear Discriminant，由Ronald Fisher發(fā)明自1936年，是在目前機器學習、數(shù)據(jù)挖掘領域經典且熱的一個算法。

LDA的思想可以用一句話概括，就是”投影后類內方差最小，類間方差最大”。簡單來說就是一種投影， 是將一個高維的點投影到一個低維空間，我們希望映射之后，不同類別之間的距離越遠越好，同類別之中的距離越近越好。

是不是很抽象哇，來舉個栗子吧。假設我們有兩類數(shù)據(jù)：分別為紅色和藍色，如下圖所際，這些數(shù)據(jù)特征是二維的， 我們希望將這些數(shù)據(jù)投影到一維的一條直線，讓每一種類別數(shù)據(jù)的投影點盡可能的接近，而紅色和藍色數(shù)據(jù)中心之間的距離盡可能的大。

從直觀上可以看出，右圖要比左圖的投影效果好，因為右圖的紅色數(shù)據(jù)和色數(shù)據(jù)各個較為集中，類別之間的距離明顯。左圖則在邊界處數(shù)據(jù)混雜。當然在實際應用中，我們的數(shù)據(jù)多個類別的，我們的原始數(shù)據(jù)一般也是超過二維的，投影后的也一般不是直線，而是一-個低維的超平面。

菌群標志物預測能力評估

受試者工作特征（ROC）曲線分析是一種常用的統(tǒng)計學分析方法，在醫(yī)學研究中主要用于評價診斷試驗的效能。在報告中，通過繪制ROC曲線，并計算ROC曲線下面積（AUC），來確定哪種菌（群）具有最佳的診斷價值

上圖以靈敏度為縱坐標，特異度為橫坐標繪制曲線。ROC曲線越靠近左上角，試驗的準確性就越高。若AUC值為1.0，反映出對兩個群組的完美區(qū)分，且不存在預測誤差。若AUC值在1.0和0.5之間，在AUC>0.5的情況下，AUC越接近于1，說明診斷效果越好。AUC在0.5~0.7時有較低準確性，AUC在0.7~0.9時有一定準確性，AUC在0.9以上時有較高準確性。AUC=0.5時，說明診斷方法完全不起作用，無診斷價值。AUC<0.5不符合真實情況，在實際中極少出現(xiàn)。

菌群基因功能預測

因為菌群功能預測軟件PICRUSt（Phylogenetic Investigation of Communities by Reconstruction of Unobserved States）的出現(xiàn)，研究者能進一步基于16S測序數(shù)據(jù)預測菌群可能參與的代謝通路（盡管并沒有測定菌群基因信息），以便能初步討論菌群組成變化與疾病或表型是如何關聯(lián)在一起的。在聯(lián)川報告中，使用最新的PICRUSt 2，相比上一版，用于預測的參考基因組數(shù)據(jù)庫已擴展超過10倍，可以獲得包括COG，EC，KO，PFAM，TIGRFAM等數(shù)據(jù)庫對菌群的基因功能注釋結果。然后，再使用STAMP軟件進行差異分析，得到在不同樣本組中顯著差異的菌群基因功能。如果要系統(tǒng)研究菌群攜帶的基因及其功能，則應該做宏基因組測序。菌群基因功能預測

• 菌群基因功能預測的原理：根據(jù)16s測序的 MARK基因和已知的參考基因數(shù)據(jù)庫進行比對，來預測宏基因組可能得功能

• 聯(lián)川生物醫(yī)學16S測序報告內容

參考：

1.https://mp.weixin.qq.com/s/MQkb1pyVV2YRTIGeE887kQ

2.https://mp.weixin.qq.com/s/SakHS9QuqFIpQTN6gCVbpg

3.https://mp.weixin.qq.com/s/ryvA1DSE4r5MqUPiSb9qKw

4.https://www.omicsclass.com/article/112

4.https://www.bilibili.com/video/BV1Mi4y1M7vY/?spm_id_from=333.337.search-card.all.click&;vd_source=738c849c9def6c13e683914ca83d858e