簡(jiǎn)介

16S rRNA位于原核細(xì)胞核糖體小亞基上，包括 10 個(gè)保守區(qū)域(Conserved Regions)和 9 個(gè)高變區(qū)域(Hypervariable Regions)，其中保守區(qū)在細(xì)菌間差異不大，高變區(qū)具有屬或種的特異性，隨親緣關(guān)系不同而有一定的差異。因此，16S rDNA被認(rèn)為是最適于細(xì)菌系統(tǒng)發(fā)育和分類鑒定的指標(biāo)，用于揭示生物物種的特征核酸序列[1]，常見的是基于第4可變區(qū)(V4)約290bp或者第3和第4(V3+V4)可變區(qū)約460bp的測(cè)序策略[2]。

根據(jù)所擴(kuò)增的區(qū)域特點(diǎn)構(gòu)建小片段文庫，進(jìn)行Illumina雙末端測(cè)序(Paired_End)。經(jīng)過Reads質(zhì)控、拼接、過濾和OTU/ASV聚類，可以進(jìn)行物種注釋及豐度分析；最后，通過α多樣性(Alpha Diversity)和β多樣性分析(Beta Diversity)，揭示樣本中物種組成和樣本間群落結(jié)構(gòu)的差異，并進(jìn)行個(gè)性化分析和深度的數(shù)據(jù)挖掘。

測(cè)序流程

從DNA樣本到最終數(shù)據(jù)獲得的過程中，樣本采集和檢測(cè)、核酸提取、PCR、純化、建庫、測(cè)序每一個(gè)環(huán)節(jié)都會(huì)對(duì)數(shù)據(jù)質(zhì)量和數(shù)量產(chǎn)生影響，而數(shù)據(jù)質(zhì)量又會(huì)直接影響后續(xù)信息分析的結(jié)果。

分析流程

測(cè)序得到的原始數(shù)據(jù)(Raw Data)，存在一定比例的干擾數(shù)據(jù)(Dirty Data)，為了使信息分析的結(jié)果更加準(zhǔn)確、可靠，需要對(duì)原始數(shù)據(jù)進(jìn)行處理、過濾和質(zhì)控來獲得有效數(shù)據(jù)(Clean Data)。

然后基于有效數(shù)據(jù)進(jìn)行聚類和物種注釋分類分析,根據(jù)OTU/ASV聚類結(jié)果，一方面對(duì)每個(gè)OTU/ASV的代表序列做物種注釋，得到對(duì)應(yīng)的物種信息和基于物種的豐度分布情況。同時(shí)，對(duì)OTU/ASV進(jìn)行豐度、Alpha多樣性計(jì)算、共有特有OTU/ASV統(tǒng)計(jì)，以得到樣本內(nèi)物種豐富度和均勻度、不同樣本或分組間的共有和特有OTU/ASV等信息。另一方面，可以對(duì)等進(jìn)行多序列比對(duì)并構(gòu)建系統(tǒng)發(fā)生樹，通過PCoA、PCA、NMDS等降維分析和樣本聚類樹展示，可以探究不同樣本或組別間群落結(jié)構(gòu)的差異。為進(jìn)一步挖掘分組樣本間的群落結(jié)構(gòu)差異，選用T-test、Simper、MetaStat、LEfSe、Anosim和MRPP等統(tǒng)計(jì)分析方法對(duì)分組樣本的物種組成和群落結(jié)構(gòu)進(jìn)行差異顯著性檢驗(yàn)。

同時(shí)，也可結(jié)合環(huán)境因素進(jìn)行CCA/RDA/dbRDA分析和多樣性指數(shù)與環(huán)境因子的相關(guān)性分析，得到顯著影響組間群落變化的環(huán)境影響因子。擴(kuò)增子的注釋結(jié)果還可以和相應(yīng)的功能數(shù)據(jù)庫相關(guān)聯(lián)，可以選用PICRUST、Tax4Fun、FAPROTAX 、BugBase軟件對(duì)生態(tài)樣本中的微生物群落進(jìn)行功能預(yù)測(cè)分析。

數(shù)據(jù)預(yù)處理

示例數(shù)據(jù)來自MicrobiomeStatPlot的github倉(cāng)庫。

抽平和標(biāo)準(zhǔn)化

由于測(cè)序過程存在一定的不均衡性，為了保證樣本結(jié)果的一致性，便于后續(xù)分析比較和解釋，往往需要對(duì)OTU/ASV表進(jìn)行隨機(jī)重抽樣,即抽平處理；但當(dāng)樣品測(cè)序量相差比較大時(shí)候，容易造成數(shù)據(jù)的極大浪費(fèi)，此時(shí)可以利用Deseq2和edgeR等基于分布的標(biāo)準(zhǔn)化方法；因此，關(guān)于差異OTU/ASV分析前的抽平和標(biāo)準(zhǔn)化一般需要結(jié)合原始稀釋曲線和實(shí)驗(yàn)設(shè)計(jì)進(jìn)行選擇。

library(vegan)
otu_table <- read.delim("16S-amplicon-analysis/otutab.txt",  header=T, sep="\t",  row.names=1,  stringsAsFactors = FALSE)
head(colSums(otu_table))
##   KO1   KO2   KO3   KO4   KO5   KO6 
## 32859 35897 38718 37755 38827 36790
# Rarefaction Species Richness
set.seed(13)
otu_rare <- t(rrarefy(t(otu_table),  min(colSums(otu_table))))
head(colSums(otu_rare))
##   KO1   KO2   KO3   KO4   KO5   KO6 
## 32859 32859 32859 32859 32859 32859

物種注釋表處理

一般流程如Qiime2注釋獲得的taxonomy數(shù)據(jù)中往往含有P_、C_等字母，需要進(jìn)行分割和總計(jì)。

taxonomy <- read.delim("16S-amplicon-analysis/taxonomy.txt",  sep = "\t",  header = FALSE)

names(taxonomy) <- c("OTUs",  "taxo")

taxonomy$taxo <- gsub("[a-z]__", "", taxonomy$taxo)

taxo_names = c("kingdom", "phylum", "class", "order", "family"， "genus", "species")
clean_taxonomy <- tidyr::separate(taxonomy, col = taxo,  into = taxo_names,  sep = ";")

參考

Bukin, Yu S, Yu P Galachyants, IV Morozov, SV Bukin, AS Zakharenko, and TI Zemskaya. 2019. “The Effect of 16S rRNA Region Choice on Bacterial Community Metabarcoding Results.” Scientific Data 6 (1): 1–14.

Drengenes, Christine, Tomas ML Eagan, Ingvild Haaland, Harald G Wiker, and Rune Nielsen. 2021. “Exploring Protocol Bias in Airway Microbiome Studies: One Versus Two Pcr Steps and 16S rRNA Gene Region V3 V4 Versus V4.” BMC Genomics 22 (1): 1–15.