一. 16S-seq實驗
所用示例數(shù)據(jù)使用Illumina公司兩步PCR的方法進(jìn)行擴(kuò)增：
1st stage PCR, 使用含adaptor的16S rRNA基因的通用引物，比如505F-80R，799F-1192等烹困；
2nd stage PCR, 使用含P5/P7蠢莺，index雳锋，和adaptor的引物，給每個樣品的序列加上特定的index
測序使用2x250bp雙端測序趁桃。
二. 安裝分析所需要的軟件
安裝R語言
(1)?RStudio | Open source & professional software for data science teams - RStudio
(2) R: The R Project for Statistical Computing (r-project.org)
安裝分析包仑最。
(1) R含有大量的分析包扔役，統(tǒng)計分析、畫圖等packages可以直接通過R.studio菜單中的tools>install packages安裝警医；或者使用以下命令：
install.packages("pkg-name")
(2) R中生物信息學(xué)常用的分析包來自bioconductor亿胸，需要先在R中安裝BiocManager，再使用BiocManager安裝分析包
install.packages("BiocManager")? ? # 安裝BiocManager軟件包
?BiocManager::install("dada2") ? ? ? # 安裝bioconductor中的dada2分析包
···
安裝中可能會出現(xiàn)以下提示：
可輸入a法严，表示update所有與dada2有關(guān)聯(lián)的分析包损敷。
(3) 本示例中需要安裝以下軟件包：dada2葫笼，phyloseq深啤，ggplot2，DESeq2
(4) 分析開始前路星，在Rstudio讀入所需的軟件包
(5) 設(shè)置好工作目錄：通過Rstudio的session>Set Working Directory>choose Directory設(shè)置
或者直接在終端中輸入以下命令：
setwd("E:/teaching/da2-pipline")? #中間的路徑需要根據(jù)自己的電腦修改
三溯街、測序文件的處理
1. 首先查看以下收到的測序文件，常用的2x250bp測序結(jié)果會含兩個文件洋丐，比如sample-R1.fq和sample-R2.fq呈昔，分別用R1和R2表示5'-和3'-測序結(jié)果，sample為原理提供的樣品名友绝。
2. 將測序文件的信息讀入R語言中堤尾。
以下示例中，現(xiàn)將測序結(jié)果所在的文件夾路徑讀入到path迁客，然后將該文件夾下面所有以“-R1.fq”
#設(shè)置path保存測序文件所在的路徑
path <- "Z:/dada2_pipline";
#以下命令會顯示path保存的文件夾路徑下所有文件名
list.files(path);
#以下命令顯示path文件夾中郭宝，所有以-R1.fq結(jié)尾的文件，并排序
sort(list.files(path, pattern="-R1.fq", full.names = TRUE))
# 明白以上命令后掷漱，可運(yùn)行下面的代碼粘室，將測序結(jié)果文件信息讀入fnFs和fnRs
fnFs <- sort(list.files(path, pattern="-R1.fq", full.names = TRUE))?
fnRs <- sort(list.files(path, pattern="-R2.fq", full.names = TRUE))?
# 在Rstudio，可查看右上角中fnFs和FnRs的內(nèi)容卜范，是否為正確的測序文件名
#或者直接在Rstudio的控制終端衔统，輸入fnFs查看
3. 查看測序質(zhì)量。
#查看測序質(zhì)量
#繪制樣品1-2測序質(zhì)量的圖片锦爵，并保存到plot_output
plot_output<-plotQualityProfile(fnFs[1:2]);
#查看圖片, 注意右側(cè)Plots窗口
plot_output;? ? ?
#繪制所有樣品測序質(zhì)量的圖
plot_output<-plotQualityProfile(fnFs);? ?
plot_output;
#保存圖片到當(dāng)前目錄
ggsave(filename="For_seq_qual.pdf",plot_output,
? ? ? width = 18, height = 6,
? ? ? units = "cm");
plotQualityProfile(fnRs);? #? 繪制Rev測序的質(zhì)量圖片舱殿，可按上述代碼保存圖片
#移除暫時不用的plot_output
rm(plot_output);
關(guān)于測序質(zhì)量圖的注解
In gray-scale is a heat map of the frequency of each quality score at each base position. The mean quality score at each position is shown by the green line, and the quartiles of the quality score distribution by the orange lines. The red line shows the scaled proportion of reads that extend to at least that position (this is more useful for other sequencing technologies, as Illumina reads are typically all the same length, hence the flat red line).
4. 過濾接頭和低質(zhì)量序列
從上述的圖中怀薛，以Q30為閾值，可以看出迷郑，F(xiàn)orward（R1）的整體質(zhì)量都還可以枝恋；Reverse（R2）的測序質(zhì)量稍低，平均來看嗡害，大概也在200bp左右出現(xiàn)了下降焚碌。
設(shè)置過濾后文件保存的位置
此示例中，是在原來path保存的目錄下面霸妹，建立一個“filtered”的文件夾十电，用來保存過濾后的序列文件，文件名按照sample_F_filt.fastq.gz的模式命名叹螟，其中后面的.gz表示該文件是壓縮文件鹃骂。
#paste0()函數(shù)將sample.names中保存的樣品名（即原來測序文件的前半部分）加上"_F_filt.fastq.gz"
paste0(sample.names, "_F_filt.fastq.gz");
#file.path()函數(shù)，將目錄path罢绽，“filtered”畏线，文件名用/連起來，
filtFs <- file.path(path, "filtered", paste0(sample.names, "_F_filt.fastq.gz"));
filtRs <- file.path(path, "filtered", paste0(sample.names, "_R_filt.fastq.gz"));
運(yùn)行以下filterAndTrim函數(shù)良价，過濾低質(zhì)量序列寝殴。out保存過濾信息。
在該示例中明垢，
fnFs和fnRs為原始序列文件蚣常，filtFs和filtRs為過濾后的序列文件；
truncLen=c(208,200)表示分別去掉208（For測序結(jié)果）和200（rev測序結(jié)果）以后的堿基;
MaxN=0痊银，切除末端的低質(zhì)量序列后抵蚊，得的的序列中不能再有N；
maxEE=c(2,2)溯革，序列的EE值不能低于maxEE贞绳，EE = sum(10^(-Q/10)). 比如200bp的序列，每個堿基的Q=30鬓照，EE=0.2.
rm.phix=TRUE熔酷，去除其中來自phix噬菌體的序列（測序反應(yīng)中加入phix DNA作為對照）；
compress=TRUE豺裆，壓縮輸出的序列文件拒秘。
如果需要切除5'端的引物序列号显，可以用加入trimLeft=c(N1,N2)。
out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(208,200),
? ? ? ? ? ? ? ? ? ? maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
? ? ? ? ? ? ? ? ? ? compress=TRUE, multithread=FALSE);
#查看過濾信息
head(out);
#查看過濾后的序列質(zhì)量?
plotQualityProfile(filtFs);?
plotQualityProfile(filtRs);
5. 推斷序列變異（Infer Sequence variant）
DADA2中躺酒，learnErrors程序“學(xué)習(xí)”輸入數(shù)據(jù)中測序產(chǎn)生錯誤押蚤，并構(gòu)建預(yù)測模型。
# 這里是DADA2比較神奇的地方, 它會從輸入的數(shù)據(jù)中學(xué)習(xí)測序的誤差羹应。理論上揽碘，輸入的數(shù)據(jù)越大，
# 學(xué)習(xí)模型越準(zhǔn)確园匹，需要的計算資源也越大
errF <- learnErrors(filtFs, multithread=FALSE);
errR <- learnErrors(filtRs, multithread=FALSE);
#繪制quality score與error frequency相關(guān)圖形
plotErrors(errF, nominalQ=TRUE);
plotErrors(errF, nominalQ=TRUE);
然后利用該模型，去除測序產(chǎn)生的堿基替換裸违、indel錯誤掖桦，來獲得“真實”的序列（denoised sequence）。
> dadaFs <- dada(filtFs, err=errF, multithread=FALSE)
Sample 1 - 8382 reads in 4284 unique sequences.
Sample 2 - 9341 reads in 4959 unique sequences.
Sample 3 - 8866 reads in 4794 unique sequences.
Sample 4 - 8097 reads in 4720 unique sequences.
> dadaRs <- dada(filtRs, err=errR, multithread=FALSE)
Sample 1 - 8382 reads in 6273 unique sequences.
Sample 2 - 9341 reads in 5980 unique sequences.
Sample 3 - 8866 reads in 6526 unique sequences.
Sample 4 - 8097 reads in 6833 unique sequences.
將雙端測序的reads合并
minOverlap
#將雙端測序的reads合并
mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs)
構(gòu)建ASV table
> #構(gòu)建sequence table供汛，即OTU table或ASV
> #包含ASV枪汪，每一個ASV
> seqtab <- makeSequenceTable(mergers);
#表格的行、列數(shù)目怔昨。本示例中雀久，有4行（4個樣品）；587列（587條ASV）趁舀；
> dim(seqtab);
[1]? ?4 587
> seqtab[,1];
root_exp1 root_exp2 soil_exp1 soil_exp2?
? ? ? 522? ? ? ?364? ? ? ? ?0? ? ? ? ?0?
> #ASV的長度信息統(tǒng)計
> table(nchar(getSequences(seqtab)));
去嵌合體（remove chimeras）
## 去PCR產(chǎn)生的嵌合體
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE);
dim(seqtab.nochim);
?#計算每個樣品中嵌合體的比例
apply(seqtab.nochim,1,sum)/apply(seqtab,1,sum);
#查看第一條ASV的序列（即第一列的標(biāo)題）
colnames(seqtab)[1]
增加序列的物種信息
##給ASV序列賦予物種分類信息
## 分類數(shù)據(jù)庫保存在當(dāng)前目錄的tax文件夾中(“./tax”)
## 如內(nèi)存不夠赖捌，可刪除中間信息
rm(errF,errR,plot_output, seqtab,dadaFs,dadaRs,derepFs,derepRs);
taxa <- assignTaxonomy(seqtab.nochim, "./tax/silva_nr_v132_train_set.fa.gz", multithread=TRUE)
#增加species的信息。
taxa <- addSpecies(taxa, "./tax/silva_species_assignment_v132.fa.gz")
到此赫编，序列的處理結(jié)束巡蘸，接下來用phyloseq包分析微生物多樣性信息奋隶。
四. 多樣性分析
1. 將序列處理結(jié)果phyloseq的格式擂送。
#構(gòu)建樣品信息表格
#本示例中。有兩個root和soil兩種材料唯欣，每個兩個重復(fù)嘹吨。
#查看sample names
sample.names;
samdf <- data.frame(rownames(seqtab.nochim),c("root","root","soil","soil"));
rownames(samdf)<- sample.names;
colnames(samdf)<-c("id","group")
samdf;
2. 整理OTU table和taxa table蟀拷，主要是將原來序列作為列標(biāo)題的表格換為人工設(shè)定的OTU+數(shù)字
##整理OTU table（或ASV table）
# 上述腳本輸出的seqtab.nochim表格中，列的標(biāo)題是ASV的序列（非常長）萍聊，在這里用簡單的OTU+數(shù)字來代替
seqtab.nochim0<-seqtab.nochim;
colnames(seqtab.nochim0)<-paste(rep("OTU",ncol(seqtab.nochim)),1:ncol(seqtab.nochim),sep="");
## 同上问芬，將taxa表格中的分類信息用OTU+數(shù)據(jù)來代替。
taxa0<-taxa
rownames(taxa0)<-paste(rep("OTU",ncol(seqtab.nochim)),1:ncol(seqtab.nochim),sep="");
head(taxa0);
## 將信息讀入phyloseq所需的格式
ps <- phyloseq(otu_table(seqtab.nochim0, taxa_are_rows=FALSE),? ?sample_data(samdf),? ? ? ? ? ? ? ? tax_table(taxa0));
3. 獲取多樣性信息
（1）alpha多樣性
#可以開始查看多樣性的信息
# 查看alpha多樣性
estimate_richness(ps, measures =c("Observed", "Chao1", "ACE", "Shannon", "Simpson", "InvSimpson", "Fisher"));
#繪制alpha多樣性的比較圖
plot_richness(ps);
（2）對數(shù)據(jù)進(jìn)行rarefaction平衡化后重新分析寿桨。
##rarefaction后重新估計多樣性數(shù)據(jù)
set.seed(7)
#參考每個樣品的序列總數(shù)
rowSums(otu_table(ps));
#最小值為3993此衅，可將sample.size設(shè)為3993以下强戴，進(jìn)行rarefaction
ps.rarefy<-rarefy_even_depth(ps, sample.size=3900, replace=T);
#利用rarefactio的數(shù)據(jù)重新估計樣品多樣性
# 查看alpha多樣性
estimate_richness(ps.rarefy, measures =c("Observed", "Chao1", "ACE", "Shannon", "Simpson", "InvSimpson", "Fisher"));
#繪制alpha多樣性的比較圖
plot_richness(ps);
#繪制熱圖
plot_heatmap(ps.rarefy, method = "NMDS", distance = "bray");
比較兩次的輸出結(jié)果
（3）beta多樣性比較
繪制heatmap
#繪制熱圖
plot_heatmap(ps.rarefy, method = "NMDS", distance = "bray");
PCA分析
#PCA分析
pslog<-transform_sample_counts(ps.rarefy, function(x) log(1+x));
out.log<-ordinate(pslog, method="NMDS", distance="bray");
plot_ordination(pslog, out.log, color="id", shape="group");
細(xì)菌組成
#繪制每個樣品門水平的豐度信息
plot_bar(ps.rarefy,x="id", fill="Phylum");
#繪制豐度高的OTU代表的細(xì)菌Family
top20 <- names(sort(taxa_sums(ps.rarefy), decreasing=TRUE))[1:20];
ps.top20 <- transform_sample_counts(ps.rarefy, function(OTU) OTU/sum(OTU));
ps.top20 <- prune_taxa(top20, ps.top20);
plot_bar(ps.top20, x="id", fill="Family") + facet_wrap(~group, scales="free_x");
#繪制豐度高的OTU代表的細(xì)菌Family
top20 <- names(sort(taxa_sums(ps.rarefy), decreasing=TRUE))[1:20];
ps.top20 <- transform_sample_counts(ps.rarefy, function(OTU) OTU/sum(OTU));
ps.top20 <- prune_taxa(top20, ps.top20);
plot_bar(ps.top20, x="id", fill="Family") + facet_wrap(~group, scales="free_x");
