易擴增子：易用、可重復和跨平臺的擴增子分析流程

1、 /exxxx DOI:10.21769/BioProtoc.xxxx55555111112000易擴增子：易用、可重復和跨平臺的擴增子分析流程EasyAmplicon：An easy-to-use, reproducible and cross-platform pipeline for amplicon analysis劉永鑫1, 2, 3, #, *，陳同4#，周欣5，白洋1, 2, 3, 6, *1中國科學院遺傳與發(fā)育生物學研究所，植物基因組學國家重點實驗室，北京；2中國科學院大學，生物互作卓越創(chuàng)新中心，北京；3中國科學院遺傳與發(fā)育生物學研究所，中國科學院英國約翰英納斯中心植物和微生物

2、科學聯(lián)合研究中心，北京；4中國中醫(yī)科學院，中藥資源中心，北京；5中國科學院微生物研究所，真菌學國家重點實驗室，北京；6中國科學院大學現(xiàn)代農(nóng)學院，北京*通訊作者郵箱: HYPERLINK mailto:yxliu yxliu, HYPERLINK mailto:ybai ybai#共同第一作者/同等貢獻摘要擴增子測序是目前微生物組研究中最廣泛的使用手段，主流的分析流程有QIIME、USEARCH和Mothur，但仍分別存在依賴關(guān)系過多導致的安裝困難、大數(shù)據(jù)收費和使用界面不友好等問題。本文搭建了一套完整的擴增子分析流程命名為易擴增子(EasyAmplicon)，可以實現(xiàn)簡單易用、可重復和跨平臺地開

3、展擴增子分析。流程計算核心采用體積小、安裝方便、計算速度快且跨平臺的軟件USEARCH，同時整合VSEARCH以突破USEARCH免費版的限制。分析選用RStudio的圖形界面對流程代碼文檔管理和運行，實現(xiàn)命令行和/或鼠標點擊操作方式開展擴增子可重復分析。同時流程還提供了數(shù)10個腳本，實現(xiàn)特征表過濾、重采樣、分組均值等常計算，以及為常用軟件如STAMP、LEfSe、PICRUSt1/2等提供標準的輸入文件的功能。易擴增子可以從 HYPERLINK /YongxinLiu/EasyAmplicon /YongxinLiu/EasyAmplicon 下載輕松部署于Windows/Mac/Linux

4、操作系統(tǒng)，在普通個人電腦上2小時內(nèi)可完成數(shù)十個樣本的分析。關(guān)鍵詞: 擴增子，分析流程，16S，ITS，USEARCH儀器設(shè)備個人電腦/服務(wù)器 (操作系統(tǒng)：Windows 10 / Mac OS 10.12+ / Linux Ubuntu 18.04+；CPU：2核+；內(nèi)存：8G+；硬盤： 10 GB，且大于10倍原始數(shù)據(jù)大小)，網(wǎng)絡(luò)訪問暢通。軟件和數(shù)據(jù)庫R語言環(huán)境，下載適合自己系統(tǒng)的4.0.2版： HYPERLINK / / R語言開發(fā)環(huán)境，用于執(zhí)行流程，下載適合自己系統(tǒng)的RStudio 1.3.1056： HYPERLINK /products/rstudio/download/ l dow

5、nload /products/rstudio/download/#download(可選)僅Windows系統(tǒng)安裝，提供Git Bash命令行環(huán)境的GitForWidnows 2.28.0： HYPERLINK / / 擴增子分析流程USEARCH v10.0.240 ADDIN EN.CITE Edgar2010172(Edgar 2010)17217217Edgar, Robert C.Search and clustering orders of magnitude faster than BLASTBioinformaticsBioinformaticsBioinformaticsB

6、ioinformatics2460-2461261920101367-4803/10.1093/bioinformatics/btq46110.1093/bioinformatics/btq4615/28/2019(Edgar 2010) HYPERLINK /usearch/download.html /usearch/download.html擴增子分析流程VSEARCH v2.15.0 ADDIN EN.CITE Rognes2016171(Rognes等， 2016)17117117Rognes, TorbjrnFlouri, TomNichols, BenQuince, Christ

7、opherMah, FrdricHrbek, TomasVSEARCH: a versatile open source tool for metagenomicsPeerJPeerJPeerJPeerJPeerJPeerJPeerJPeerJe258442016/10/27ClusteringChimera detectionSearchingMaskingShufflingParallellizationMetagenomicsAlignmentSequencesDereplication20162016/10/182167-835927781170/10.7717/peerj.2584P

8、MC507569710.7717/peerj.2584(Rognes等， 2016) HYPERLINK /torognes/vsearch/releases /torognes/vsearch/releases易擴增子流程EasyAmplicon v1.10 ADDIN EN.CITE ADDIN EN.CITE.DATA (Zhang等， 2018; Chen等， 2019; Huang等， 2019; Zhang等， 2019; Liu等， 2020; Qian等， 2020; Qian等， 2020)： HYPERLINK /YongxinLiu/EasyAmplicon /Yongx

9、inLiu/EasyAmplicon 核糖體數(shù)據(jù)庫RDP v16 ADDIN EN.CITE Cole201490(Cole等， 2014)9090017Cole, James R.Wang, QiongFish, Jordan A.Chai, BenliMcGarrell, Donna M.Sun, YanniBrown, C. TitusPorras-Alfaro, AndreaKuske, Cheryl R.Tiedje, James M.Ribosomal Database Project: data and tools for high throughput rRNA analysi

10、sNucleic Acids ResearchNucleic Acids ResearchNucleic Acids Res.Nucleic Acids ResD633-D64242D120140305-1048/10.1093/nar/gkt124410.1093/nar/gkt12445/6/2019(Cole等， 2014)： HYPERLINK / / (可選)核糖體數(shù)據(jù)庫GreenGene數(shù)據(jù)庫(gg) 13_8 ADDIN EN.CITE McDonald2011305(McDonald等， 2011)30530517McDonald, DanielPrice, Morgan N.

11、Goodrich, JuliaNawrocki, Eric P.DeSantis, Todd Z.Probst, AlexanderAndersen, Gary L.Knight, RobHugenholtz, PhilipAn improved Greengenes taxonomy with explicit ranks for ecological and evolutionary analyses of bacteria and archaeaThe ISME JournalThe ISME JournalISME J.ISME J6106201112/01/onlineThe Aut

12、hor(s)Original Article/10.1038/ismej.2011.13910.1038/ismej.2011.139(McDonald等， 2011)： HYPERLINK ftp:/greengenes.microbio.me/greengenes_release ftp:/greengenes.microbio.me/greengenes_release (可選)核糖體數(shù)據(jù)庫SILVA v123 ADDIN EN.CITE ADDIN EN.CITE.DATA (Quast等， 2013)： HYPERLINK http:/www.arb-silva.de http:/w

13、ww.arb-silva.de (可選)轉(zhuǎn)錄間隔區(qū)(ITS)數(shù)據(jù)庫UNITE v8.2 ADDIN EN.CITE Nilsson201993(Nilsson等， 2019)9393017Nilsson, Rolf HenrikLarsson, Karl-HenrikTaylor, Andy FSBengtsson-Palme, JohanJeppesen, Thomas SSchigel, DmitryKennedy, PeterPicard, KathrynGlckner, Frank OliverTedersoo, LehoSaar, IrjaKljalg, UrmasAbarenkov

14、, KessyThe UNITE database for molecular identification of fungi: handling dark taxa and parallel taxonomic classificationsNucleic Acids ResearchNucleic Acids ResearchNucleic Acids Res.Nucleic Acids ResD259-D26447D120190305-1048/10.1093/nar/gky102210.1093/nar/gky10225/7/2019(Nilsson等， 2019)： HYPERLIN

15、K https:/unite.ut.ee/ https:/unite.ut.ee/ (可選)Windows版下載工具wget： HYPERLINK /packages/wget.htm /packages/wget.htm 軟件安裝和數(shù)據(jù)庫部署注：以下的軟件安裝和使用均在64位Windows 10系統(tǒng)中演示，Linux/Mac中不同的地方會有說明，流程代碼提供有Mac版本(pipeline_mac.sh)。Windows系統(tǒng)需要安裝GitForWindows( HYPERLINK / /)提供Git bash環(huán)境支持常用Shell命令。Linux/Mac系統(tǒng)自帶Bash命令行工作環(huán)境。以64位

16、Windows 10系統(tǒng)為例，我們先安裝R/RStudio軟件，再把本流程(EasyAmplicon/目錄)保存于C盤中，然后根據(jù)需要下載數(shù)據(jù)庫至指定目錄即完成部署。注：代碼行添加灰色底紋背景，其中需要根據(jù)系統(tǒng)環(huán)境修改的部分標為藍色。流程運行環(huán)境R和RStudio依次安裝適合系統(tǒng)的最新版R語言( HYPERLINK )和RStudio( HYPERLINK /products/rstudio/download/ l download /products/rstudio/download/)。注意操作系統(tǒng)用戶名不要使用中文，否則會影響R語言使用。批量安裝依賴R包流程會調(diào)用數(shù)百個R包，使用時可自動

17、安裝。但由于網(wǎng)絡(luò)或系統(tǒng)等個性原因經(jīng)常出現(xiàn)下載或安裝失敗，可以使用中根據(jù)提示手動安裝缺失R包。本文推薦直接下載我們預(yù)編譯好的R包合輯( HYPERLINK /datadownload /datadownload)，替換至R包所在目錄即可，詳見常見問題1。易擴增子流程EasyAmplicon訪問 HYPERLINK /YongxinLiu/EasyAmplicon /YongxinLiu/EasyAmplicon，選擇Code Download ZIP下載并解壓，如保存于C盤并確保目錄名為EasyAmplicon。如在RStudio的Terminal中可使用git下載流程：git clone gi

18、t:YongxinLiu/EasyAmplicon.git(可選)擴增子流程依賴軟件EasyAmplicon依賴的Windows/Mac/Linux版本軟件已經(jīng)保存于EasyAmplicon中的win/mac/linux目錄中，如果出現(xiàn)問題，可按如下方法手動安裝。USEARCH下載頁 HYPERLINK /usearch/download.html /usearch/download.html，選擇適合自己系統(tǒng)的10.0.240版本 (不要下載最新版，因為有更多功能使用受限)，如Windows版本保存至EasyAmplicon目錄中的win目錄中，解壓后改名為usearch.exe。Linux

19、/Mac系統(tǒng)需下載到環(huán)境linux/mac目錄，解壓后改名為usearch，并添加可執(zhí)行權(quán)限(chmod +x usearch)。VSEARCH下載頁面 HYPERLINK /torognes/vsearch/releases /torognes/vsearch/releases，選擇適合自己系統(tǒng)的最新版下載，接下來操作與USEARCH類似。Windows系統(tǒng)還需下載wget程序( HYPERLINK /packages/wget.htm /packages/wget.htm)至win目錄。(可選)16S核糖體基因物種注釋數(shù)據(jù)庫16S擴增子測序分析，常用RDP/SILVA/GreenGene數(shù)

20、據(jù)庫進行物種注釋，可以從上述數(shù)據(jù)庫官網(wǎng)下載并整理為USEARCH使用的格式，此處推薦從USEARCH官網(wǎng)( HYPERLINK /sintax /sintax)下載USEARCH兼容格式的數(shù)據(jù)庫。默認流程使用體積小巧的RDP v16數(shù)據(jù)庫(rdp_16s_v16_sp.fa.gz)，并已保存于usearch目錄中?？蛇xGreenGenes 13.5(gg_16s_13.5.fa.gz)和SILVA(silva_16s_v123.fa.gz)數(shù)據(jù)庫，可根據(jù)需要下載并保存于usearch目錄中。此外，如果要開展PICRUSt和Bugbase功能預(yù)測分析，還需要使用GreenGenes數(shù)據(jù)庫13.5

21、中按97%聚類的OTU序列(己保存于流程gg目錄中97_otus.fasta.gz)。可選手動下載GreenGenes官方數(shù)據(jù)庫( HYPERLINK ftp:/greengenes.microbio.me/greengenes_release ftp:/greengenes.microbio.me/greengenes_release)，解壓后選擇其中的97_otus.fasta保存于gg目錄下即可。(可選)ITS物種注釋數(shù)據(jù)庫如果研究真菌或真核生物采用轉(zhuǎn)錄間隔區(qū)(Intergenic Transcribed Spacer)測序，需要使用UNITE數(shù)據(jù)庫，目前最新版已經(jīng)保存于usearch目

22、錄(utax_reference_dataset_all_04.02.2020.fasta.gz)。如流程中數(shù)據(jù)庫沒有及時更新，可在UNITE官網(wǎng)( HYPERLINK https:/unite.ut.ee/ https:/unite.ut.ee/)下載最新版適合USEARCH的注釋數(shù)據(jù)庫。官方數(shù)據(jù)庫存在格式問題，詳細常見問題2。實驗步驟開始新項目分析前，我們需要在項目目錄(如c:/test)中準備三類起始文件：1. EasyAmplicon流程中復制分析流程文件(pipeline.sh)；2. 編寫樣本元數(shù)據(jù)(metadata.txt)；3. seq目錄存放測序數(shù)據(jù)(*.fq.gz)。然后使

23、用RStudio打開pipeline.sh，設(shè)置分析流程(EasyAmplicon, ea)和工作目錄(work directory, wd)位置，添加依賴可執(zhí)行程序至環(huán)境變量，并切換至工作目錄。注：用戶請根據(jù)操作系統(tǒng)類型、軟件和工作目錄實際位置自行修改。ea=/c/EasyAmpliconwd=/c/testPATH=$PATH:$ea/wincd $wd準備輸入數(shù)據(jù)(測試數(shù)據(jù)是流程正對照)建議下載測試元數(shù)據(jù)和測序數(shù)據(jù)作為實驗的正對照，首先完成全部流程分析，以確定流程部署成功。將來使用自己的數(shù)據(jù)出現(xiàn)問題，可以與測試數(shù)據(jù)分析中對應(yīng)結(jié)果比較，以便確定問題產(chǎn)生的原因。下載示例元數(shù)據(jù)用于參考編寫格式

24、(表1)。wget -c HYPERLINK 10/github/EasyAmplicon/data/metadata.txt 10/github/EasyAmplicon/data/metadata.txt表1.元數(shù)據(jù)格式示例SampleIDGroupDateSiteCRACRRKO1KO2017/6/30ChaoyangCRA002352 CRR117575KO2KO2017/6/30ChaoyangCRA002352 CRR117576KO3KO2017/7/2ChangpingCRA002352CRR117577可用Excel編寫，保存存為制表符分隔的的文本文件。注意：有行列標題，行為

25、樣品名(字母開頭+數(shù)字組合)，列為分組信息(至少1列，可多列)、地點和時間(提交數(shù)據(jù)必須)、及其它屬性，詳見附表1，或下載的metadatat.txt文件。通常測序公司會返回原始數(shù)據(jù)，如Illumina雙端測序的文件，每個樣本有一對文件。本文使用的數(shù)據(jù)來自發(fā)表于Science雜志關(guān)于擬南芥根系微生物組研究的文章 ADDIN EN.CITE Huang2019116(Huang等， 2019)116116017Huang, Ancheng C.Jiang, TingLiu, Yong-XinBai, Yue-ChenReed, JamesQu, BaoyuanGoossens, AlainNtz

26、mann, Hans-WilhelmBai, YangOsbourn, AnneA specialized metabolic network selectively modulates Arabidopsis root microbiotaScienceScienceScienceScienceeaau638936464402019/content/364/6440/eaau6389/content/sci/364/6440/eaau6389.full.pdf10.1126/science.aau6389(Huang等， 2019)，GSA項目號為PRJCA001296。為方便演示流程的使用

27、，我們從中選取三個組（每組包括6個生物學重復共18個樣本），并且隨機抽取了50000對序列作為軟件的測序數(shù)據(jù)，該數(shù)據(jù)可以從中國科學院基因組研究所的原始數(shù)據(jù)歸檔庫（Genome Sequence Archive，GSA， HYPERLINK /gsa/ /gsa/） ADDIN EN.CITE Wang2017632(Wang等， 2017)63263217Wang, YanqingSong, FuhaiZhu, JunweiZhang, SisiYang, YadongChen, TingtingTang, BixiaDong, LiliDing, NanZhang, QianBai, Zho

28、uxianDong, XunongChen, HuanxinSun, MingyuanZhai, ShuangSun, YubinYu, LeiLan, LiXiao, JingfaFang, XiangdongLei, HongxingZhang, ZhangZhao, WenmingGSA: Genome Sequence Archive*Genomics Proteomics BioinformaticsGenomics Proteomics BioinformaticsGenom. Proteom. Bioinf.Genom Proteom Bioinf14-18151Genome S

29、equence ArchiveGSABig dataRaw sequence dataINSDC20172017/02/01/1672-0229/science/article/pii/S167202291730002510.1016/j.gpb.2017.01.001(Wang等， 2017)中按批次編號CRA002352搜索并下載。本文使用wget根據(jù)樣本元數(shù)據(jù)中批次和樣本編號批量下載至seq目錄，代碼如下。mkdir -p seqawk system(wget -c /gsa/$5/$6/$6_f1.fq.gz -O seq/$1_1.fq.gz) (tail -n+2 metadata

30、.txt)awk system(wget -c /gsa/$5/$6/$6_r2.fq.gz -O seq/$1_2.fq.gz) temp/all.fq圖1. 典型擴增子測序雙端序列合并結(jié)構(gòu)圖引物切除和質(zhì)量控制采用等長的方式切除引物，引物外側(cè)如果有標簽(Barcode)，標簽的長度需要計算在內(nèi)(圖1)。如本示例的結(jié)構(gòu)為10 bp左端標簽 + 19 bp正向引物V5共計29 bp，18 bp反向引物V7，分別使用-fastq_stripleft和-fastq_stripright傳遞給程序。注意：務(wù)必清楚實驗設(shè)計中引物和標簽長度，如果引物已經(jīng)去除，可在下方參數(shù)處填0表示無需去除。此外，-fas

31、tq_maxee_rate指定質(zhì)量控制的錯誤率，0.01代表要求錯誤率小于1%。質(zhì)控控制后，-fastaout輸出體積更小的無質(zhì)量值fasta(fa)格式文件。vsearch -fastx_filter temp/all.fq -fastq_stripleft 29 -fastq_stripright 18 -fastq_maxee_rate 0.01 -fastaout temp/filtered.fa序列去冗余并挑選代表序列(OTU/ASV)序列去冗余可以使總數(shù)據(jù)量降低數(shù)量級(降維)，減少下游計算資源消耗和縮短等待時間，更重要的是提供序列的出現(xiàn)頻次對鑒定真實特征序列至關(guān)重要。擴增子分析中特

32、征序列有可操作分類單元(OTUs)或擴增序列變體(ASVs)。參數(shù)-miniuniqusize設(shè)置使用序列的最小出現(xiàn)頻次，默認為8，此處設(shè)置為10，推薦最小值為總數(shù)據(jù)量的百萬分之一(如1億條序列至少需要過濾掉頻次100以下的序列噪音)，可實現(xiàn)去除低豐度或測序噪音并極大地增加計算速度。-sizeout輸出頻次, -relabel設(shè)置輸出序列前綴(示例輸出序列ID：Uni1;size=17963)。vsearch -derep_fulllength temp/filtered.fa -output temp/uniques.fa -relabel Uni -minuniquesize 10 -si

33、zeout按97%聚類生成OTUs和去噪挑選ASVs是目前挑選特征序列的兩種主流方法。通常兩種分析方法的結(jié)果整體上比較類似，細節(jié)上會略有不同。下面按方法1或2分別介紹，用戶可根據(jù)實際需求選擇。方法1. 使用UPARSE ADDIN EN.CITE Edgar2013173(Edgar 2013)17317317Edgar, Robert C.UPARSE: highly accurate OTU sequences from microbial amplicon readsNature MethodsNature MethodsNat. MethodsNat Methods996-998101

34、02013/08/21AlgorithmsDatabases, GeneticHumansMetagenomicsMicrobiota/*genetics*PhylogenyRNA, Ribosomal, 16S/*geneticsResearch DesignSensitivity and SpecificitySoftware201308/18/onlineNature Publishing Group, a division of Macmillan Publishers Limited. All Rights Reserved.1548-709123955772Brief Commun

35、ication/10.1038/nmeth.2604/articles/nmeth.2604#supplementary-information10.1038/nmeth.2604(Edgar 2013)算法按97%的序列相似度聚類OTU。usearch -cluster_otus temp/uniques.fa -otus temp/otus.fa -relabel OTU_此方法累計使用次數(shù)最多、分析速度快，適合大數(shù)據(jù)或ASV方法結(jié)果規(guī)律不明顯時嘗試。方法2. 使用UNOISE算法 ADDIN EN.CITE ADDIN EN.CITE.DATA (Edgar and Flyvbjerg

36、2015)去噪生成ASV 。此方法是當前的主流，推薦優(yōu)先使用。類似于按100%的序列相似度聚類，或不聚類的方法。采用更先進的方法來鑒定測序過程中可能的錯誤，因此也需要消耗更多的計算時間。UNOISE算法雖然慢于UPARSE算法，但也比同類去噪算法Deblur和DADA2分別快10倍和100倍 ADDIN EN.CITE Amir201717(Amir等， 2017)1717017Amir, AmnonMcDonald, DanielNavas-Molina, Jose A.Kopylova, EvgueniaMorton, James T.Zech Xu, ZhenjiangKightley,

37、 Eric P.Thompson, Luke R.Hyde, Embriette R.Gonzalez, AntonioKnight, RobGilbert, Jack A.Deblur rapidly resolves single-nucleotide community sequence patternsmSystemsmSystemse00191-16222017/content/2/2/e00191-16.abstract/pmc/articles/PMC5340863/pdf/mSystems.00191-16.pdf10.1128/mSystems.00191-16(Amir等，

38、 2017)。此處-unoise3去噪結(jié)果默認前綴為Zotu，我們修改為主流使用的ASV。usearch -unoise3 temp/uniques.fa -zotus temp/zotus.fased s/Zotu/ASV_/g temp/zotus.fa temp/otus.fa(可選) 基于參考去嵌合。全頭(de novo)去嵌合時要求親本豐度為嵌合體16倍以上防止造成假陰性，而參考數(shù)據(jù)庫無豐度信息，只需1條序列在參考數(shù)據(jù)中沒有相似序列且與23條序列相似即判定為嵌體，因此容易引起假陰性(真實序列被當作假序列丟棄)，不推薦使用。如果必須要使用，由于已知序列不會被去除，選擇越完整的數(shù)據(jù)庫可降

39、低假陰性率。方法1. VSEARCH結(jié)合RDP數(shù)據(jù)庫去嵌合(快但容易假陰性)，推薦SILVA去嵌合(silva_16s_v123.fa)，但計算極耗時 (本例用時3小時，是RDP計算時間的30倍以上)。vsearch -uchime_ref temp/otus.fa -db $ea/usearch/rdp_16s_v16_sp.fa -nonchimeras result/raw/otus.faWindows系統(tǒng)下vsearch結(jié)果會添加了windows換行符M必需刪除，否則會出現(xiàn)換行混亂的問題。Mac/Linux系統(tǒng)無須執(zhí)行此命令。sed -i s/r/g result/raw/otus.f

40、a方法2. 不去嵌合。cp -f temp/otus.fa result/raw/otus.fa特征表生成和篩選使用vsearch的-usearch_global命令比對擴增子序列(temp/filtered.fa)至特征序列(result/raw/otus.fa)即可生成特征表，-threads設(shè)置整數(shù)使用計算機可用的多線程加速計算。vsearch -usearch_global temp/filtered.fa -db result/raw/otus.fa -otutabout result/raw/otutab.txt -id 0.97 -threads 4sed -i s/r/ res

41、ult/raw/otutab.txt基于物種注釋，(可選)可以篩選質(zhì)體和非細菌和非古菌去除并統(tǒng)計比例。使用USEARCH的sintax算法進行物種注釋。選擇RDP物種注釋(rdp_16s_v16_sp)數(shù)據(jù)庫具有體積小、分類速度極快(本例耗時15秒)的特點，但缺少真核數(shù)據(jù)無法注釋真核污染物來源詳細。SILVA數(shù)據(jù)庫(silva_16s_v123.fa)可以更好地注釋真核質(zhì)體序列來源，但速度較慢(本例耗時約3小時)。-sintax_cutoff設(shè)置分類結(jié)果的可信度閾值，范圍01之間，文章中常用0.6/0.8，取值越大注釋結(jié)果越可靠同時注釋比例也越低。注意，結(jié)果中第三列方向正常應(yīng)全為正向(+)，如

42、果全為反向(-)，請參考常見問題5中方法將第3步結(jié)果序列取反向互補。vsearch -sintax result/raw/otus.fa -db $ea/usearch/rdp_16s_v16_sp.fa -tabbedout result/raw/otus.sintax -sintax_cutoff 0.6為去除16S rDNA測序中的非特異性擴增和質(zhì)體污染，我們編寫了R腳本otutab_filter_nonBac.R實現(xiàn)選擇細菌和古菌(原核生物)、以及去除葉綠體和線粒體并統(tǒng)計比例，輸出篩選后并按豐度排序的特征表。輸入文件為原始特征表(result/raw/otutab.txt)和物種注釋(

43、result/raw/otus.sintax)，輸出篩選并排序的特征表(result/otutab.txt)、統(tǒng)計污染比例文件(result/raw/otutab_nonBac.txt)和過濾細節(jié)(otus.sintax.discard)，特征表的格式詳見文件或附表2。注：真菌ITS數(shù)據(jù)，請改用otutab_filter_nonFungi.R腳本，只篩選注釋為真菌的序列。查看腳本幫助，請運行Rscript $ea/script/otutab_filter_nonBac.R -h。Rscript $ea/script/otutab_filter_nonBac.R -input result/ra

44、w/otutab.txt -taxonomy result/raw/otus.sintax -output result/otutab.txt -stat result/raw/otutab_nonBac.stat -discard result/raw/otus.sintax.discard特征表篩選后，對應(yīng)的代表序列(otus.fa或附表3)、物種注釋信息也需要對應(yīng)進行篩選。cut -f 1 result/otutab.txt | tail -n+2 result/otutab.idusearch -fastx_getseqs result/raw/otus.fa -labels resu

45、lt/otutab.id -fastaout result/otus.faawk NR=FNRa$1=$0NRFNRprint a$1 result/raw/otus.sintax result/otutab.id result/otus.sintaxsed -i s/t$/td:Unassigned/ result/otus.sintax此外，如果上述篩選方案不適合你的研究，如去除的Chloroplast為你研究的對象，可以跳過此步不進行篩選，運行cp result/raw/otu* result/。接下來對最終的特征表進行統(tǒng)計，結(jié)果有助于優(yōu)化前面的分析方案，以及選擇下游分析合適的參數(shù)。us

46、earch -otutab_stats result/otutab.txt -output result/otutab.statcat result/otutab.stat圖2. USEARCH的特征表統(tǒng)計結(jié)果示例統(tǒng)計結(jié)果顯示了總讀長數(shù)量、樣本量、特征數(shù)量，可以了解特征表的總體數(shù)據(jù)量和維度信息；接下來是特征表中單元格數(shù)量、為0、1和10的數(shù)量和比例，了解特征中頻次分布情況；然后是特征在100%、90%和50%樣品中發(fā)現(xiàn)的數(shù)量，展示了特征在樣本中的流行頻率；最后是樣本測序量的分位數(shù)，對選擇合適的重采樣閾值非常有幫助。我們根據(jù)特征表統(tǒng)計結(jié)果，將選擇合適的參數(shù)對特征標進行等量重抽樣方式的標準化，以減

47、小由于樣本測序量不一致引起的多樣性差異，可實現(xiàn)更準確地多樣性分析。重采樣使用otutab_rare.R腳本調(diào)用vegan包 ADDIN EN.CITE Oksanen2007288(Oksanen等， 2007)28828817Oksanen, JariKindt, RoelandLegendre, PierreOHara, BobStevens, M Henry HOksanen, Maintainer JariSuggests, MASSThe vegan packageCommunity ecology packageCommunity ecology packageCommun. Ec

48、ol. Pack.Commun Ecol Pack631-637102007(Oksanen等， 2007)實現(xiàn)，并計算了6種常用alpha多樣性(richness、chao1、ACE、shannon、simpson和invsimpson)指數(shù)(vegan.txt或表2)。重采樣深度(-depth)參考特征表統(tǒng)計結(jié)果(圖2)選擇，一般默認按最小值重采樣。提高采樣深度可以保留樣本中更大的測序量，但也會剔除低于閾值的樣本。因此如果樣本測序量波動極大，盡量選擇合適的閾值重采樣以最大化保留測序量。mkdir -p result/alphaRscript $ea/script/otutab_rare.R

49、 -input result/otutab.txt -depth 32086 -seed 1 -normalize result/otutab_rare.txt -output result/alpha/vegan.txtusearch -otutab_stats result/otutab_rare.txt -output result/otutab_rare.statcat result/otutab_rare.stat結(jié)果顯示所有樣本重采樣后讀長數(shù)量均為32086。這樣特征表可以最大化減少測序量的影響，以便更準確評估多樣性。表2. Alpha多樣性指數(shù)示例SampleIDrichness

50、chao1ACEshannonsimpsoninvsimpsonKO123502692.0082686.8696.1328350.990308103.1788KO223162664.352661.866.174060.991875123.0733KO319352278.2522283.3435.8284520.98958295.98662由otutab_rare.R調(diào)用vegan包計算的6種常用alpha多樣性指數(shù)，圖中僅展示結(jié)果前4行作為示例。Alpha多樣性計算前面在特征表重采樣標準化時，計算了6種常用alpha多樣性指數(shù)。此外，USEARCH的-alpha_div命令可以快速計算18種a

51、lpha多樣性指數(shù)(alpha.txt)，各種指數(shù)的計算公式和描述詳見： HYPERLINK /usearch/manual/alpha_metrics.html /usearch/manual/alpha_metrics.html。這些結(jié)果我們常用于結(jié)合樣品元數(shù)據(jù)開展組間比較，或箱線圖展示組間異同。usearch -alpha_div result/otutab_rare.txt -output result/alpha/alpha.txt由于測序數(shù)據(jù)深度對多樣性影響較大，有時我們也關(guān)注不同測序樣量多樣性的變化，即可以判斷組間差異是否在不同測序深度下穩(wěn)定存在，同時確定測序量是否飽和并反映出結(jié)

52、果較真實的多樣性。USEARCH的-alpha_div_rare命令實現(xiàn)快速無放回百分數(shù)重采樣計算各樣本的豐富度(richness/observed_feature，詳見alpha_rare.txt或附表4)。結(jié)果可進一步可視化為樣本稀釋曲線，或分組帶誤差棒的稀釋曲線或箱線圖。usearch -alpha_div_rare result/otutab_rare.txt -output result/alpha/alpha_rare.txt -method without_replacementAlpha多樣性豐富度指數(shù)相似只代表物種數(shù)量相近，然而其中的物種種類可能完全不同。我們需要制作記錄每個

53、組指大于定豐度的物種是否存在的數(shù)據(jù)格式(表3和附表5)，用于組間比較物種共有和特有的情況。可以使用ImageGP在線( HYPERLINK /ImageGP/ /ImageGP/)選擇維恩圖(Venn diagram)、集合圖(Upset view)或桑基圖(Sankey diagram)等方式展示。表3. 用于比較各組特征共有/特有的數(shù)據(jù)示例特征ID分組IDASV_1AllASV_1KOASV_1OEASV_1WTASV_2KO我們通常結(jié)合元數(shù)據(jù)計算各組的豐度均值，如以Group列為分組信息計算原始計數(shù)的相對豐度并求組均值。Rscript $ea/script/otu_mean.R -inp

54、ut result/otutab.txt -design metadata.txt -group Group -thre 0 -output result/otutab_mean.txt因為特征的數(shù)量較大，而且低豐度的特征是否存在偶然性較大，準確性不高且與測序噪音無法區(qū)分。因此篩選大于某一豐度閾值結(jié)果，可實現(xiàn)數(shù)據(jù)降維并保留數(shù)據(jù)的主體，然后用于組間比較共有和特有的情況。如以平均豐度0.1%為閾值，可選0.5%或0.05%，得到每個組中符合條件的特征(表3)。awk BEGINOFS=FS=tif(FNR=1) for(i=2;i=NF;i+) ai=$i; else for(i=2;i0.1)

55、print $1, ai; result/otutab_mean.txt result/alpha/otu_group_exist.txtBeta多樣性計算Beta多樣性是群落整體結(jié)構(gòu)的降維分析方法，需要基于特征表計算樣本間的各種距離矩陣。常用的Unifrac算法 ADDIN EN.CITE Lozupone2010309(Lozupone等， 2010)30930917Lozupone, CatherineLladser, Manuel E.Knights, DanStombaugh, JesseKnight, RobUniFrac: an effective distance metric

56、 for microbial community comparisonThe Isme JournalThe ISME JournalISME J.ISME J1695201009/09/onlineInternational Society for Microbial EcologyCommentary/10.1038/ismej.2010.133/articles/ismej2010133#supplementary-information10.1038/ismej.2010.133(Lozupone等， 2010)考慮物種間的進化距離，這里我們使用usearch的-cluster_agg

57、命令基于代表序列獲得進化樹，然后再使用-beta_div基于特征表和進化樹計算多種矩陣矩陣，包括bray_curtis, euclidean, jaccard, manhatten, unifrac等，每類矩陣還分為有或無(binary)權(quán)重兩種(附表6展示Bray-Curtis距離矩陣示例)。mkdir -p result/beta/usearch -cluster_agg result/otus.fa -treeout result/otus.treeusearch -beta_div result/otutab_rare.txt -tree result/otus.tree -filen

58、ame_prefix result/beta/物種注釋分類匯總前在特征表篩選時已經(jīng)對特征序列完成了物種注釋，并根據(jù)注釋進行篩選。物種注釋存在命名混亂、分類級不完整和名稱缺失等問題。我們先對格式進行調(diào)整方便開展分析。調(diào)整物種注釋為特征ID和7級分類注釋的兩列格式(表4和附表7)。注意7級分類可能存在不完整的情況，可能是該特征沒有相近種的報導，也可能是參考物種注釋自身不完善。cut -f 1,4 result/otus.sintax |sed s/td/tk/;s/:/_/g;s/,/;/g;s/g;s/Chloroplast/ result/taxonomy2.txt表4. 物種注釋2列格式示例

59、特征ID分組IDASV_1k_Bacteria;p_Actinobacteria;c_Actinobacteria;o_Actinomycetales;f_ThermomonosporaceaeASV_2k_Bacteria;p_Proteobacteria;c_Betaproteobacteria;o_Burkholderiales;f_Comamonadaceae;g_Pelomonas;s_Pelomonas_puraquaeASV_3k_Bacteria;p_Proteobacteria;c_Betaproteobacteria;o_Burkholderiales;f_Comamona

60、daceae;g_Pelomonas;s_Pelomonas_puraquae物種注釋的7級分類，也提供了特征表多角度降維分析的可能，可以把特征表按物種注釋信息合并為門、綱、目、科、屬級別的表，以進行更容易與現(xiàn)在科學發(fā)現(xiàn)結(jié)合的討論。首先將物種注釋轉(zhuǎn)化為7級分類的8列表格(表5和附表8)，其中缺失的分類級別填充為未分類(Unassigned)。awk BEGINOFS=FS=tdelete a; ak=Unassigned;ap=Unassigned;ac=Unassigned;ao=Unassigned;af=Unassigned;ag=Unassigned;as=Unassigned; sp