温馨提示：请使用火狐或者Chrome的网页浏览器来查看报告

真菌框架图分析报告

---

一、 概述

真菌基因组研究，是通过基因组测序和组装获得真菌全基因组序列，并对其进行结构和功能研究的方法。 真菌基因组测序为真菌的研究提供强有力的支撑，可用于预测真菌的重要基因和蛋白以了解其功能和可能机制。

真菌框架图分析对二代测序数据进行组装，获得真菌全基因组序列，并使用基因组组分分析、 基因功能分析、比较基因组分析、群体进化分析等分析方法，对真菌基因组相关信息进行全面、详尽地了解。

二、 工作流程概述

2.1 实验流程

2.1.1 文库构建及库检

（1）DNA的提取

采用SDS或STE的方法对样本的基因组DNA 进行提取，之后利用琼脂糖凝胶电泳检测DNA 的纯度和完整性，利用Qubit进行定量。

（2）Pacbio平台建库及库检

采用SMRT bell TM Template kit（version 1.0）试剂盒构建20K SMRT Bell文库，将经电泳检测合格的DNA样品用Covaris g-TUBE打断成构建文库所需大小的目的片段，经DNA损伤修复及末端修复，使用DNA黏合酶将发卡型接头连接在DNA片段两端，并使用AMpure PB磁珠对DNA片段进行纯化，使用BluePipin片段筛选特定大小的片段，使用AMpure PB磁珠对SMRT Bell文库进行浓度筛选，随后修复DNA损伤，再次使用AMpure PB磁珠对SMRT Bell文库纯化，将构建好的文库经Qubit浓度定量，并利用Agilent 2100检测插入片段大小，最后用PacBio平台进行测序。

（3）Illumina平台建库及库检

经电泳检测合格的 DNA 样品用 Covaris 超声波破碎仪随机打断成长度约为350bp的片段。处理完成后的DNA片段，使用 NEBNext®Ultra™ DNA Library Prep Kitfor Illumina(NEB, USA)试剂盒，经末端修复、加A尾、加测序接头、纯化、PCR扩增等步骤完成整个文库制备。

文库构建完成后，先使用 Qubit 2.0 进行初步定量，稀释文库至 2ng/ul，随后使用 Agilent 2100 对文库的插入片段进行检测，insert size符合预期后，使用 Q-PCR 方法对文库的有效浓度进行准确定量，以保证文库质量。

2.1.2 建库类型

2.1.3 上机测序

库检合格后，把不同文库按照有效浓度及目标下机数据量进行PacBio Sequel和Illumina NovaSeq PE150测序。

2.2 生物信息分析流程

图1 生物信息分析流程图

信息分析分以下几个步骤：

1 原始下机数据处理：此步骤过滤测序质量值低的reads，保留高质量reads，过滤后的数据称为Clean Data；

2 样品组装：进行基因组组装，得到能反映样品基因组基本情况的序列文件，并对组装结果进行评价；

3 基因组组分分析：组装完成后，分析样品基因组的成分，包括编码基因、非编码RNA、重复序列等基因组成分的预测；

4 功能注释：针对编码基因序列进行不同数据库的功能注释，包括常用的KEGG、KOG数据库、针对致病性的数据库；

5 基因组比较分析：此步骤从基因组、基因两层面分别比较样品与参考基因组的差异，包括共线性分析、基因家族构建、共有/特有基因分析、物种进化树构建；

6 群体进化分析：此步骤包括共有基因及特有基因、基因家族分析和群体进化分析等内容；

三、 分析结果

3.1 数据概况

customer_files/01reads_info/
├── 01.01trim_pic
│   └── [样本].svg			[各样本二代测序数据reads数据过滤统计图]
├── 01.02ErrorRate
│   └── ErrorRate.[样本].svg			[各样本二代测序数据reads错误率分布图]
├── 01.03BaseDistribution
│   └── GCContentDistribution_1.[样本].svg		[各样本二代测序数据reads碱基含量分布图]
└── 01.04QualityDistribution
    └── QualityDistribution.[样本].svg		[各样本二代测序数据reads测序质量分布图]

3.1.1 二代测序数据概况

测序获得的原始数据中包含少量带有测序接头或测序质量较低的reads，为保证数据分 析的质量及可靠性，需要对原始数据进行过滤。本分析使用Trimmomatic[1]软件对 测序数据进行过滤，过滤前后各部分reads所占比例均在饼图中呈现。本项目测序数 据过滤情况统计图见结果文件 customer_files/01reads_info/01.01trim_pic/。下图为其中一个样本作为示例：

data.filter.sallet_1data.filter.sallet_2data.filter.sallet_3

注：
   Both Surviving：该部分为过滤后read1和read2均被保留的数据
   Forward Only Surviving：该部分为仅在read1中被保留的数据
   Reverse Only Surviving：该部分为仅在read2中被保留的数据
   Dropped：该部分为因接头或数据质量等原因被丢弃的数据

3.1.2 测序错误率分布

测序过程本身存在机器错误的可能性，测序错误率分布检查可以反映测序数据的质量， 序列信息中每个碱基的测序质量值保存在fastq文件中。如果测序错误率用e表示， Illumina的碱基质量值用Qphred表示，则有：Qphred=-10log10(e)。 Illumina Casava 1.8版本碱基识别与Phred分值之间的简明对应关系见下表。

表3.1 Illumina Casava 1.8版本碱基识别与Phred分值之间的简明对应关系

当前RNA-seq测序技术，测序错误率分布存在以下两个特征。

测序错误率随着测序序列(Sequenced Reads)长度的增加而升高。 这是由测序过程中化学试剂的消耗导致的，为Illumina高通量测序平台所具有的特征。

前6个碱基具有较高的测序错误率，此长度恰好为RNA-seq建库过程中反转录所需的随机引物长度。 前6个碱基测序错误率较高是因为随机引物和RNA模版的不完全结合。 此特征为illumina高通量测序平台的共有特征。

在该部分分析中，若样品80%的测序序列错误率在0.1%以下即为合格。 本项目测序数据的错误分布图见结果文件 customer_files/01reads_info/01.02ErrorRate/。 下图为其中一个样本作为示例：

ErrorRate.sallet_1ErrorRate.sallet_2ErrorRate.sallet_3

注：
pos：横坐标为reads碱基位置，其中从0-150为read1碱基位置，151-300为read2碱基位置。
errorRatio：纵坐标为碱基错误率。

3.1.3 GC含量分布

核苷酸序列中鸟嘌呤(G)和胞嘧啶(C)所占的比例称为GC含量。 GC含量在物种间存在一定特异性，但由于反转录过程中所使用的6bp随机引物， 会引起前几位碱基在核苷酸组成上有一定偏好性，产生正常波动，随后则趋于稳定。 对于NEB普通建库方法，由于序列的随机性打断和双链互补等原则， 理论上测序读段在每个位置的GC及AT含量应分别相等，且在整个测序过程基本稳定不变， 呈水平线。而对于链特异性建库而言，由于只保留了单链信息， 可能会出现AT分离或GC分离现象。

本项目各个样本的GC含量分布见结果文件 customer_files/01reads_info/01.03BaseDistribution。 下图为其中一个样本作为示例：

GCContentDistribution_1.sallet_1GCContentDistribution_1.sallet_2GCContentDistribution_1.sallet_3

注：
base position：横坐标为reads碱基位置；
base count：纵坐标为该碱基所占比例。

3.1.4测序数据质量分布

测序数据的质量主要分布在 Q30（≥80%）以上，这样才能保证后续分析的正常进行， 根据测序技术的特点，测序片段末端的碱基质量一般会比前端的低。

本项目各个样本测序数据质量分布见结果文件 customer_files/01reads_info/01.04QualityDistribution。 下图为其中一个样本作为示例：

QualityDistribution.sallet_1QualityDistribution.sallet_2QualityDistribution.sallet_3

注：
pos：横坐标为碱基位置；
quality：碱基的测序质量值。

3.2 基因组概况

customer_files/02assemble/
├── 02.01genome_survey/
│   ├── k_mer_table.txt		[各样本kmer信息统计表格]
│   └── [sample].linear_plot.png		[各样本GenomeScope图]
└── 02.02genome_assembly/
    ├── [sample].scaffold.fasta			[各样本scaffold序列]
    ├── [sample].contig.fasta			[各样本contig序列]
    ├── scaffold.table.txt		[scaffold信息统计]
    ├── contig.table.txt		[contig信息统计]
    └── [样本].GC_depth_count.svg		[各样本GC_depth图]

3.2.1 基因组大小估计

组装前采用基于K-mer统计的分析方法来估计基因组大小。用15-mer对测序数据进行分析，估算得到样品的基因组大小见表2，K-mer统计分布图如图2所示：

表2 K-mer分析数据统计

sampleID	K-mer	K-mer number	K-mer depth	Genome size	Revised gsize	Heterozygous Rate	Repeat Rate
sallet_1	15	43243309	45.73633319781333	7,061,672	7,063,929	0.188	0.00218
sallet_2	15	43593600	46.05903192211701	7,096,864	7,100,711	0.193	0.00255
sallet_3	15	47824612	50.495170436510804	7,094,176	7,102,095	0.199	0.00282

说明：
sample ID：样本名称；
K-mer：K-mer值；
K-mer number：从reads中提取用于统计的kmer个数；
K-mer depth：统计的kmer平均深度；
Genome size：估计的基因组大小；
Revised gsize：纠正后估计的基因组大小；
Heterozygous Rate: 估计的基因组杂合比例；
Repeat Rate：估计的基因组重复比例。

sallet_1.linear_plotsallet_2.linear_plotsallet_3.linear_plot

3.2.2 组装分析

从各样品质控后的 Clean Data 出发，进行基因组组装，得到能反映样品基因组基本情况的序列文件，并对组装结果进行评价。

基因组组装的具体处理步骤如下：

1）经过预处理后得到 Clean Data，使用 SOAP denovo组装软件进行组装；

2）首先选取不同的 K-mer（默认选取35、47、59、71、83、95、107、119）进行组装，根据项目类型选择最优的kmer（细菌项目选取最少scaffold，真菌项目选取最大的N50）；利用最优kmer并调节其他参数（-d -u -R -F等）再次筛选得到初步组装结果；

3）采用krskgf、gapclose等软件对初步组装结果进行优化和补洞，从而得到最终的组装结果；

4）过滤掉 500bp以下的片段，并进行评估和统计分析以及后续基因预测。

结果统计：

表3 组装结果Scaffold统计表

sample name	Totle Num	Totle Length	N50 Length	N90 Length	Max Length	Min Length	Sequence GC(%)
sallet_1	149	13788254	261769	449669	811501	502	53.715669873792585
sallet_2	134	13719254	303919	407752	1066020	534	53.7240144398522
sallet_3	113	13720804	343015	500465	1402145	530	53.73521843180619

表4 组装结果Contig统计表

sample name	Totle Num	Totle Length	N50 Length	N90 Length	Max Length	Min Length	Sequence GC(%)
sallet_1	176	13780007	239054	276070	811501	502	53.71772307517696
sallet_2	154	13715172	274956	407752	1066020	534	53.72612899058065
sallet_3	137	13719598	315691	441489	1402145	530	53.736064278268216

注：
表3中统计的是长度大于500bp的Scaffold，表格4中统计的是大于500bp的Scaffold在N处打断得到的Contig序列。

将reads比对到组装好的基因组序列上，通过统计组装序列的GC含量和reads覆盖深度，总结基因组的GC偏向性和重复序列情况。分析结果如下图所示：

sallet_1.GC_depth_countsallet_2.GC_depth_countsallet_3.GC_depth_count

图3 样品GC含量与测序深度（Depth）关联分析统计图

说明：
横坐标GC参与组装的read的GC含量，纵坐标depth为参与组装的read的评价测序深度。上方的柱状图为GC含量对应的柱状图，右侧柱状图为测序深度对应的柱状图。

3.3 基因组组分分析

customer_files/03genomeComponent/
├── 03.01gene_predict/
│   ├── [物种].[样本].augustus.protein.fa		[各样本预测基因对应蛋白质序列]
│   ├── [物种].[样本].augustus.gff			[各样本预测基因gff文件]
│   ├── [物种].[样本].coding_gene_length.svg	[各样本预测基因长度分布图]
│   ├── [物种].[样本]_coding_gene_length.txt	[各样本预测基因长度信息]
│   └── coding_gene_info_count.txt		[各样本预测基因长度汇总信息]
├── 03.02repeat_seq/
│   ├── tandem_repeat_count.svg			[串联重复统计图]
│   ├── tandem_repeat_count.txt			[串联重复统计表]
│   └── [物种].[样本].interspersed_repeat_count.txt	[各样本散在重复序列统计表]
└── 03.03ncRNA_predict/
    ├── ncRNA_count.txt				[ncRNA统计表]
    ├── [物种].[样本].purged.tblout.dealed.table	[ncRNA预测结果]
    ├── [物种].[样本].purged.tRNA.out		[tRNA预测结果]
    └── [物种].[样本].purged.tRNA.structure	[tRNA预测的结构文件]

微生物基因组包含的功能区域非常丰富，除编码基因区域，更有非编码区域实现转录调控、转录后调控、翻译调控、表观遗传调控等功能，部分功能区域还与物种进化的多样性存在关系。 通过多种方法，对编码基因、重复序列、非编码RNA等进行预测，获取目标基因组的组成情况。

3.3.1 编码基因

根据获得的组装序列，我们使用Augustus软件对真菌样本的编码基因进行预测。该软件基于HMM（隐马尔科夫模型）和贝叶斯理论，根据序列信息对其中的编码基因进行预测。基因预测结果统计信息如下表所示：

表4 编码基因预测结果统计表

sample	genome length	gene number	gene totle length	gene average length	gene length/genome(%)
sallet_1	13788254	4698	8121189	1728.648148148148	58.89932837036509
sallet_2	13719254	4671	8075761	1728.914793406123	58.86443242467848
sallet_3	13720804	4717	8064343	1709.63387746449	58.77456598024431

说明：
Genome size：全基因组总长度；
Gene number：预测到的编码基因个数；
Gene total length：所有编码基因的总长度；
Gene average length：编码基因的平均长度；
Gene length / Genome：编码区总长度占全基因组的比例。

绘制基因长度统计图如下：

Diutina_catenulata.sallet_1_coding_gene_lengthDiutina_catenulata.sallet_2_coding_gene_lengthDiutina_catenulata.sallet_3_coding_gene_length

图4 基因长度统计图

说明：
横坐标gene length为预测的基因长度区间，纵坐标gene number为对应长度区间内的基因数量。

3.3.2 重复序列

重复序列是基因组中不同位置出现的相同或互补性片段，是基因调控网络的组成成分。 根据重复的序列在基因组上的分布，分为散在重复序列、串联重复序列。

散在重复序列又分短分散重复序列（Short interspersed nuclear elements，SINEs）以及长散在重复序列（Longinterspersed nuclear elements，LINEs），其中长散在重复序列常具有转座活性。串联重复序列（Tandem Repeat，TR），即相邻的、重复两次或多次特定核酸序列模式的重复序列。分为Minisatellite DNA（小卫星DNA）和Microsatellite DNA（微卫星DNA）。串联重复单元具有种属组成特异性，可作为物种的遗传性状，进行进化关系的研究。

通过RepeatMasker软件进行散在重复序列预测，TRF（Tandem Repeats Finder）搜寻DNA序列中的串联重复序列。

预测结果如下表所示:

表5 重复序列信息统计

说明：
samples：样本名称；
Minisatellite DNA number：Minisatellite DNA数量；
Minisatellite DNA length range：Minisatellite DNA长度范围；
Minisatellite DNA totle length：Minisatellite DNA总长度；
Minisatellite DNA in genome(%)：Minisatellite DNA长度占基因组长度百分比；
Microsatellite DNA number：Microsatellite DNA数量；
Microsatellite DNA length range：Microsatellite DNA长度范围；
Microsatellite DNA totle length：Microsatellite DNA总长度；
Microsatellite DNA in genome(%)：Microsatellite DNA长度占基因组长度的百分比；
other tandem repeat number：other tandem repeat数量；
other tandem repeat length range：other tandem repeat长度范围；
other tandem repeat totle length：other tandem repeat总长度；
other tandem repeat in genome(%)：other tandem repeat长度占基因组长度的百分比。

重复序列统计可视化如下：

图5 重复序列统计图

说明：
横坐标samples为不同的样本，纵坐标count为各个样本中各类重复序列的数量。

3.3.3 非编码RNA

非编码RNA（ncRNA）是一类执行多种生物学功能的RNA分子，其本身并不携带翻译为蛋白质的信息，直接在RNA水平对生命活动发挥作用。对于微生物而言，研究较为普遍的包括sRNA、rRNA、tRNA。

tRNA：通过tRNAscan-SE软件对tRNA进行预测。

rRNA：即核糖体RNA，rRNA在相邻物种中高度保守。rRNA的预测方法有两种，一是通过与近缘参考序列的rRNA库比对找到rRNA，二是用rRNAmmer软件预测rRNA。

sRNA：小RNA，首先进行Rfam database比对注释，接着用cmsearch程序（参数默认）确定最终的sRNA。

snRNA：（small nuclearRNA，小核RNA），它是真核生物转录后加工过程中RNA剪接体(spilceosome）的主要成分。

miRNA：MicroRNA(miRNA)是在真核生物中发现的一类内源性的具有调控功能的非编码RNA，前体全长约90bp，其成熟miRNA大小长约20~25个核苷酸（nt）。miRNA广泛存在于真核生物中，是一组不编码蛋白质的短序列RNA，它本身不具有开放阅读框（ORF），对基因的表达具有调控作用。

sRNA、snRNA、miRNA的预测原理类似，首先用Rfam软件进行Rfam database比对注释，接着用其cmsearch程序（参数默认）确定最终的sRNA、snRNA、miRNA。

表6 ncRNA去冗余后的统计结果

说明：
samples：样本名称；
* number：ncRNA数量；
* average length：ncRNA平均长度；
* totle length：ncRNA总长度。

3.4 基因功能分析

customer_files/04function_annotation/
├── 04.01common_database/
├────── 04.01.00all_result_count
│   	├── [物种].[样本].all_anno_count.svg		[所有常见数据库注释结果汇总图]
│   	└── [物种].[样本].all_anno_count.txt		[各样本预测基因gff文件]
├────── 04.01.01GO_result
│   	├── [物种].[样本].dataPlot.go.alldata.txt		[各样本GO数据库注释结果统计]
│   	├── [物种].[样本].dataPlot.go.txt			[各样本GO数据库注释结果统计(前20个结果)]
│   	└── [物种].[样本].GO_classes.svg			[各样本GO数据库注释结果统计图]
├────── 04.01.02KEGG_result
│   	├── [物种].[样本].dataPlot.kegg.txt			[各样本KEGG数据库注释结果统计]
│   	└── [物种].[样本].Kegg_Classes.svg			[各样本KEGG数据库注释结果统计图]
├────── 04.01.03eggNOG_result
│   	├── [物种].[样本].Eggnog.count.txt			[各样本Eggnog数据库注释结果统计]
│   	├── [物种].[样本].Eggnog.result.emapper.annotations	[各样本Eggnog数据库注释结果(该结果也包含GO和KEGG注释结果)]
│   	└── [物种].[样本].Eggnog_plot.svg			[各样本Eggnog数据库注释结果统计图]
├────── 04.01.04KOG_result
│   	├── [物种].[样本].diamond.KOG.result		[各样本KOG数据库注释结果]
│   	├── [物种].[样本].KOG.fig.count.txt			[各样本KOG数据库注释结果统计]
│   	└── [物种].[样本].KOG_plot.svg			[各样本KOG数据库注释结果统计图]
├────── 04.01.05NR_result
│   	├── [物种].[样本].nr_count.txt			[各样本NR数据库注释结果统计]
│   	├── [物种].[样本].nr_plot.svg			[各样本NR数据库注释结果统计图]
│   	└── [物种].[样本].diamond.nr.result			[各样本NR数据库注释结果]
├────── 04.01.06TCDB_result
│   	├── [物种].[样本].tcdb_plot.svg			[各样本tcdb数据库注释结果统计图]
│   	└── [物种].[样本].diamond.tcdb.result		[各样本tcdb数据库注释结果]
├────── 04.01.07Pfam_result
│   	└── [物种].[样本].Pfam.tblout.result		[各样本Pfam数据库注释结果]
├────── 04.01.08Swiss_Prot_result
│   	└── [物种].[样本].diamond.swissprot.result		[各样本swissprot数据库注释结果]
├────── 04.01.09CAZy_result
│   	├── [物种].[样本].cazy_plot.svg			[各样本cazy数据库注释结果图]
│   	└── [物种].[样本].CAZyme.table			[各样本cazy数据库注释结果]
├── 04.02Effector/
│   ├── [物种].[样本].antismash_result				[各样本次级代谢基因簇分析结果]
│   ├── [物种].antiSMASH_count_type_number.txt				[各样本次级代谢基因簇及基因数量统计]
│   ├── [物种].Effector_T3SS_count.txt				[T3SS效应蛋白预测结果统计]
│   ├── [物种].Effector_TNSS_count.txt				[TNSS效应蛋白预测结果统计]
│   ├── [物种].[样本].antiSMASH_count_type_number.svg				[各样本各类PKS基因簇及基因数量统计]
│   ├── [物种].[样本].protein.oneline_summary.signalp5		[各样本分泌蛋白预测结果]
│   ├── [物种].[样本].diamond.p450.result			[各样本P450数据库注释结果]
│   ├── [物种].[样本].tmhmm.txt				[各样本分泌蛋白预测结果]
│   ├── [物种].signal_table.txt				[分泌蛋白预测结果]
│   └── [物种].[样本].TNSS_count.txt				[各样本分泌系统蛋白及T3SS效应蛋白预测结果]
└── 04.03Pathogenicity_analysis/
    ├── [物种].[样本].diamond.DFVF.result	[各样本DFVF数据库注释结果]
    ├── [物种].[样本].diamond.PHI.result		[各样本PHI数据库注释结果]
    ├── [物种].[样本].PHI_count.svg		[各样本PHI数据库注释结果统计图]
    └── [物种].[样本].PHI_count.txt		[各样本PHI数据库注释结果统计]

目前提供注释的通用功能数据库主要有GO、KEGG、KOG、NR、Pfam和Swiss-Prot。

功能注释基本步骤如下：

1）将预测基因的蛋白序列与各功能数据库进行Diamond 比对（evalue ≤ 1e-5）；

2）比对结果过滤：对于每一条序列的比对结果，选取 score 最高的比对结果进行注释。

本项目进行的编码基因的注释结果统计如下图所示：

Diutina_catenulata.sallet_1.all_anno_countDiutina_catenulata.sallet_2.all_anno_countDiutina_catenulata.sallet_3.all_anno_count

图6 基因功能分析结果统计图

说明：
横坐标database为参与注释的各个数据库，纵坐标number of gene为各个数据库注释出来的基因数量。

3.4.1 常用数据库

3.4.1.1 GO数据库注释

GO的全称是Gene Ontology，是一套国际标准化的基因功能描述的分类系统。GO分为三大类：1）细胞组分（Cellular Component）：用于描述亚细胞结构、位置和大分子复合物，如核仁、端粒和识别起始的复合物；2）分子功能（Molecular Function）：用于描述基因、基因产物个体的功能，如与碳水化合物结合或 ATP 水解酶活性等；3）生物过程（Biological Process）：用来描述基因编码的产物所参与的生物过程，如有丝分裂或嘌呤代谢等。

GO数据库三大分类统计结果如下图：

sallet_1.GO_classessallet_2.GO_classessallet_3.GO_classes

图7 GO数据库分类统计图

说明：
横坐标GO Term为注释到的GO Term，结果数量太多的只显示前20个。纵坐标Count为注释到各个GO Term的基因数量。图中不同颜色分别对应GO数据库中的BP、CC、MF三个分类。

3.4.1.2 KEGG数据库注释

KEGG全称为Kyoto Encyclopedia of Genes and Genomes。系统分析基因产物和化合物在细胞中的代谢途径以及这些基因产物的功能的数据库。它整合了基因组、化学分子和生化系统等方面的数据,包括代谢通路(KEGG PATHWAY)、药物(KEGG DRUG)、疾病(KEGG DISEASE)、功能模型(KEGG MODULE)、基因序列(KEGG GENES)及基因组(KEGG GENOME)等等。KO(KEGG ORTHOLOG)系统将各个KEGG注释系统联系在一起，KEGG已建立了一套完整KO注释的系统，可完成新测序物种的基因组或转录组的功能注释。详见http://www.genome.jp/kegg/。

绘制KEGG数据库中注释基因数目统计图如下：

sallet_1.Kegg_Classessallet_2.Kegg_Classessallet_3.Kegg_Classes

图8 KEGG数据库分类统计图

说明：
横坐标count(%)为注释到各个KEGG pathway的基因占所有基因中的百分比，纵坐标为注释到的KEGG pathway。

3.4.1.3 eggNOG数据库注释

eggNOG数据库，全称是evolutionary genealogy of genes: Non-supervised Orthologous Groups，是一个蛋白聚类数据库，带有功能描述和功能类别说明，由EMBL（欧洲分子生物实验室）维护。包含1,133个物种、721,801个直系同源组、41个不同水平的直系同源组分类，整合了5,214,234个蛋白序列。分别更新了4,873个COG数据库信息和4,850个KOG数据库信息。

eggNOG数据库按照功能一共可以分为二十五类，其统计结果如下图：

sallet_1.Eggnog_plotsallet_2.Eggnog_plotsallet_3.Eggnog_plot

图9 eggNOG数据库分类统计图

说明：
横坐标eggNOG classes为eggNOG数据库的不同分类，纵坐标count为注释到eggNOG数据库不同分类的基因数量。

3.4.1.4 KOG数据库注释

KOG数据库，属于COG数据库的一个针对真核生物的直系同源数据库。

COG，全称是Cluster of Orthologous Groups of proteins，由NCBI创建并维护的蛋白数据库，根据细菌、藻类和真核生物完整基因组的编码蛋白系统进化关系分类构建而成。通过比对可以将某个蛋白序列注释到某一个COG中，每一簇COG由直系同源序列构成，从而可以推测该序列的功能。COG数据库按照功能一共可以分为二十五类，详见http://www.ncbi.nlm.nih.gov/COG/。

COG数据库按照功能一共可以分为二十五类，其统计结果如下图：

sallet_1.KOG_plotsallet_2.KOG_plotsallet_3.KOG_plot

图10 KOG数据库分类统计图

说明：
横坐标KOG classes为KOG数据库不同分类，纵坐标count(%)为注释到不同KOG classes基因占所有基因的百分比。

3.4.1.5 NR数据库注释

NR全称为Non-Redundant Protein Database，是一个非冗余的蛋白质数据库，由NCBI创建并维护，其特点在于内容比较全面，同时注释结果中会包含有物种信息，可作物种分类用。根据基因注释到的物种情况，统计注释到的物种及基因数目，其统计结果如下图:

sallet_1.nr_plotsallet_2.nr_plotsallet_3.nr_plot

图11 NR数据库物种比对统计图

说明：
横坐标species为基因在NR数据库注释到的不同物种，纵坐标count为注释到不同物种的基因数量。

3.4.1.6 TCDB数据库注释

TCDB，全称是Transporter Classification Database，转运蛋白分类数据库，是膜转运蛋白，包括离子通道（ion channels）的分类系统（TC system）。TCDB数据库转移系统以5个级别进行分类，第一级统计结果如下图：

sallet_1.tcdb_plotsallet_2.tcdb_plotsallet_3.tcdb_plot

图12 TCDB数据库功能分类统计图

说明：
横坐标Function classes为TCDB数据库的不同功能分类，纵坐标Number of matched genes为注释到不同功能分类的基因数量。

3.4.1.7 Pfam数据库注释

蛋白质一般由一个或多个功能区构成，这些区通常被称为域。结构域的不同组合方式产生的蛋白质在自然界中各种不同。因此蛋白结构域的鉴别对分析蛋白质的功能来说尤其重要。Pfam数据库有两个组成部分：Pfam-A和Pfam-B，其中Pfam-A经过人工筛选，质量较高。详见http://pfam.xfam.org/。

该分析对预测基因对应的蛋白质序列进行Pfam注释比对，获得每个蛋白序列的结构域信息，结果在customer_files/04function_annotation/04.01common_database/04.01.07Pfam_result文件夹下，文件名格式为：[物种].[样本].Pfam.tblout.result。

3.4.1.8 Swiss-Prot数据库注释

Swiss-Prot是一个精选的蛋白质序列数据库，它提供一个高水平的注释结果，例如一个蛋白质功能、其域结构、翻译后修饰、变异等的描述。详见http://www.ebi.ac.uk/uniprot/。

我们使用diamond软件对预测的蛋白质序列和Swiss-Prot数据库提供的蛋白质序列进行比对，从而对每个蛋白质序列进行注释，结果在customer_files/04function_annotation/04.01common_database/04.01.08Swiss_Prot_result文件夹下，文件名格式为：[物种].[样本].diamond.swissprot.result。

3.4.1.9 碳水化合物活性酶（CAZy）数据库注释

CAZy全称为Carbohydrate-Active enZYmes Database，碳水化合物酶相关的专业数据库，内容包括能催化碳水化合物降解、修饰、以及生物合成的相关酶系家族。其包含五个主要分类：糖苷水解酶 （Glycoside Hydrolases, GHs）、糖基转移酶（GlycosylTransferases, GTs）、多糖裂解酶（Polysaccharide Lyases, PLs）和糖类酯解酶（Carbohydrate Esterases, CEs）、氧化还原酶（Auxiliary Activities, AAs）。

碳水化合物结合结构域是一种非催化结构域，能折叠成特定的三维空间结构，具有结合碳水化合物的功能。近年来研究表明：碳水化合物结合结构域能通过结合碳水化合物活性酶的底物，提高碳水化合物活性酶的催化结构域作用于底物的活性。

CAZy数据库分类注释结果个数统计图展示如下：

sallet_1.cazy_plotsallet_2.cazy_plotsallet_3.cazy_plot

图13 CAZy数据库分类统计图

说明：
横坐标Cazy Classes为Cazy数据库不同分类，纵坐标Number of matched genes为注释到各个分类的基因数量。

3.4.2 效应子

3.4.2.1 分泌蛋白预测

分泌蛋白是指在细胞内合成后，在信号肽的引导下穿过细胞膜分泌到细胞外起作用的蛋白质。分泌蛋白中有许多是生命活动所需的重要酶类。分泌蛋白的N端是由15～30个氨基酸组成的信号肽，对分泌蛋白的分泌起主导作用。

使用SignalP、TMHMM工具进行预测，检测是否含有信号肽及跨膜结构，综合预测蛋白序列是否是分泌蛋白。

表7 分泌蛋白预测结果

sample ID	Signal Protein	TMHMM Protein	Secreted Protein
sallet_1	372	633	70
sallet_2	372	619	69
sallet_3	367	632	75

说明：
sample ID：样本名称；
Signal Protein：Signal软件预测的含有信号肽的蛋白质数量；
TMHMM Protein：TMHMM软件预测的含有跨膜结构的蛋白质数量；
Secreted Protein：综合预测为分泌蛋白的数量。

3.4.2.2 分泌系统蛋白及T3SS效应蛋白预测

病原菌通过分泌系统TNSS（type N secretion systems，目前确定的有7种，I型-VII型）将该类蛋白分泌至胞外或是宿主细胞中，通过控制免疫应答反应以及细胞衰亡引起病理反应，而其中革兰氏阴性菌的T3SS通常用来从分子水平研究病原菌，感染机制，毒力作用等，是研究得比较多的分泌系统。

对于TNSS系统，通过蛋白序列功能数据库注释结果中，提取分泌系统相关蛋白进行注释。对于革兰氏阴性菌，另外采用EffectiveT3 软件预测T3SS效应蛋白。

表8 TNSS结果统计

sample ID	totle gene	T1SS num	T2SS num	T3SS num	T4SS num	T5SS num	T6SS num	T7SS num
sallet_1	4698	0	0	0	0	0	0	0
sallet_2	4671	0	0	0	0	0	0	0
sallet_3	4717	0	0	0	0	0	0	0

说明：
sample ID：样本名称；
totle gene：所有基因数量；
T1SS num：T1SS数量；
T2SS num：T2SS数量；
T3SS num：T3SS数量；
T4SS num：T4SS数量；
T5SS num：T5SS数量；
T6SS num：T6SS数量；
T7SS num：T7SS数量；

表9 T3SS效应蛋白预测结果统计

sample ID	Totle gene	T3SS effective true	T3SS effective false
sallet_1	4698	0	4698
sallet_2	4671	0	4671
sallet_3	4717	0	4717

说明：
sample ID：样本名称；
Totle gene：所有基因数量；
T3SS effective true：预测为T3SS效应蛋白的数量；
T3SS effective false：预测不是T3SS效应蛋白的数量；

3.4.2.4 次级代谢基因簇分析

次级代谢产物是微生物在一定的生长时期，以初级代谢产物为前体合成的对微生物的生命活动无明确功能，并非生长繁殖所必需的物质。采用antiSMASH程序对基因组进行预测。

PKS可分为三种类型：I型也成为模块类PKS，是由s多个结构域组成的多功能酶复合物。II型也成为芳香类PKS，主要合成芳香类化合物。III型也成查尔酮型PKS。使用antiSMASH程序对基因组进行预测

表10 次级代谢基因簇及基因数量统计

说明：
sample ID：样本名称；
Totle gene：所有基因数量；
T3SS effective true：预测为T3SS效应蛋白的数量；
T3SS effective false：预测不是T3SS效应蛋白的数量；

Diutina_catenulata.sallet_1.antiSMASH_count_type_numberDiutina_catenulata.sallet_2.antiSMASH_count_type_numberDiutina_catenulata.sallet_3.antiSMASH_count_type_number

图14 各类PKS基因簇及基因数量统计

说明：
横坐标PKS type为不同PKS类型，纵坐标count为注释到不同PKS类型的基因数量。

3.4.2.3 P450数据库注释

细胞色素P450（cytochromeP450或CYP450,简称CYP450）为一类亚铁血红素—硫醇盐蛋白的超家族，它参与内源性物质和包括药物、环境化合物在内的外源性物质的代谢。

我们使用diamond软件对预测的蛋白序列进行P450数据库注释，获得每个蛋白序列对应的P450信息，结果在04.02Effector/文件夹下，文件名格式为：[物种].[样本].diamond.p450.result。

3.4.3 致病性分析

3.4.3.1 病原与宿主互作数据库（PHI）注释

PHI全称为Pathogen Host Interactions Database，病原与宿主互作数据库，主要来源于真菌、卵菌和细菌病原，感染的宿主包括动物、植物、真菌以及昆虫。该数据库对寻找药物干预的靶基因研究有重要作用，同时该数据库还包括抗真菌化合物和相应的靶基因。数据库中的每个基因都包含核酸和氨基酸序列，以及感染宿主过程中预测的蛋白功能的详细描述。

病原体PHI表型突变类型基因数目的统计情况如下图所示：

sallet_1.PHI_countsallet_2.PHI_countsallet_3.PHI_count

图15 病原体PHI表型突变类型统计

说明：
横坐标mutation type为不同的病原PHI表型突变类型，纵坐标number of genes为注释到不同突变类型的基因数量。

3.4.3.2 真菌毒力因子数据库（DFVF）注释

DFVF数据库全称为database of fungal virulence factors（真菌毒力因子数据库），是一个综合的已知真菌毒力因子数据库，收集了来自85个属的228个真菌菌株所产生的2058个致病基因。每个基因详细描述引起的疾病和作用的宿主，更与Pfam功能域注释和GO注释信息进行了关联。

我们使用Diamond软件，把目标物种的氨基酸序列，与DFVF数据库进行比对，把目标物种的基因和其相对应的功能注释信息结合起来，得到注释结果，结果在04.03Pathogenicity_analysis/文件夹下，文件名格式为：[物种].[样本].diamond.DFVF.result。

3.5 比较基因组分析

customer_files/05Comparative_Genomics_Analysis
├── 05.01Collinearity_analysis/
│   ├── [样本].minimap.paf		[各样本比较基因组分析结果]
│   ├── [样本].dot.png		[各样本共线性分析点图]
│   ├── [样本].[svg|png]		[各样本共线性图]
│   └──	[样本].chr.[svg|png]		[各样本共线性图]
├── 05.02SNP/
│   ├── snp_count.txt		[各个区域SNP统计]
│   └── [样本].snp						[各样本SNP识别结果]
├── 05.03Indel/
│   ├── [样本]_indel.result		[各样本Indel识别结果]
│   └── Indel_count.txt			[各区域Indel统计]
└── 05.04SV/
    ├── [样本]_sv.SVs_all.tsv		[各样本SV识别结果]
    ├── YGTC_sv_plot_top.svg		[各样本SV识别结果circos简约图]
    └── [样本]_sv_plot.svg					[各样本SV识别结果circos图]

3.5.1 全基因组共线性分析

共线性指的是遗传学中的基因连锁关系， 是不同物种基因组上同源基因以相同顺序排列的现象。两个物种之间的共线性程度可以作为衡量他们之间进化距离的尺度，可以知道物种间的亲缘关系。

使用MUMmer软件对样本基因组和参考基因组进行比对，确定样本基因组和参考基因组之间的比对关系。然后使用SyRI软件识别样本基因组和参考基因组之间的共线性区域。

样品基因组和参考基因组之间的共线性展示结果如下图：

Diutina_catenulata.sallet_1.plotDiutina_catenulata.sallet_2.plotDiutina_catenulata.sallet_3.plot

图16 基因组共线性线图

说明：
横坐标和纵坐标分别表示比对的两个基因组位置，图中的连线为两个基因组共线性部分。

Diutina_catenulata.sallet_1.plotsrDiutina_catenulata.sallet_2.plotsrDiutina_catenulata.sallet_3.plotsr

图17 基因组共线性比对图

说明：
横坐标Chromosome postion(in Mbp)表示基因组上的位置，单位为Mb。纵坐标Reference Chromosome ID为参考基因组的不同contig，每个contig ID对应的上下两条横线分别表示参考基因组和对比基因组的contig。

3.5.2 SNP统计与注释

SNP（单核苷酸多态性）主要是指在基因组水平上由单个核苷酸的变异所引起的DNA序列多态性，包括单个碱基的转换、颠换等。

采用MUMmer比对软件进行个体SNP的检测，依照SNP和基因之间的位置关系及相互作用，对SNP功能进行注释。

表11 SNP结果统计

说明：
sample：样本名称；
snp in start_codon：位于起始密码子的SNP数量；
snp in intron：位于内含子的SNP数量；
snp in CDS：位于CDS区的SNP数量；
snp in stop_codon：位于终止密码子的SNP数量；
snp not in gene：位于非基因区的SNP数量。

3.5.3 InDel统计与注释

InDel是指基因组中小片段的插入和缺失序列。

利用LASTZ软件检测长度小于50bp的小片段插入与缺失（InDel），对InDel类型进行统计，结果见下表：

表12 InDel类型统计

说明
sample：样本名称；
indel in start_codon：位于起始密码子的indel数量；
indel in intron：位于内含子的indel数量；
indel in CDS：位于CDS区的indel数量；
indel in stop_codon：位于终止密码子的indel数量；
indel not in gene：位于非基因区的indel数量。

3.5.4 SV统计与注释

基因组结构变异（Structural Variation，SV）一般是指基因组上大长度的序列变化和位置关系变化。基因组结构变异包含很多种类型，通常定义是长度大于50bp的插入（Insertion）、缺失（Deletion）、串联重复（Tandem repeate）、染色体倒位（Inversion）、染色体内部或染色体之间的序列易位（Translocation）、拷贝数变异（CNV）以及形式更为复杂的嵌合性变异。其中，占比最大的就是大片段的插入删除。

使用MUMmer软件对样本基因组和参考基因组进行比对，确定样本基因组和参考基因组之间的比对关系。然后使用SyRI软件识别样本基因组中的结构变异信息。

全基因组结构变异类型配对图如下所示:

sallet_1.svsallet_2.svsallet_3.sv

图18 结构变异可视化

说明
红色圆环：样本基因组；
绿色圆环：参考基因组；
橙色连线：大片段插入；
粉色连线：易位；
浅绿色连线：大片段删除；
浅褐色连线：串联重复收缩；
紫色连线：串联重复复制；
灰色连线：倒置；

3.6 群体进化分析

customer_files/06population_evolution/
├── 06.01cor_pan_genes/
│	├── [物种].cd_hit.result.clstr	[各样本蛋白质聚类结果]
│	└── Flower_plot.svg			[各样本共有基因及特意基因花瓣图]
└── 06.02Hcluster/
	├── gene_family_bar_count.txt		[各样本基因家族柱状图数据]
	├── gene_family_count.txt			[各样本基因家族统计]
	├── orthologs_gene_bar_plot.svg		[同源基因柱状图]
	├── Phylogenetic_Tree.svg			[系统发生树]
	└── Phylogenetic_Tree_dis.svg		[系统发生树带距离]

3.6.1 共有和特有基因分析

所有样本中均存在的同源基因称为共有基因（Core gene），除去共有基因，其他的基因称为非共有基因（Dispensable gene），特有基因（Specific gene）是只有在某个样品中所特异拥有的基因。其中共有基因（Core gene）和特有基因（Specific gene）很可能与样品的共性和特性相对应，可以作为样本间功能差异的研究依据。

使用cd-hit软件对需要分析的多个样品的蛋白序列进行聚类，并用R进行绘图。

通过比较多个样品的蛋白质序列，我们统计了它们的共有基因和特有基因，以花瓣图形式展示：

图19 共有基因和特有基因花瓣图

说明：
花瓣图中不同花瓣表示各个样本中特异的基因数量，花瓣中心为所有样本共有基因数量。

3.6.2 基因家族分析

基因组进化中，一个基因通过基因重复产生了两个或更多的拷贝，这些基因即构成一个基因家族。基因家族是具有共同祖先的一组基因，家族内不同基因往往具有相似的结构和功能。

使用diamond软件对多个目标基因组的蛋白进行两两比对，过滤比对不可信的结果，再使用Solar（Version 0.9.6）去除冗余，用Hcluster-sg按照比对相似度对蛋白进行聚类，得到基因家族聚类结果。

表14 基因家族鉴定的统计结果

samples	totle gene number	gene number in families	uncluster genes number	family number	unique family number
sallet_1	4698	4663	35	4525	4422
sallet_2	4671	4649	22	4527	4437
sallet_3	4717	4513	204	4351	4230

说明：
samples：样本名称；
totle gene number：所有基因数量；
gene number in families：可归到基因家族中基因数量；
uncluster genes number：未归到聚类中基因数量；
family number：样本中基因家族数量；
unique family number：唯一家族数量

图20 同源基因数目统计条形图

说明：
横坐标samples为各个样本名称，纵坐标Number of genes为各个样本不同基因家族中的基因数量。

图21 多个物种/样本的同源基因家族数目Venn图

3.6.3 物种进化分析

系统发生树（英文：phylogenetic tree或evolutionary tree）是表明被认为具有共同祖先的各物种相互间演化关系的树。 它用来表示系统发生研究的结果，用它描述物种之间的进化关系。

用Treebest（Version 1.9.2）（Neighbor-Joining,NJ）或 PHYML（Maximum likelihood,ML)（Version v3.0）软件构建进化树，物种之间的进化树见下图：

图23 物种进化系统发生树

三、所用软件的版本

软件	版本
Trimmomatic	0.39
FastQC	0.11.9
jellyfish	2.2.10
bedtools	2.30.0
canu	2.2
samtools	1.7
pbmm2	1.9.0
minimap2	2.15
RepeatModeler	2.0.3
tRNAscan-SE	2.0.9
diamond	0.8.22
SignalP	5.0b
TMHMM	2.0
antiSMASH	6.1.1
MUMmer4	4.0.0
lastz	1.04.15
OrthoFinder	2.5.4
ipdSummary	3.0

四、 参考文献

1. Trimmomatic: a flexible trimmer for Illumina sequence data. (Trimmomatic)

2. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. (jellyfish)

3. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. (canu)

4. The Sequence Alignment/Map format and SAMtools. (samtools)

5. Minimap2: pairwise alignment for nucleotide sequences. (Minimap2)

6. tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. (tRNAscan-SE)

7. Sensitive protein alignments at tree-of-life scale using DIAMOND. (diamond)

8. SignalP 5.0 improves signal peptide predictions using deep neural networks. (SignalP)

9. Predicting Transmembrane Protein Topology with a Hidden Markov Model: Application to Complete Genomes. (TMHMM)

10. antiSMASH: rapid identification, annotation and analysis of secondary metabolite biosynthesis gene clusters in bacterial and fungal genome sequences. (antiSMASH)

11. MUMmer4: A fast and versatile genome alignment system. PLoS computational biology. (MUMmer4)

12. Improved pairwise alignment of genomic DNA. (LASTZ)

13. OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy. (OrthoFinder)