matplot畫圖控制marker點的個數_微生物組16S rRNA資料分析小結：從OTU table到marker物種...

導讀：由二代測序産生的序列資料(.fastq)到OTU table， 距離矩陣，物種多樣性指數，序列的進化樹及物種注釋資訊的可分析資料，為正常分析流程。可以使用usearch, vsearch, qiime等分析軟體實作，在必要的時候需要根據序列資訊的具體情況編寫腳本予以實作。

本筆記描述的是微生物組16S rRNA的資料分析階段，大概包括哪些内容，如何通過python和R來實作，大概會遇到什麼問題以及如何解決。注意本分析是基于OTU table， 樣本的距離矩陣，物種多樣性指數，序列的進化樹及物種注釋資訊而非基于序列資料(.fastq等)。

包括以下兩個内容：

input檔案介紹

OTU(Operational taxonomic unit)，可操作分類單元

各OTU的代表序列(.fasta)

物種注釋資訊

物種進化樹

各樣本的多樣性指數：α-diversity的input

距離矩陣：β-diversity的input

分析内容

物種構成及優勢物種

α-diversity

β-diversity

biomarker

功能分析：PICRUSt

input檔案介紹

OTU(Operational taxonomic unit)，操作分類單元

在二代測序中，每個sample都會測到許多許多序列：

sample1: seq1, seq2, seq3, seq4, seq5...

sample2: seq1, seq2, seq3, seq4, seq5...

sample3: seq1, seq2, seq3, seq4, seq5...

...

每個序列都會有一小段barcode标記，以示它是來自哪個sample。經過一些預先處理，包括去除barcode, 低品質序列，污染序列，嵌合體，等。使用序列聚類算法将相似度(similarity)為97%以上的序列放在一起，組成一個OTU。是以一個OTU内所有的序列均為相似度97%以上的，相似度不足97%的則分到其他的OTU中去。于是我們可以得到OTU_table，一個給出每個sample中每個OTU包含多少reads數目的矩陣：

即每個sample對應每個OTU中的序列reads數目。如sample1在OTU1中有2個序列reads數目。如下所示的OTU_table即豐度。相對豐度則以每個sample(每行)為100%，計算各OTU的reads數目占一個sample中所有的reads數目的百分比。

OTU是對相似性序列進行聚類，将海量測序序列聚類成數量較少的分類單元，并且每個OTU提供一個代表序列，基于它進行後續物種注釋及分析，更加簡便和清晰。

各OTU的代表序列(.fasta)

以下為一個代表序列内容示例：

物種注釋資訊

将OTU代表序列分别與資料庫進行比對，給每個OTU追溯到其物種來源。劃分到界(Kingdom)、門(Phylum)、綱(Class)、目(Order)、科(Family)、屬(Genus)、種(Species)。雖然不同軟體及流程output的物種注釋格式可能不一緻，但内容大同小異，均包含了上述資訊。

你可以結合OTU table和物種注釋資訊，将相同level的物種豐度相加，整理出每個level的物種豐度檔案。比方說将family level中，相同family的物種豐度相加，形成一個family level的物種豐度檔案。其作用為可以通過直接比較不同分組的物種豐度，進而找出哪些物種的豐度在組間存在差異，即挑選可以區分不同組的marker.

其形式一般為每個OTU對應以下一條物種注釋資訊。有些公司測序并初步分析給出OTUtable, 在每行OTU後面直接附上了注釋資訊。

k__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__Robinsoniella; s__peoriensis

物種進化樹

為了研究OTU序列所代表的物種進化關系，我們通過OTU代表序列之間的相似性建構物種進化樹，代表每個OTU的進化關系。其字尾名為.tree, 可以用figtree軟體打開。

各樣本的多樣性指數：α-diversity的input

α-diversity是用于回答“一個樣本中有多少個物種？”的。是以α多樣性指數是針對每個樣本的。最簡單的一種指數為richness，即每個樣本中OTU的個數。

距離矩陣: β-diversity的input

β-diversity是用于回答“兩個樣本之間的相似程度如何？”這樣的問題。它比較兩個樣本之間的相似度或者差異程度，并給這兩個樣本計算一個值，通常在0-1之間，以距離矩陣的形式呈現。可以觀察到它是對稱的，因為兩兩樣本之間相似性的值是一樣的。

分析内容

物種構成及優勢物種

以下為一篇文章中的示例：

Phylum level microbial compositions of faeces, lavage and tissue samples. Tong M, Li X, Wegener Parfrey L, Roth B, Ippoliti A, et al. (2013) A Modular Organization of the Human Intestinal Mucosal Microbiota and Its Association with Inflammatory Bowel Disease. PLoS ONE 8(11): e80702. doi:10.1371/journal.pone.0080702

可以大概了解為根據以上5個組的phylum level相對豐度繪制barplot, 觀察到相對豐度最大的物種則為優勢物種。在上圖中為Firmicutes和Bacteroidetes。

α-diversity

包括rarefaction curve, rank abundance curve, 各項多樣性指數的組間差異等。

rarefaction curve(稀釋曲線)

可參考α diversity分析https://www.jianshu.com/p/7cb452fede5a

在每個樣本中不斷抽樣，每次都随機抽取一定數量的序列，以抽取到的序列建構OTU。其核心在于resampling。随着抽取的序列數目不斷增加，其建構的OTU個數從迅速增加到趨于平坦，則說明抽樣的數目合理，更多的序列不會再增加更多信OTU個數。即測序深度達到了要求。其橫軸為每次抽取的read counts, 縱軸為以抽取的read counts建構的OTU個數。qiime可以生成rarefaction curve, 也可以用R實作。

rank abundance curve

rank abundance curve比較簡單，其橫軸為按照相對豐度從大到小排序的OTU的ID(或者其他物種ID)，縱軸為相對豐度，如下所示。以下這個例子沒有标出X軸的OTU ID, 隻注明了其rank(排序)。

來自wiki

通過它可以了解優勢物種有哪些。如果rank abundance curve很陡峭(即一開始很高，然後一個大跳水降很低)，說明在樣本中有明确的優勢物種，且占了很大的比例。拖尾的物種豐度比較稀少。如果它下降的比較平緩，說明各物種都占有一定比例。當然你也可以根據樣本分組在一個圖中畫多個rank abundance curve. 還可以根據一組樣本，以相對豐度均值為縱軸畫圖，如下所示：

出自https://rpubs.com/marschmi/133626

各項多樣性指數的組間差異

根據各樣本的多樣性指數：α-diversity的input(即每個樣本對應一種多樣性指數的數值)和分組資訊檔案(即根據科學問題，将樣本分成不同的類别，而分析的意義在于尋找不同類别之間的差異)。

可以用R将結果可視化。有時候是分多個組，需要注意使用哪種統計方法。可以參考統計檢驗簡單小結

fake data

β-diversity

目前常以距離矩陣作為input，使用PERMANOVA做組間差異比較，判斷組間物種構成是否有差異，并以主坐标軸分析Principal Coordinates Analysis (PCoA, = Multidimensional scaling, MDS) 産生的新坐标繪圖進行可視化。

PCoA(MDS)是根據樣本之間的距離，将高維資料投射到低維，達到降維的目的。在PcoA散點圖中，聚集在一起的樣本相似程度越高，在微生物β-diverisity分析中，即為聚集在一起的樣本物種構成更相似。

用R中的ade4及ggplot2繪制PcoA散點圖及置信橢圓，具體代碼實作見R使用筆記: scatterplot with confidence ellipses;envfit的實作及釋意。輸出例子如下所示：

PERMANOVA(Permutational (nonparametric) MANOVA)是适用于距離矩陣的非參數統計檢驗方法, 可以通過R中vegan包的adonis()實作。具體統計原理不做詳叙，可以參考這個連結。其在R中使用要點有三個：

1. input檔案必須為距離矩陣，不是的話需要使用dist()轉換。

2. MAN 實作

3. 如果結果p <0.05, 則說明兩個組别的物種構成存在差異，且根據R^2值可以知曉這種差異多少比例是由組别的因素貢獻的。

對16s微生物組資料而言，組間物種構成的差異以PERMANOVA的統計檢驗結果為準，PcoA(MDS)所作的二維或三維散點圖為可視化手段，為更直覺的展現組間差異。

biomarker

marker, 在微生物組物種層面上，如果某個物種的相對豐度在兩組(或多組)間存在有統計學意義的差異，即該物種的豐度高低可以有效區分不同組别，則這個物種為這兩組(或多組)的biomarker.

目前我們常用傳統統計學方法，boruta及lefse來挑選marker物種。傳統的統計學及boruta方法靠python或者R代碼完成并可視化，lefse為包裝好的軟體，可以使用web端版本直接輸入input檔案，選擇參數運作即可得結果，也可以在linux系統中安裝lefse軟體，在終端運作并得到結果。

傳統統計學方法

傳統統計學方法則簡單的以物種豐度表為Input，按照不同分組，對每個物種的組間進行統計檢驗，如果組間差異有統計學意義，則該物種為biomarker。以下為一個input及結果展示示例，可以用heatmap及boxplot來了解不同組間物種豐度的差異方向，并使用統計檢驗(本例中使用t.test)了解差異是否具有統計學意義，并需要做多重統計檢驗矯正(p.adjust)。用R實作的代碼參考如下：

fake data example

test 'test_marker.csv', header = ### heatmaplibrary(gplots)library(RColorBrewer)library(ggplot2)tax = c('Tax_1','Tax_2','Tax_3')test_data grouping as.factor(test$group)plot_color 5, heatmap.2(as.matrix(test_data),         trace="none",         scale="column",         RowSideColors = plot_color,         dendrogram = "none",         symbreaks = TRUE,         col=rev(brewer.pal(10, "RdBu")), breaks = seq(-1.5,1.5,0.3),         density.info=c("none"),         margins=c(8,16), cexRow = 0.8, srtCol = 45, Colv=FALSE, Rowv = FALSE)### boxplot and test for comparint between groupstest_m 'test_marker.csv', header = colnames(test_m)[1] = 'ID'test_m 'group')ggplot(test_m, aes(x = variable, y = value, fill = group)) +      geom_boxplot(position = position_dodge(0.8)) +      geom_point(position = position_jitterdodge()) +      scale_fill_brewer(palette = "Set2")test_for_marker function(tax_ra){ P.value = c() tax = c() tt = colnames(tax_ra)[2:length(colnames(tax_ra))] for (i in tt){   t    P.value = c(P.value, t$p.value)   tax = c(tax, i) }  p  return(p)}test_for_marker(test)
           

當然不一定非要用heatmap和boxplot來展示結果，可視化方法是靈活的，需要根據具體科學問題來變通。

傳統統計學方法可能會稍微有些麻煩，但是可以很細緻的把每層物種都篩查一遍尋找marker；相比包裝好的軟體，自主調節一些分析和作圖的細節也更加靈活。

boruta和lefse

可以參考boruta：用于特征物種的挑選， 它大概是将每個feature(物種)的數值打亂随機排列，然後看該feature(物種)的原始評分是否比打亂的随機數值更加有規律(顯著)。它的好處在于不必拘泥于傳統統計學檢驗方法的先驗條件，對于變量特别多，容易過拟合的資料也很友好。

Lefse有線上使用的網站，也可以下載下傳軟體到linux系統中在終端使用。它的原理是先使用非參數檢驗Kruskal-Wallis test檢測所有物種在指定組别中的顯著性，然後用wilcoxon rank-sum test兩兩比較。LEfSe采用線性判别分析(LDA)來估算每個物種的豐度對差異效果影響的大小，應該是線性模型的原理。其網頁版基本上跟着以上連結上傳檔案，按照提示指引輸入分組資訊及各項參數，就可以得到biomarker分析結果。其優點在于友善省心，圖好看且使用廣泛。

對與marker的挑選沒有一個固定的方法标準，比方說什麼情況下适合用哪種方法之類的。就使用了16s資料的各種文獻來看，用哪種方法的都有。

功能分析：PICRUSt

PICRUST也是可以安裝在linux系統上的軟體，用于物種功能預測分析。在16s資料中并不十分推薦，因為其在16s中的功能預測有一定局限性。但是目前也有很多文獻将功能分析納入了分析結果中。功能分析僅供參考，需要結合科學問題慎重考慮。

另外也有人做出了分析網站http://www.microbiomeanalyst.ca/faces/home.xhtml可以直接上傳符合格式要求的OTU table, meta檔案及物種資訊就可以輸出一系列上文提到的各種分析結果，操作輕松，不用寫代碼。但是在資料形式不尋常，比方說對于時間序列，單個樣本的重複測量結果，分組複雜的資料，需要進行資料預處理的資料(比如需要取對數等)并不能完全依賴包辦型的分析網站。不論是代碼、軟體還是包辦的分析網站都隻是分析工具，需要根據具體科學問題，合理使用它們。

Ref:

R作圖，PERMANOVA的adonis()實作，rank abundance等參考

PcoA(mds)的R實作

生态學資料的多元分析：包括PERMANOVA的學習網站，含原理。

推薦閱讀

1.ggplot2繪制曼哈頓圖示例 2.phyloseq | 用 R 分析微生物組資料及可視化 3.R語言PCA分析教程 | Principal Component Methods in R 4.vegan包進行微生物群落主坐标分析(PCoA)及ggplot2作圖 5.R語言 | 相關性分析 6.排序分析(Ordination Analysis) 7.微生物組16S rRNA資料分析小結 8.微生物組16S rRNA資料分析小結：從fastq測序資料到OTU table 9.dplyr包的copy版--poorman包應用執行個體

連結：

https://www.jianshu.com/p/920a5ce3a7a0

著作權歸作者所有。

我就知道你在看?