筆記內(nèi)容：
由二代測(cè)序產(chǎn)生的序列數(shù)據(jù)（fastq格式）到物種豐度的OTU table位迂， 樣本群落距離矩陣娄昆，物種多樣性指數(shù)晃痴，序列的進(jìn)化樹及物種注釋信息的分析結(jié)果，為常規(guī)分析流程缰趋∨跎迹可以使用usearch, vsearch, qiime等分析軟件實(shí)現(xiàn)，在必要的時(shí)候需要根據(jù)樣本信息的具體情況編寫腳本予以實(shí)現(xiàn)秘血。以下為雙端測(cè)序(pair-end)產(chǎn)生的fastq格式數(shù)據(jù)常規(guī)分析流程味抖，包括：
- fastqc及multiqc
- (qiime1)multiple_join_paired_ends.py
- fastq_screen
- (qiime1)multiple_split_libraries_fastq.py
- cutadapt
- usearch 流程
需要先了解barcode, primer的概念，可以通過下圖快速了解一下：barcode是每個(gè)樣本對(duì)應(yīng)的一段與樣本一一對(duì)應(yīng)的短序列灰粮，像超市貨物條形碼一樣區(qū)別每個(gè)樣本蔫敲。Fwd和Rev primer(引物)是用于16S擴(kuò)增與建庫的引物序列，在測(cè)序技術(shù)中起到定位的作用掌猛。有些下機(jī)數(shù)據(jù)會(huì)已經(jīng)將樣本按照barcode分好绰疤，每個(gè)樣本測(cè)得的序列片段都放在一個(gè)與樣本對(duì)應(yīng)的文件里。
input: fastq格式的下機(jī)數(shù)據(jù)（pair-end）:
已經(jīng)按照各樣本的barcode分好柑肴，格式為每個(gè)樣本兩個(gè)序列文件霞揉，比方說sampleID_1.fastq.gz, sampleID_2.fastq.gz。1為正向測(cè)序一次得到的結(jié)果晰骑，2為反向測(cè)序一次得到的結(jié)果适秩。這是pair-end雙向測(cè)序的原理，所以每個(gè)樣本得到兩個(gè)序列文件。
（注：標(biāo)注了(qiime1)表明需要在qiime環(huán)境中運(yùn)行秽荞。）
-
fastqc及multiqc
可以先把raw data骤公，即下機(jī)數(shù)據(jù)做一遍fastqc以了解其質(zhì)量情況。在安裝好fastqc, multiQC之后扬跋，在終端輸入：
$ fastqc -c dir/*.fastq.gz -o fastqout/
即在fastqout/文件夾中得到所有樣本的測(cè)序質(zhì)量報(bào)告阶捆。可以將這些報(bào)告整合成一個(gè)大的report钦听，為.html格式趁猴，雙擊在瀏覽器內(nèi)打開即可。
$ multiqc fastqout/ - (qiime1)multiple_join_paired_ends.py
把正反兩次測(cè)序結(jié)果join在一起彪见，正常的運(yùn)行結(jié)果應(yīng)為每個(gè)樣本一個(gè)文件夾儡司，文件夾內(nèi)為三個(gè)文件，一個(gè)是該樣本join到一起的序列文件（fastqjoin.join.fastq）余指，另外兩個(gè)為1捕犬，2兩個(gè)序列中無法join到一起的序列文件（fastqjoin.un1.fastq, fastqjoin.un2.fastq）。
multiple_join_paired_end.py為join_paired_ends.py的多個(gè)input文件用法酵镜，在qiime環(huán)境下的linux終端中使用碉碉。可以設(shè)置各種參數(shù)調(diào)整淮韭，比如常用--read1_indicator和--read1_indicator用于根據(jù)文件名稱定位哪個(gè)fastq文件為read1和read2垢粮。default設(shè)置為'_R1_'和'_R2_'。其他參數(shù)可詳見官網(wǎng)靠粪。
介于運(yùn)行結(jié)果為下圖左邊所示蜡吧，每個(gè)sample對(duì)應(yīng)一個(gè)文件夾，儲(chǔ)存了上述三個(gè)文件占键，且文件名沒有改變昔善。這一步可以用python寫個(gè)腳本用于文件名的調(diào)整，去掉沒有join到一起去的兩個(gè)文件畔乙，及統(tǒng)計(jì)每個(gè)樣本join起來的reads數(shù)占總reads數(shù)的比例君仆，以了解雙端測(cè)序join的情況。
fastq_screen
fastq_screen是可以用于檢查是否存在污染序列的軟件牲距，需要先install再在終端中使用返咱。FastQ Screen Documentation非常詳細(xì)。在終端中使用也不難牍鞠，input為上一步multiple_join后咖摹，修改了文件名的.fastq文件。
（需要將包含了污染序列的tag的文件去除皮服，留下screen之后的文件楞艾，自己寫腳本實(shí)現(xiàn)）-
(qiime1)multiple_split_libraries_fastq.py
multiple_split_libraries的目的為把上一步得到的一系列，每個(gè)樣本的fastq文件龄广，整合成一個(gè)大fastq文件硫眯，且每行開頭以“@樣本ID”標(biāo)注，如下所示择同。紅色虛線框里是sample1的第一個(gè)序列片段两入，第一行為樣本ID及第幾條序列，以及一些其他的敲才，后續(xù)可以去除的信息裹纳，接下來就是序列數(shù)據(jù)。加號(hào)單獨(dú)占一行紧武，接下來為一些質(zhì)量信息剃氧，以上也為fastq文件的格式的一些說明。
如下圖所示阻星，樣本ID整合到fastq文件中的結(jié)果：
cutadapt
使用cutadapt去掉序列頭尾的引物（primer）
測(cè)序公司會(huì)把頭尾兩段引物序列給你朋鞍。cutadapt有兩個(gè)用于指定這段引物序列的參數(shù)-a和-g。你可以先在終端用grep和less查看一下手頭序列信息有沒有去掉引物序列妥箕。grep的簡單用法可以參考答生信從業(yè)人員的N個(gè)問題（更新非常慢): 亂入一個(gè)grep怎么用
grep -n 'ATGCATGC' seqs.fastq把含有primer的行顯示出來
grep -c 'ATGCATGC' seqs.fastq顯示有多少行有primer
如果顯示出來大量序列在頭尾均包含primer序列滥酥，那說明primer應(yīng)該沒被去掉。同理你可以在去掉之后這樣看看畦幢，到底有沒有去掉坎吻。在grep -n的時(shí)候你可以留意一下primer是出現(xiàn)在頭部還是尾部。全部出現(xiàn)在頭部則把這段序列指定在-g之后宇葱，在尾部則指定在-a之后瘦真。注意這里不能反了。這兩個(gè)參數(shù)是用于決定把這段序列之前還是之后的序列切掉的黍瞧。比方說把尾部的序列填充到了-g, 則把尾部之前的序列都切掉吗氏，那output就沒有東西了。
其他cutadapt參數(shù)為質(zhì)量控制的參數(shù)雷逆，比如設(shè)置一個(gè)長度閾值弦讽，在切掉primer后如果序列長度小于這個(gè)值則扔掉這個(gè)序列等。詳細(xì)可見官網(wǎng)文件膀哲。
另外primer中可能使用了簡并堿基往产，比如M，R某宪，K仿村，Y等這樣不是ATGC四個(gè)密碼子的值。在終端grep的時(shí)候需要使用正則把他們對(duì)應(yīng)的ATGC寫出來兴喂，但是在cutadapt中不用蔼囊，直接填充M焚志，R等字母即可。usearch流程
fastq_filter
這個(gè)命令是一個(gè)序列質(zhì)量控制的過程畏鼓，使用-fastq_maxee 來設(shè)置每個(gè)reads所引起的expected errors酱酬。比方說設(shè)置成0.5，則會(huì)刪除expected errors大于0.5的reads.
bmp-Qiime2Uparse.pl
這是BMP的腳本云矫，該腳本的目的為將上一步fastq_filter得到的fasta文件header重新格式化膳沽，以便于后續(xù)的分析。
fastx_uniques
fastx_uniques是一個(gè)查找重復(fù)序列的過程让禀，將重復(fù)序列歸攏在一起挑社，從而只顯示unique序列（里面同時(shí)包含了每條unique序列出現(xiàn)的次數(shù)信息）。
-fastout 則output為fasta文件巡揍，且以各unique序列豐度降序排列痛阻。
-sizeout則在output文件的header中加上其重復(fù)reads的數(shù)目。
cluster_otus
這一步將上述得到的unique進(jìn)行聚類腮敌，按照97%相似度聚為各OTU录平。并輸出每個(gè)OTU的代表序列，以及每個(gè)OTU包含的序列數(shù)目缀皱。
Pick OTU一般有三種策略斗这，詳見16s微生物組第十條，這里usearch用的是UPARSE算法啤斗，本質(zhì)上應(yīng)該是一種denove的策略表箭。
這一步還去除了chimera(嵌合體序列)。 chimera是因?yàn)樵跀U(kuò)增中出現(xiàn)了問題钮莲，導(dǎo)致兩個(gè)來自不同DNA模板的序列被擴(kuò)增為同一條序列免钻。可以在一個(gè)chimera reference數(shù)據(jù)庫中找到相關(guān)序列的chimera并去除崔拥。
以下為一個(gè)從fastq_filter到cluster_otus的過程中极舔，樣本與序列信息的結(jié)果小結(jié)：
otutab
這一步將根據(jù)前面qiime中得到的seqs.fastq，及OTU代表序列链瓦，將各個(gè)樣本的序列與OTU代表序列相比對(duì)拆魏，生成OTU table。
-sample_delim 是一個(gè)指定樣本ID和序列ID分隔符的參數(shù)慈俯，比方說qiime生成的seq.fastq是以“sample1_0”,“sample1_1”命名的渤刃，_ 之前的內(nèi)容才是樣本名稱。則可以把sample_delim指定為 _贴膘。
-otutabout
這一步得到了最初的OTU table卖子，即每個(gè)sample對(duì)應(yīng)在每個(gè)OTU中有多少個(gè)reads。接下來可以查看一下樣本中OTU的豐度如何刑峡，都分布在什么范圍內(nèi)洋闽。
biom summarize-table -i otu_raw.tab
隨后可以用otutab_trim來對(duì)raw OTU table進(jìn)行過濾玄柠。
otutab_trim
設(shè)置參數(shù)用于去除OTU table中豐度太低的樣本或者OTU。
-min_sample_size 5000
即去除reads數(shù)小于5000的樣本
-min_otu_size 2
去除reads數(shù)小于2的OTU
-min_otu_freq 0.001
去除 單個(gè)OTUreads數(shù)/總OTUreads數(shù)目 小于0.001的OTU
otutab_norm
指定參數(shù)-sample_size诫舅，將每個(gè)樣本的reads數(shù)目抽平到指定的reads數(shù)目羽利。比如 -sample_size 5000即把每個(gè)樣本的reads總數(shù)目抽平到5000.這一步不涉及對(duì)樣本或者OTU的刪除或規(guī)整，只是歸一化骚勘。
-cluster_agg
以O(shè)TU代表序列的fasta格式文件為Input铐伴，生成otu.tree(OTU的進(jìn)化樹)撮奏。
-alpha_div
以trim及norm后的OTU table為input, 生成alpha diversity的各種多樣性指數(shù)俏讹。
-beta_div
以otu.tree及OTU table作為input, 生成beta diversity中的各種距離矩陣。注意想要生成Unifrac距離矩陣必須用指定-tree參數(shù)(tree file)畜吊。
-sintax
這一步用于根據(jù)各OTU的代表序列進(jìn)行物種預(yù)測(cè)泽疆，并給每個(gè)OTU添加物種注釋。以O(shè)TU代表序列的fasta文件及物種參考數(shù)據(jù)庫為input玲献。參考的物種數(shù)據(jù)庫一般為一個(gè)包含詳細(xì)物種信息的fasta格式的文件殉疼。usearch官網(wǎng)上推薦使用一個(gè)小而權(quán)威的參考數(shù)據(jù)庫，可以參考RDP或SILVA LTP等16S數(shù)據(jù)庫捌年。
這一步得到的output為以下格式的文檔：
| OTU_ID | tax_with_condifence | strand | tax |
|---|---|---|---|
| OTU1 | d:Bacteria(1.0000),p:Firmicutes(0.9300),c:Clostridia(0.8400),o:Clostridiales(0.8400),f:Ruminococcaceae(0.3200),g:Clostridium_III(0.2400) | - | d:Bacteria,p:Firmicutes,c:Clostridia,o:Clostridiales |
第二列中括號(hào)里為該物種注釋的置信度瓢娜，第三列為strand, 即正負(fù)鏈。一般為全是正鏈(+)或全是負(fù)鏈(-)礼预。第四列為根據(jù)-sintax_cutoff參數(shù)設(shè)置界值眠砾，比如設(shè)置為0.8則只納入第二列中置信度大于0.8的物種注釋。
經(jīng)過上述步驟托酸，得到trim及標(biāo)準(zhǔn)化后的OTU table, otu.tree, alpha多樣性指數(shù)文件褒颈，beta diversity各距離矩陣及物種注釋信息sintax.txt。后續(xù)分析參考微生物組16S rRNA數(shù)據(jù)分析小結(jié)：從OTU table到marker物種
附：本流程中使用version10的usearch, qiime的版本參數(shù)如下(在qiime1環(huán)境中使用print_qiime_config.py查看版本信息)
