美格基因宏基因组结题报告

概 3.1数据预处 3.1.2数据质 3.1.3数据质量统 ORF..................................................................................................................................................... 3.3.1序列聚 3.3.2丰度计 3.3.3数目Venn KEGG注 eggNOG注 CAZy注 物种-功能-丰度聚类Heatmap 组间Wilcox秩和检 两组间T检 附 6.1................................................................................................................................................................ 宏组（tagnome，也称微生物环境组或元组，是由ndlmn等1998tomsofthetotlmirobiotafoundinntur物的，目前主要指环境样品中细菌、古菌和真菌的组总（metgnomi就是一种以环境样品中的新的微生物研究方法。最开始的宏组学技术主要被用来从环境微生物中发掘新生物分子，包括基于功能性筛选（基于代谢功metbolicfution（mrkrgn）A的克隆以及小片段或大片段载体的构建，构建好的载体需要转入适当的宿主菌中，一般为Escherichiacoli。的研究利用宏组文库构建加上短序列（shotgunsequencing）方法研究了酸性矿山废水生物膜上及Sargasso海洋中微生物群落的组成及功能构成。然司的SOLID仪、Illumina的MiSeq、HiSeq仪的应用，带来的通量幅度的提高，例如SorcererII全球海洋采样计划（TheSorcererIIGlobalOceanSamplingproject，GOS）的开展。另一方面，宏组避免了由体外扩增系genes（bias组的重建成为可能，同时基于生境中主要代谢网络的构建对于微生物实验分析流程实验流DNADNA样行DNA重新提取或对现有DNA进行纯化等。检测合格的DNA样品加入fragmentationbuffer，采用超声破碎仪进行随机打断，将打断后得到的短片段DNA用于文库构建，建好的文库需进行文库质检，对于质检合格的文库将采用Reads以FASTQ（简称为fq）文件格式，其中包含序列（Reads）的序列信息1生物信息分析流程质在获得每个样品宏组数据之后，首先需要对的数据质量进行流程采用为Trimmomactic，参数为(LEADING:3TRAILING:3SLIDINGWINDOW:5:20MINLEN:50)。序列组装及原始的序列剪切成更短的片段K-mers来组合拼接，也更适合二代结果，比如MEGAHIT，SOAPdenovo，Velet，IDBA，TruSPAdes等用的就是采用这种方法。在此我们组装上利用 /megahit）CleanDatadenovo拼接（K-merk-min35，k-max95，k-step20。并将各样品未被利用上的Reads放在一起进行混合组装，以期发现样品中的低丰度物种信息。拼接完成后，组装得到的Scaffolds从N连接的Scaftigs进行后续分析。从单样品和混合组装后的Scaftigs出发，采用MetaGeneMark来开放阅读框（ORF，评估出来的ORFs数目、长度非冗余集构建及丰度计将各样品和混合组装的ORFs结果放在一起，采用Linclust -e0.001min-seq-id0.9c0.80）Cluster取最长的序列为代表性序列到非冗genecatalogu（Unigenesbbmap将质控后CleanData比对至genecatalogue，来计算得到各Unigene在各样品中的丰度信genecatalogue各在各个样品中的有无，进行core-pan分析、组间数目差异分析、数目Venn图分析；基于genecatalogue中各在各样品中的物种及功能注释将非冗余集的Unigenes序列与NCBI-NR数据库比对进行物种注释，获得Unigenes的物种注释信息，并结合丰度表，获得各个分类层级的物种组成和丰度信息。另外，利用metaphlan2将CleanReads快速分类到物种KEGGeggNOGCAZym物种及功能组成分基于物种丰度表和功能丰度，可以绘制样品间共有特有物种（）Venn样本比较分析Beta多样性分析来比较不同样品在物种（功能）组成方面存在的差异样品聚类热图及样PCAPCoA图；当有分组重时利Anosim，Adonis等组间群落差异显著性检验，揭示组间物种（功能）的差异是否显著。物种及功能差异分能分类层级上每个物种（功能）在组间的差异情况，寻找具有统计学差异的环境因子关联分析CCA/RDA、db-RNA、VPA、Man检验、偏Man检验等分析，找到与样品关联与模型分数据预处理原始序列数据Reads的序列信息及其对应的质量信息。PE文库的数据结果中，每个样品均有两端的Reads，并且Reads的顺序是严格一致、相互对应的。@HWI-D00524:166:C6AJYANXX:4:1101:1710:19731:N:0:CACTCA@HWI-D00524:166:C6AJYANXX:4:1101:1710:19731:N:0:CACTCA列名称由Illumina仪产生；第2行是read的碱基序列，第4行是read中每个碱基量系统的Phredqualityscore（33或64），即为该碱基的质量值（Phred值。raw_data_list.txt: 数据质控其是末端（3’端）呈现明显的下降的趋势，因此RawReads会存在一定比例的低质量数据，同时也存在建库过程中产生的PCR重复。为了保证后续分析的准确可靠，质控流程采用软件为Trimmomactic，参数为(LEADING:3TRAILING:3SLIDINGWINDOW:5:20碱基质量小于20，则切除；如果样品来源于宿主（比如人或动物的粪便通过将reads比对宿主组，去除比对到组的reads来进行宿主过滤，推荐使用BWA（ Read1Read2CleanData，使用FastQC（）对 的clean_data_list.txt: 数据质量统计表1数据质量统计表（展示部分，详细资料请序列组装及序列组装利用MEGAHIT（ k-k-20，由于Scaffolds里可能存在嵌合体序列，我们严格的将Scaffolds从一个或者多个连续的N处reads500bpScaftigsN50等指标评估组N50：指Scaftigs长度覆盖50%所有核苷酸的最大序列群长度，把Scaftigs按长度从大到小排序，ScaftigsScaftigs的长度。N90含义与N50相似。相比序列平均长度，N50（N90）更能准确表示此次的序列拼接效果；3Scaftigs示对应长度的Scaftigs数占总Scaftigs数的百分比。每条所有ORFScaftigs（500bp） 90bpORFbbmapcleanreadsORFORFcleanreads0cleanreads长度≥ORFORFORF表3有效ORF信息统计（展示部分，详细资料请predict_gene.summary.xls）SampleGeneNumberTotallen.(bp)Averagelen.(bp)GCpercent(%)W1S1L 注：x轴表示长度，柱状图对应左侧坐标，表示对应长度的数量；曲线对应右侧坐标，表示对应/*/protein.faa:CDS子表翻译得到的氨基酸序列*/gene.gff:ORFs征的描述文件(gff)tab9非编码RNA(ncRNA)是指一类本身不携带可以翻译为蛋白质的遗传信息的，可以执行多种生物学功能的RNA分子，主要包括有rRNA、tRNA和microRNA。rRNA的使用Barrnap，tmRNA的则利用Aragorn。 /用的reads混合拼接所得的序列 结非冗余集构建及丰度计序列聚类将各样品及混合组装的结果放在一起，采用Linclust /articles/s41467-018-04964-5)进行聚类（默认参数为-e0.001--min-seq-id0.9-c0.80）及去冗余，将每个Cluster中最长的序列作为代表性序列，以获得非冗余的genecatalogue（Unigenes。表5聚类结果文件(展示部分，详细资料请Unigene ORF 表6非冗余集的代表序列及其长度统计表(展示部分，详细资料请Unigene 表7非冗余集数目和长度统计表(展示部分，详细资料请 Gene GC Unigenes:非冗余集GeneNumber:非冗余集的代表序列数Averagelen.:所有非冗余集序列的平均长度 丰度计算

=𝑟𝑘

1 𝑖=1注:rk指比对到k的reads数目;Lk指k的长……的代表序列 Venn丰度样品间相关性分析要求较高。根据丰度结果，我们使用R语言计算出所有样本两两之间的相关系数6 物种与功能注释物种注释使用BLASTP（Version2.2.31+，）将非冗余的UnigenesNCBI-NRBLAST比对（e-value≤0.0001不同的分类级别，为了保证其生物意义，采取LCA算法（应用于MEGAN的系统分Gene……对丰度;taxonomy为代表此Gene被注释的物种分类Unigenes_gene_nr.txt:UnigeneNRKronaKrona分析结果来展示各样品中不同分类学水平物种的相对丰度，Krona示例图7使用Krona多详细的信息请参考Krona展示结果详解。KEGG，从分子水平信息，尤其是大型分子数据集生成的组和其他高通量实验技术的实用数的Pathway；第四层为每个代谢通路图的具体注释信息。我们统计各样品在第一和第二层次的KEGG注释结果，并提供注释上的所有代谢通路注释图。使用diamond将非冗余集序列与KEGG的数据库进行的同源性比对(e-value<=0.001)。由于每一条序列比对结果可能不止一条，为保证后续研究的表10的相对丰度及KEGG数据库注释表（展示部分，详细资料请 W1S1L… Functio

Enzy

nunigene_1

7385297

map00910(nitrate/nitrite K1transpo map02010(6.3.-00438system binding)Function：KO对应的功能描述及编码的酶的Hyperlink：KO的网页8命名规则：K+num（ID号，表示在所有同源物种中具有相似结构或功能的一类同源蛋白图9KEGGLevel1水平数目柱状图KOko_sum.xls:koeggNOGeggNOG（evolutionarygenealogyofgenes:Non-supervisedOrthologousGroups，）数据库是利用Smith-Waterman比对算法对构建的直系同源分类学功能注释。目前该数据库（v4.5.1）1902031个物种value≤0.001 …

…bactNOG03029|bproNOG00822|COG

carbonatetransport

Class：OG所属的功能大类注：条形图对应的纵坐标指注释到的数目，横坐标代表eggNOG各功能类型的代码，简写全称见右eggNOG eggNOGOrthologousCAZy GHsPLs水化合物结合模块（Carbohydrate-BindingModules,CBMs。使用CAZy数据库的对应工具hmmscan将非冗余集与CAZyD（ersion:2016.7.15）进行比对，比对参数设置期望值e-value≤1e-5，获得对应的碳水化合物活性酶注释信息。我们统计各样品在功能大类层级的CAZy注释结果。表12相对丰度及CAZy数据库注释表(展示部分，详细资料请 name cazyclass 注：Gene_ID：unigenes的代表序列ID；Gene_ID与name中间的各列是在各个样品中的丰度值；name：序列注释到的CAZy；ncbi-cdd:pfam资料库单蛋白质;cazyclass：CAZy注：条形图的纵坐标指匹配的数目，横坐标代表CAZy各功能类型的代码，简写全称见右侧的图CAZyARDB- 抗生素抗性注ARDB(AntibioticGenesDatabase，)又称耐药数素抗性类型所抑制的抗生素名称(Profile)等信息。BacMet(AntibacterialBiocide&MetalGenesDatabase)是一个抗菌剂和金属抗性数据库，其收录有753个经实验验证的抗性和155,512个的抗性(截止于11March2018)。信息，及其所抑制的抗生素名称，结合unigenes代表序列的相对丰度获得相对丰度及BacMet数据库注释表。表13ARBD数据库注释表(展示部分，详细资料请Require AminoglycosideN-

gentamicin,sisomicin-描述；Profile：抗生素抗性类型所抑制的抗生素名称;RequireUnigenes:。图12 Profile数目条形图 注释Unigenes_abundance_bacmetFunc.tsv:BacMet注释表Unigenes.ardb.output.txt:ARDB注释结果Unigenes.ardb.bltOut.txt:ARDB注释结果 VFDB致病菌毒力因子注释VFDB数据库（htt /VFs/main.htm）全称为VirulenceFactorsofPathogenicBacteria，是专门用于研究致病细菌、衣原体和支原体致病因子的数据库。其中，1796，30178表14相对丰度及VFDB数据库注释表(展示部分，详细资料请

注：Gene_ID：unigenes的代表序列ID；Gene_ID与VFDBID中间的各列是在各个样品中的丰度值；VFDBID：毒力因子的；Virulencefactor：VFDB数据库中的毒力信息描述；Function：VFDB数据库中的功能描述；Sourcespecies：毒力因子的物种来源。13每个样品VFDB3.4-3.4-物种与功能注释/3.4.6-VFDB致病菌毒力因子注释/Unigenes.VFs.tsv:VF注释结果Unigenes.vf_sum.txt:VFDB二级分类统计表Unigenes.vf_pieplot.png:VFDB二级分类丰度饼图PHI病原与宿主互作注释PHI（PathogenHostInteractions，）病原与宿主互作数据原与宿主互作关系，264种病原体，176种宿主等信息。对丰度获得PHI数据库功能注释信息。3.4-3.4-物种与功能注释/3.4.7-PHI病原与宿主互作注释/Unigenes.phi.annotation.tsv:PHI注释结果TCDB转运蛋白分类注释谢产物以及信号转导等广泛的细胞活动中起着重要的作用。转运蛋白类别数据库多种非冗余转运系统，分为1000多个转运蛋白。该分类系统采用类似于酶分类的ECfamilies、subfamilies和transportsystems，每个层级对应TC里的一个字母或数字，每表15相对丰度及TCDB转运蛋白资料库注释表(展示部分，详细资料请 …n

The

binding unigene_1738 Q5513.A.1.16 (ABC)

GO:00151 注：Gene_ID：unigenes的代表序列ID；Gene_ID与Proteinid之间各列为在各样品中的丰度；Proteinid及TCDBid：转运蛋白在TCDB数据库中的对应的系统分类；IDDescription：转运蛋白描述信DataBank中对应PDBID；Pfam：PfamPfamID；GO：GO数据库中对GOID。3.4-3.4-物种与功能注释/3.4.8-TCDB转运蛋白分类注释/Unigenes.tcdb.annotation.txt:TCDB注释结果物种与功能组成分析物种-功能-丰度柱状图 phylum/*.pdf:门阶层class/*.pdf:纲阶层order/*.pdf:目阶层family/*.pdf:科阶层;genus/*.pdf:属阶层;speices/*.pdf:种阶层物种-功能-VennCAZy的注释结果。图1512物种数目图15.2功能数目的Unigene个数。 /Venn/:物种-功能-丰度聚类Heatmap对于物种丰度Heatmap图，基于各样品(组)物种分类的结果及其相对丰度，默认在phylumspeciesHeatmapKEGG的KOeggNOGOG图16.1物种相对丰度heatmap图（门 图16.2物种heatmap图（属图16.3功能heatmap图(KEGG 图16.4功能heatmap图(eggnog某个分类上的Z值为样品在该分类上的相对丰度和所有样品在该分类的平均相对丰度的差除以所有样品Heatmap图phylum/*.pdf:门阶层class/*.pdf:纲阶层order/*.pdf:目阶层family/*.pdf:科阶层;genus/*.pdf:属阶层;speices/*.pdf:种阶层Heatmap图/KEGG/Heatmap图/NOG/Heatmap图物种-功能-Bubblephylumspecies水平、KEGG、eggNOG及CAZy各分类结果的丰度分布，以帮助寻找优势及差异的物种、KEGG、NOG及CAZy功能等信息，以及分析其变化趋势。图17.1物种Bubble图(屬 图17.2功能Bubble图(KEGGLevel3.5-3.5-物种与功能组成分析/3.5.4-物种-功能-Bubble图/-物种Bubble图/top50*_[L2/L3/L4/L5/L6/L7].pdfL2~L7分别以门、纲、目、科、属、种丰度所绘制top50_*_[ko/L1/L2].pdf:各分类阶层绘制top50_*_[L1/L2].pdf:各分类阶层绘制:Ternary三元相5的物种-功能-ternaryplot（三元相图-功能-在不同分类水18三元相图(门phylum/*.pdf:门阶层class/*.pdf:纲阶层order/*.pdf:目阶层family/*.pdf:科阶层;genus/*.pdf:属阶层;speices/*.pdf:种阶层注:*-VS*-VS-*为任选三样本(组)比较分析，可根据您的需求找寻想要比较的组合Single为单样本,组名样本比较分析分别基于物种、KEGGKO、eggNOGOG、CAZY的第二层级功能类（Level2）等相对丰度数据，利用beta多样性距离矩阵等统计方法，分析样品间的差异情况。样本间相关性分析重复间相关系数要求较高。根据丰度结果，我们使用R语言计算出所有样本两两之间的相关系数（Sperman相关系数spearmam_*_[ko/L1/L2].pdf:各分类阶层绘制spearmam_*_[L1/L2].pdf:各分类阶层绘制spearmam_*_[L1/family].pdf:各分类阶层绘制 聚类与柱状图组合分Bray-Curtis距离矩阵以及层次聚pairgroupmethodwitharithmeticmean）算法构建树状结构，得到树状关系形式用于可视化注：左侧是UPGMA聚类树结构，每个分支代表一个样本，右侧是各样本在种水平、KEGG(Level2)、eggNOG(Level2)OGCAZY(Level1)的物种-功能-相对丰度柱状图，不同颜色代表不同物种-功能-。图3.6-3.6-样本比较分析/3.6.2-UPGMA聚类与柱状图组合分析/-UPGMA物种聚类与柱状图:3.6-样本比较分析/3.6.2-UPGMA聚类与柱状图组合分析/-UPGMA功能聚类与柱状图:3.6-样本比较分析/3.6.2-UPGMA聚类与柱状图组合分析/-UPGMA功能聚类与柱状图:3.6-样本比较分析/3.6.2-UPGMA聚类与柱状图组合分析/-UPGMA功能聚类与柱状图:PCA主成分分PCA（PrincipalComponentysis，主成分分析）是一种传统的数据降维方法，常用注：一个点代表一个样本，点与点空间距离表示物种-功能-组成结构的差异程度。主成分1（PC1）和主成分2（PC2）是造成样品间差异的第一和第二大特征值，百分比表示两个主成分对样品差异的方差贡献3.6-3.6-样本比较分析/3.6.3-PCA主成分分析/-物种PCA主成分分析/[组名]PCA_*_[L2/L3/L4/L5/L6/L7].pdfL2~L7分别以门、纲、目、科、属、种丰度所绘制:PCA_*_[L1/L2].pdf:各分类阶层绘制:注:各命名中如出现(1)ellipase表加椭圆区;(2)label表加样本名称PCoA主坐标分PCA是基于原始的物种-功能-两种距离来进行PCoA分析，并选取贡献率最大的两个主成分进行作图展示。注：一个点代表一个样本，点与点空间距离表示物种-功能-组成结构的差异程度。主成分1（PC1）和主成分2（PC2）是造成样品间差异的第一和第二大特征值，百分比表示两个主成分对样品差异的方差贡献物种:功能:功能:功能:注:各命名中如出现(1)ellipase表加椭圆区;(2)label表加样本名称NMDS非度量尺度分NMDS（NonmetricMultidimensionalScaling，非度量尺度法）常用于比较样本组（PCA、PCoA）的缺点，能更好地反映生态学数据的非线性结构。样品的空间距离物种:功能:功能:功能:注:各命名中如出现(1)ellipase表加椭圆区;(2)label表加样本名称MDS标度法分多元数据分析技术,是一种将空间的研究对象(样本或变量)简化到低进行定位、物种:功能:功能:功能:注:各命名中如出现(1)ellipase表加椭圆区;(2)label表加样本名称组间比较分析PCA、PCoA等分析之后，虽然肉眼能大概分辨分组是否可以清组间群落差异显著性检验又称DissimilarityTest（非相似性检验（Anosim析（MRPP）这三种方法，并默认基于物种-功能-Bray-Curtis、Jaccard这三种距离矩阵分别进行分析。这些方法不单可以输出组间群落结构差异程度（R值能得到检验的显著性（P值）等结果。 置换多元方差分析Adonis（PermutationalMmultivariateysisofVariance，PERMANOVA）又称置换多得到R和P值。表18Adonis分析结果(展示部分，详细资料请

1vssub-

l

)

9

注：表中dis表所使用计算distancematrice方法[bray/jaccard];method表[CCA/DCA/RDA/NMDS/MDS/PCoAtest表用[adonis/anosim/mrppSub-condition1vssub-condition2表任两组子分组做比较;Df为度；SumsOfSqs为总方差，又称离差平方和；MeanSqs为均方（差），即SumsOfSqs/Df；F.Model为F检验值；R2为不同分组对样品差异的解释度，即分组方差与总方差的比值，R2值接近1表示分组对差异的解释度越高；PrP值，小于0.05说明检验的可性度高。物种:功能功能功能Anosim组间相似性分间的距离矩阵，根据分组信息分别计算组内及组间的差异大小，得到R值，使用置换检验对分组进行显著性分析得到P值。下表以基于物种的Bray-CurtisJaccard距离矩阵的Anosim分析结果为例进行展示。

Leaf.W1vs 注：表中dis表所使用计算distancematrice方法[bray/jaccard];method表[CCA/DCA/RDA/NMDS/MDS/PCoAtest表用[adonis/anosim/mrppSub-condition1vssub-condition2表任两组子分组做比较;R值理论上范围为-11，但一般介于01之间，R值越接1表示组间差异越大，组间R值越大说明采样或分组越合理，P值小于0.05时说明检验的度高。物种:功能功能功能MRPP多响应置换过程分析表20MRPP分析结果(展示部分，详细资料请sub-condition1vs

Leaf.W1vs 注：表中dis表所使用计算distancematrice方法[bray/jaccard];method表[CCA/DCA/RDA/NMDS/MDS/PCoAtest表用[adonis/anosim/mrppSub-condition1vssub-condition2表任两组子分组做比较A为不同分组对样品差异的解释度，A0说明组间差异大于组内差异，A值小于0说明组内差异大于组间差异，A值越大表示分组对差异的解释度越高；observeddelta值越大说明组内差异越大；expectdelta值越大说明组间差异大；Significance值小于0.05说明检验的可性度高。物种:功能功能功能物种与功能差异分析Speces组间Wilcox秩和检Wilcox秩和检验（Wilcoxonrank-sumtest，也叫曼-U检验（Mann–WhitneyUtest丰度的两组间差异分析。以下展示Wilcox秩和检验的分析结果： 注：文件名则显示在特定分组(Condition1)下，子分组(A)与子分组(B)之间的分析,L2表菌种阶层;ID:第一栏为各阶层物种-功能-名Wilcoxon_p.vlaueWilcoxp-value值。Wilcox_sig0代表无显着差异(p>-0.051表p<0.05且呈现含量下降(BA含量下降1表p<0.05且呈现含量上升(B较A含量上升)。--注：文件名则显示在特定分组(Condition1)下，子分组(A)与子分组(B)之间的分析L7表菌种阶层enriched:列出p<0.05且含量上升(B相较A)的物种-功能-名depleted:列出p<0.05且含量下降(B相较于A)的物种-无25组间差异显著物种-功能-组间WilcoxAll_pvalue_ttest_[组名]_[子组名]-[子组名]_[L2/L3/L4/L5/L6/L7].txt:列出检验下的p-value值heatmap_sig_ttest_[组名]_[子组名]-[子组名]_[L2/L3/L4/L5/L6/L7].txt:组间差异显著物种All_pvalue_ttest_[组名]_[子组名]-[子组名]_[ko/L1].txt:列出检验下的p-value值sig_ttest_[组名]_[子组名]-[子组名]_[ko/L1].txt:列出检验下p<0.05的功能heatmap_sig_ttest_[组名]_[子组名]-[子组名]_[ko/L1].txt:组间差异显著功能热图All_pvalue_ttest_[组名]_[子组名]-[子组名]_[L1/L2].txt:列出检验下的p-value值sig_ttest_[组名]_[子组名]-[子组名]_[L1/L2].txt:列出检验下p<0.05的功能heatmap_sig_ttest_[组名]_[子组名]-[子组名]_[L1/L2].txt:组间差异显著功能热图3.8-物种与功能差异分析/3.8.1-wilcoxon/-wilcoxon/CAZy/[组名]All_pvalue_ttest_[组名]_[子组名]-[子组名]_[L1/family].txt:p-value值sig_wilcoxon_[组名]_[子组名]-[子组名]_[L1/family].txt:列出检验下p<0.05的功能两组间T检test据满足正态分布，由于大部分实验设计属于独立样本，我们默认采用独立样本t检验。验进行功能丰度的两组间差异分析。以下展示t检验的分析结果： 注：文件名则显示在特定分组(Condition1)下，子分组(A)与子分组(B)之间的分析,L7表菌种阶层;ID:第一栏为各阶层物种-功能-名t-test_p.vlaue