筆記內(nèi)容：
由二代測序產(chǎn)生的序列數(shù)據(jù)（.fastq）到OTU table， 距離矩陣，物種多樣性指數(shù)，序列的進(jìn)化樹及物種注釋信息的可分析數(shù)據(jù)，為常規(guī)分析流程?？梢允褂胾search, vsearch, qiime等分析軟件實現(xiàn)，在必要的時候需要根據(jù)序列信息的具體情況編寫腳本予以實現(xiàn)。

本筆記描述的是微生物組16S rRNA的數(shù)據(jù)分析階段，大概包括哪些內(nèi)容，如何通過python和R來實現(xiàn)，大概會遇到什么問題以及如何解決。注意本分析是基于OTU table， 樣本的距離矩陣，物種多樣性指數(shù)，序列的進(jìn)化樹及物種注釋信息而非基于序列數(shù)據(jù)（.fastq等）。
包括以下兩個內(nèi)容：

input文件介紹

OTU(Operational taxonomic unit)，可操作分類單元
各OTU的代表序列（.fasta）
物種注釋信息
物種進(jìn)化樹
各樣本的多樣性指數(shù)：α-diversity的input
距離矩陣：β-diversity的input

分析內(nèi)容

物種構(gòu)成及優(yōu)勢物種
α-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標(biāo)記，以示它是來自哪個sample。經(jīng)過一些預(yù)先處理，包括去除barcode, 低質(zhì)量序列，污染序列，嵌合體，等。使用序列聚類算法將相似度（similarity）為97%以上的序列放在一起，組成一個OTU。所以一個OTU內(nèi)所有的序列均為相似度97%以上的，相似度不足97%的則分到其他的OTU中去。于是我們可以得到OTU_table，一個給出每個sample中每個OTU包含多少reads數(shù)目的矩陣：

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

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

各OTU的代表序列（.fasta）

以下為一個代表序列內(nèi)容示例：

物種注釋信息

將OTU代表序列分別與數(shù)據(jù)庫進(jìn)行比對，給每個OTU追溯到其物種來源。劃分到界（Kingdom）、門（Phylum）、綱（Class）、目（Order）、科（Family）、屬（Genus）、種（Species）。雖然不同軟件及流程output的物種注釋格式可能不一致，但內(nèi)容大同小異，均包含了上述信息。
你可以結(jié)合OTU table和物種注釋信息，將相同level的物種豐度相加，整理出每個level的物種豐度文件。比方說將family level中，相同family的物種豐度相加，形成一個family level的物種豐度文件。其作用為可以通過直接比較不同分組的物種豐度，從而找出哪些物種的豐度在組間存在差異，即挑選可以區(qū)分不同組的marker.

其形式一般為每個OTU對應(yīng)以下一條物種注釋信息。有些公司測序并初步分析給出OTUtable, 在每行OTU后面直接附上了注釋信息。
k__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__Robinsoniella; s__peoriensis

物種進(jìn)化樹

為了研究OTU序列所代表的物種進(jìn)化關(guān)系，我們通過OTU代表序列之間的相似性構(gòu)建物種進(jìn)化樹，代表每個OTU的進(jìn)化關(guān)系。其后綴名為.tree, 可以用figtree軟件打開。

各樣本的多樣性指數(shù)：α-diversity的input

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

距離矩陣: β-diversity的input

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

分析內(nèi)容

物種構(gòu)成及優(yōu)勢物種

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

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

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

α-diversity

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

rarefaction curve（稀釋曲線）
可參考α diversity分析https://www.jianshu.com/p/7cb452fede5a

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

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

β-diversity

目前常以距離矩陣作為input，使用PERMANOVA做組間差異比較，判斷組間物種構(gòu)成是否有差異，并以主坐標(biāo)軸分析Principal Coordinates Analysis (PCoA, = Multidimensional scaling, MDS) 產(chǎn)生的新坐標(biāo)繪圖進(jìn)行可視化。
PCoA(MDS)是根據(jù)樣本之間的距離，將高維數(shù)據(jù)投射到低維，達(dá)到降維的目的。在PcoA散點圖中，聚集在一起的樣本相似程度越高，在微生物β-diverisity分析中，即為聚集在一起的樣本物種構(gòu)成更相似。
用R中的ade4及ggplot2繪制PcoA散點圖及置信橢圓，具體代碼實現(xiàn)見R使用筆記: scatterplot with confidence ellipses;envfit的實現(xiàn)及釋意。輸出例子如下所示：

PERMANOVA(Permutational (nonparametric) MANOVA)是適用于距離矩陣的非參數(shù)統(tǒng)計檢驗方法, 可以通過R中vegan包的adonis()實現(xiàn)。具體統(tǒng)計原理不做詳敘，可以參考這個鏈接。其在R中使用要點有三個：

1. input文件必須為距離矩陣，不是的話需要使用dist()轉(zhuǎn)換。

2. MAN <- adonis(y.std ~ grp, permutations = 1000) 實現(xiàn)

3. 如果結(jié)果p <0.05, 則說明兩個組別的物種構(gòu)成存在差異，且根據(jù)R^2值可以知曉這種差異多少比例是由組別的因素貢獻(xiàn)的。

對16s微生物組數(shù)據(jù)而言，組間物種構(gòu)成的差異以PERMANOVA的統(tǒng)計檢驗結(jié)果為準(zhǔn)，PcoA(MDS)所作的二維或三維散點圖為可視化手段，為更直觀的展現(xiàn)組間差異。

biomarker

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

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

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

test <- read.csv('test_marker.csv', header = TRUE, row.names = 1, check.names = FALSE)
### heatmap
library(gplots)
library(RColorBrewer)
library(ggplot2)

tax = c('Tax_1','Tax_2','Tax_3')
test_data <- test[,tax]
grouping <- as.factor(test$group)

plot_color <- c(brewer.pal(5, "Set1")[4],brewer.pal(5, "Set1")[3])[grouping]  
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 groups
test_m <- read.csv('test_marker.csv', header = TRUE, check.names = FALSE)
colnames(test_m)[1] = 'ID'
test_m <- melt(test,id.vars = '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 <- t.test(tax_ra[[i]] ~ group, data = tax_ra)
   P.value = c(P.value, t$p.value)
   tax = c(tax, i)
 }
 
 p <- data.frame(tax, P.value)
 return(p)
}

test_for_marker(test)

當(dāng)然不一定非要用heatmap和boxplot來展示結(jié)果，可視化方法是靈活的，需要根據(jù)具體科學(xué)問題來變通。
傳統(tǒng)統(tǒng)計學(xué)方法可能會稍微有些麻煩，但是可以很細(xì)致的把每層物種都篩查一遍尋找marker；相比包裝好的軟件，自主調(diào)節(jié)一些分析和作圖的細(xì)節(jié)也更加靈活。

boruta和lefse
可以參考boruta：用于特征物種的挑選， 它大概是將每個feature(物種)的數(shù)值打亂隨機排列，然后看該feature(物種)的原始評分是否比打亂的隨機數(shù)值更加有規(guī)律（顯著）。它的好處在于不必拘泥于傳統(tǒng)統(tǒng)計學(xué)檢驗方法的先驗條件，對于變量特別多，容易過擬合的數(shù)據(jù)也很友好。

Lefse有在線使用的網(wǎng)站，也可以下載軟件到linux系統(tǒng)中在終端使用。它的原理是先使用非參數(shù)檢驗Kruskal-Wallis test檢測所有物種在指定組別中的顯著性，然后用wilcoxon rank-sum test兩兩比較。LEfSe采用線性判別分析（LDA）來估算每個物種的豐度對差異效果影響的大小，應(yīng)該是線性模型的原理。其網(wǎng)頁版基本上跟著以上鏈接上傳文件，按照提示指引輸入分組信息及各項參數(shù)，就可以得到biomarker分析結(jié)果。其優(yōu)點在于方便省心，圖好看且使用廣泛。

對與marker的挑選沒有一個固定的方法標(biāo)準(zhǔn)，比方說什么情況下適合用哪種方法之類的。就使用了16s數(shù)據(jù)的各種文獻(xiàn)來看，用哪種方法的都有。

功能分析：PICRUSt

PICRUST也是可以安裝在linux系統(tǒng)上的軟件，用于物種功能預(yù)測分析。在16s數(shù)據(jù)中并不十分推薦，因為其在16s中的功能預(yù)測有一定局限性。但是目前也有很多文獻(xiàn)將功能分析納入了分析結(jié)果中。功能分析僅供參考，需要結(jié)合科學(xué)問題慎重考慮。

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