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

微科盟细菌基因组完成图分析结题报告

---

一 概述

细菌基因组完成图是采用三代测序平台，一次性获得完整的细菌基因组图谱，实现没有gap和N碱基的组装，同时可以获得质粒基因组信息。二代测序技术难以解决细菌基因组中重复区域与高GC区域的组装的问题，因此往往需要借助PCR扩增、Sanger测序等方式填补组装后基因组中的gap。而基于三代测序技术充分发挥了超长读长的优势，轻松跨越细菌基因组中的这些复杂结构，不但在组装流程上更为简洁，还能够对组装结果更加可靠稳定。

---

二 项目流程

2.1 测序实验流程

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

经电泳检测合格的DNA样品用Covaris超声波破碎仪随机打断成长度约为350bp的片段。处理完成后的DNA片段，使用NEBNext®Ultra™ DNA Library Prep KitforIllumina(NEB, USA)试剂盒，经末端修复、加A尾、加测序接头、纯化、PCR扩增等步骤完成整个文库制备。文库构建完成后，先使用Qubit 2.0进行初步定量，稀释文库至2ng/ul，随后使用Agilent 2100对文库的插入片段进行检测，insert size符合预期后，使用q-PCR方法对文库的有效浓度进行准确定量，以保证文库质量。

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

2.2 生信分析流程

详细过程见结果解读部分

---

三 结果文件夹解读

3.1 质控

01.Data_qc/
├──illumina_fastp.xls   [二代测序数据质控结果] 
├──pacbio_sum.xls[[三代测序数据longQC质控结果]

原始序列数据

高通量测序得到的图像经 Base Calling 转化为原始测序序列（Reads），我们称之为 Raw Data 或 Raw Reads，结果以 FASTQ（简写为 fq）文件格式储存，它包含测序序列（Reads） 的序列信息及其对应的测序质量信息。PE 文库的数据结果中，每个样品均有测序两端的 Reads， 并且 Reads 的顺序是严格一致、相互对应的。 FASTQ 文件中每条 read 都由四行构成，文件格 式如下：

其中，第 1 行和第 3 行为 read 名称（后来为了节省储存空间省略掉第 3 行“+”后面的序列 名称），由 Illumina 测序仪产生；第 2 行是 read 的碱基序列，第 4 行是 read 中每个碱基对应 的测序质量分数，由 Base Calling 转化而来，每个字母对应的 ASCII 值减去相应测序质量系 统的 Phred quality score（33 或 64），即为该碱基的测序质量值（Phred 值）。

在获得每个样品二代测序数据之后，首先需要对测序的数据质量进行评估并对低质 量的数据进行去除，以保证后续分析结果的可信度，以上步骤称之为原始数据的质量控制。 质量控制获得的高质量序列则用于下游的数据分析。质控流程采用软件fastp^[2]，具体处理步骤如下：

1. 去除带接头（adapter）的 paired reads；
2. 当单端测序 read 中含有的低质量（ sQ <= 20）碱基数超过该条 read 碱基总数的 20%时，去除此 paired reads;
3. 去除 PCR 扩增产生的重复 reads。

测序数据质控前后的质量信息统计表格如下：

三代测序平台Pacbio Sequel II测序仪的下机数据，我们称之为subreads，它们是已经去除了adapter和barcode序列的数据。
我们用bam2fastq^[1]软件，从测序仪下机数据subreads.bam文件提取出fastq格式的数据文件，里面所有碱基的碱基质量值都为0，这是Pacbio Sequel II平台的测序特点，即由于出现测序错误的碱基几乎是随机分布的，不像Illumina测序那样具有位置倾向性，无法提供有实际意义的单个碱基的测序质量值。
我们这里的质控处理，主要是用fastp^[2]软件把长度偏小的subreads给剔除掉，而不做局部的序列裁剪。在三代测序文库的构建中，一般每一条待测的生物样本序列会在10kb以上，如果subread测出来的序列太短，则有可能不是想测的生物样本序列，而是实验过程中引入的其它杂质序列。所以我们在质控步骤中把长度小于0.5kb的subreads给剔除掉，这样保留下来的数据我们称之为clean data。
最后我们用LongQC^[3]软件对clean data做质控分析报告。

---

3.2 基因组组装

02.Assemble/
├──all_sample.assembly_stat.xls  [所有样本组装统计结果]
├──all_sample_genomesize_stat.xls [所有样本基因组大小结果]
├──all_sample_gc_depth_stat.xls [所有样本GC_depth统计结果]
├──*.final.assembly.fasta [样本组装序列]
├──*.gfastats.tab  [样本N50等组装评估结果]
└──*.rfplasmid.plasmid.csv    [RFPlasmid结果文件]

使用unicycler^[5]软件（默认参数：--keep 0 --min_fasta_length 1000）,以质控后的clean data作为输入数据,做基因组组装。 组装出来的基因组序列文件请查看*final.assembly.fasta文件。 为了鉴别组装出来是质粒还是染色体，使用RFPlasmid^[6]软件（默认参数：--species Generic --jelly）来区分质粒和染色体。 对于每个样本最终得到的contigs, 我们使用quast^[4]软件进行质量评估

注：有些细菌组装出来多条染色体，可能原因有

• 1.该菌株确实有多条染色体（或巨质粒）；
• 2.该菌株是复杂基因组；
• 基因组大小 >7M；染色体基因组个数 + 质粒个数 > 4； 转座子或其它重复序列在基因组中的比例异常高 ( > 10% ) 最长重复序列大于9K；GC含量异常高等可以称为复杂基因组
• 3.该菌株有污染等

基因组大小是指单倍体细胞核中的所含的DNA的总量。预测未知基因组大小的方法可以通过 Illumina 测序数据的 k-mer 分析进行估计。 我们使用利用Jellyfish^[8]对reads进行处理，获得不同频率的k-mer信息，然后利用相应的软件进行基因组大小预估。而Kmerfreq利用自身算法进行K-mer剪切（K-mer=17）统计，利用K-mer频率分布使用GenomeScope^[9]进行基因组大小预估。

对组装的基因组序列以500bp为windows，无重复计算片段的平均GC含量和平均深度并作图。基于每一个windows对应的平均GC和平均深度进行绘图。可以根据此图分析测序数据是否存在GC偏向性以及样本是否存在污染。 我们使用bwa^[7]和相应脚本对组装后的基因组序列和质控后测序序列进行分析，详细统计结果见下表。

---

3.3 基因预测

03.Gene_predict/
├── all_sample.prokka.stat.txt [所有样本prokka注释结果]
├── COGs-Functional_categories-ncbi_file-delete_numbers.html [COG解释]
├── cog
   ├──sample*.cog.txt            [COG比对结果]
   ├──sample*.merged_summary.txt [COG汇总结果]
   ├──sample*.summed_up_summary.txt [COG汇总结果]
   ├──plot               [COG图片]
      ├──cog_summed_up_summary_cog_color.pdf
      ├──cog_summed_up_summary_cog_color.png
      ├──cog_summed_up_summary_ggplot2_color.pdf
      ├──cog_summed_up_summary_ggplot2_color.png
├── EC_number
   ├──sample*.ec_number.txt  [ec_number结果]
└─  prokka [prokka结果]
   ├──sample*.err
   ├── sample*.faa
   ├──sample*.faa.idx
   ├──sample*.ffn
   ├──sample*.fna
   ├──sample*.fsa
   ├──sample*.gbk
   ├──sample*.gff
   ├──sample*.pdf
   ├──sample*.png
   ├──sample*.sqn
   ├──sample*.tbl
   ├──sample*.tsv
   ├──sample*.txt
   ├──sample*_summarize.txt

prokka注释结果

使用细菌基因组功能注释工具prokka,对组装生成的assembly.fasta基因组序列文件,做基因预测和功能注释,然后用R语言分别对不同的功能注释做统计。

COG,即Clusters of Orthologous Groups of proteins（直系同源蛋白簇）。COG是由NCBI创建并维护的蛋白数据库,根据细菌、藻类和真核生物完整基因组的编码蛋白系统进化关系分类构建而成。每一簇COG由直系同源蛋白序列构成,通过比对可以将某个蛋白序列归到某一个已知功能的COG中,从而可以推测该未知功能的蛋白的功能信息。 我们对基因组prokka注释结果中的.tsv文件,提取COG注释信息,然后用R语言做分类统计。

说明：
ftype,表示function_type。
CDS,表示CoDing_Sequence蛋白编码区。
gene,表示注释到具体的基因英文单词。
COG,表示Clusters of Orthologous Groups直系同源基因簇。
EC_number,表示Enzyme_Commission_number酶学委员会编号。
tRNA,表示转运RNA。
rRNA,表示核糖体RNA。
repeat_region,表示重复序列区
tmRNA,表示既有tRNA的功能性质又有mRNA的功能性质的RNA。

EC_number是酶学委员会（Enzyme Commission）为酶所设计的一套编号分类法,每一个酶的编号都以字母“EC”开头,接着以四段号码来表示,这些号码表示对酶作出不同层级的分类。例如三肽胺基-蛋白酶的编号为EC3.4.11.4,当中的“EC3”是指水解酶（即需要用水来将其它分子分解的酶）；“EC3.4”是指与肽键作用的水解酶；“EC3.4.11”是指从多肽中分开胺基末端的水解酶；“EC3.4.11.4”则是指从三肽中分开胺基末端的水解酶。 我们对基因组prokka注释结果中的.tsv文件,提取EC_number注释信息,整理成单独的表格文件.

---

3.4 菌种鉴定

04.Bac_is_species/
    ├──all_bac_taxid_16s.txt [所有样本16s统计结果]
    ├──16s_ani_sum.tsv  [所有样本基因组16sani结果]
    ├──sample* 
        ├──16s_ANI[ANI结果]
            ├──aai_heatmap.svg  [aai结果图]
            ├──aai_matrix.tsv
            ├──aai_summary.tsv
            ├──ANIb_percentage_identity.pdf
            ├──ANIb_percentage_identity.png
            ├──ANIb_percentage_identity.tab
            ├──TETRA_correlations.pdf
            ├──TETRA_correlations.png
            └──TETRA_correlations.tab

16s序列物种鉴定

使用Prokka预测的所有16s序列结果与SILVA-16S数据库进行blastn得到结果。SILVA数据库有超过10万个手动检查的序列，是细菌分类经典数据库。16S rRNA进行blastn推荐物种界限是98.6%，这里我们使用的为99%。 另外，如果有16S序列起源于超过1物种，则被视为受污染的基因组。

ANI(average nucleotide identity)、AAI(amino acid identity)和TETRA(tetra-nucleotide signature)是计算比较基因组学常用算法，可被用于区分物种。ANI是平均核苷酸相似度，是在核苷酸水平比较两个基因组亲缘关系的指标，其在近缘物种之间有较高的区分度。类似的，AAI是氨基酸一致性，是在氨基酸水平比较两个基因组亲缘关系的指标。而TETRA是统计四核苷酸序列(tetra-nucleotide)的频率。当两个基因组序列相似时，这些四核苷酸频率的相关性越高。因此，两个基因组序列之间四核苷酸频率相关性可以粗略的用于确定两个基因组的基因组相关性。 三者物种界限推荐为：95%(ANI)、95%(AAI)和99%(TERAE).

为了进一步鉴定细菌物种，使用物种分类号（默认：16s物种分类）从NCBI上下载同分类号的基因组序列（下载最多十个）。将组装好的序列与下载的基因组进行比较ANI,AAI和TETRA,以确定组装出来的基因组物种归属。 ANI和TETRA是通过pyani^[12]得到的(使用默认参数：-m ANIb 和-m TETRA)，AAI是通过CompareM^[11]得到的(使用默认参数：aai_wf)

---

3.5 基因组组分分析

3.5.1 噬菌体分析

整合在宿主基因组上的温和噬菌体的核酸序列称之为前噬菌体(Prophage),能随宿主细菌DNA进行同步复制或分裂传代。
能与细菌DNA发生重组交叉的特异位点的DNA序列称为attB（att源自attachment）。插入序列左边的一个attB记为attL,插入序列右边的一个attB记为attR。

我们使用PhiSpy^[16]分析工具,对细菌基因组prokka注释结果中的.gbk格式文件做分析,得到前噬菌体序列分析结果。

---

3.5.2 crispris序列分析

cripris
├── all_sample.crispr.xls [所有样本cripris注释结果]
├── sample
    ├── *_crisprcas.gff3
    ├── *.CRISPR_information.xls
    ├── *.CRISPR_sequence.xls
    ├── *.CRISPR_stat.xls
    ├── Chr_1.gff
    ├── done
    ├── rawCas.fna
    ├── rawCRISPRs.fna
    ├── result.json
    ├── TSV
    │   ├── Cas_REPORT.tsv
    │   ├── Cas_REPORT.xls
    │   ├── CRISPR-Cas_summary.tsv
    │   ├── CRISPR-Cas_summary.xls
    │   ├── Crisprs_REPORT.tsv
    │   └── Crisprs_REPORT.xls
    └── Visualization
        ├── crispr.css
        └── index.html

CRISPRCas序列注释结果

我们使用CRISPRCasFinder^[15]分析工具,对unicycler的组装结果assembly.fasta文件做分析,得到CRISPR-Cas序列分析结果。

---

3.5.3 插入序列分析

插入序列是编码转座所需的酶的一种转座子，它的两侧是短反向末端重复序列。转座子插入的靶序列在插入过程中被复制，在转座子两端先形成两个短正向重复序列。正向重复序列（DR，direct repeat）的长度为5-9bp，是任一转座子的特征。IS是细菌染色体和质粒的正常组成成分。标准的大肠杆菌含有任何一种常见的IS，每一种都有不到10份拷贝。当描述插入特定位置的IS如插入入噬菌体时，可以记为λ：：IS。IS都能编码自身转座所需的酶。多数IS元件在宿主DNA内有多个插入位点，但也有些在不同程度上偏爱的特定热点。 我们使用digIS^[18]分析工具,对细菌基因组做分析,得到IS插入序列预测结果。

---

3.5.4 重复序列预测

重复序列是基因组中不同位置出现的相同或互补性片段,是基因调控网络的组成成分。 根据重复的序列在基因组上的分布,分为散在重复序列、串联重复序列。
串联重复序列(Tandem Repeat,TR),即相邻的、重复两次或多次特定核酸序列模式的重复序列,分为微卫星(MicroSatellites),1-6个碱基为一个重复单元的简单重复序列(simple sequence repeats),以及10-60个碱基的长序列为一个重复单元的小卫星重复序列(MiniSatellites)。串联重复单元具有种属组成特异性,可作为物种的遗传性状,做进化关系的研究。
散在重复序列主要是转座子序列(transposable elements, TEs)。按长度分为短分散重复序列(Short Interspersed Nuclear Elements,SINEs)以及长散在重复序列(Long Interspersed Nuclear Elements,LINEs),其中长散在重复序列常具有转座活性。

我们使用trf(Tandem Repeats Finder)^[13]分析工具,对unicycler的组装结果assembly.fasta文件做分析,得到串联重复序列分析结果。使用RepeatMasker^[14]分析工具,得到散在重复序列分析结果。

---

3.5.5 基因岛预测

基因岛（Genomics Islands,GIs）是一些细菌、噬菌体或质粒中通过基因水平转移整合到微生物基因组中的一个基因组区段。一个基因岛能与病原机理、生物体的适应性等多种生物功能相关。通过比较基因组分析手段可研究具有特殊功能的微生物功能的特异性和功能来源。基于序列组成,采用IslandPath-DIOMB软件预测基因岛,通过检测序列中二核苷酸偏向性（phylogentically bias）和移动性基因（mobility genes,如转座酶或整合酶）以判定基因岛以及潜在的水平基因转移。

我们使用IslandPath-DIMOB^[17]分析工具,对细菌基因组prokka注释结果中的.gbk格式文件做分析,得到基因岛预测结果,然后用R语言绘制基因岛染色体位置示意图。

---

3.6 基因组功能预测

---

3.6.1 功能注释

GO(Gene Onotology),是生物学领域公认的,在分子和细胞层面的英文描述词条的参考规范。比如一个蛋白具有某种功能,尽管这是一种具体的功能,但是不同的人可能会有不同的描述,此时如果大家都采用GO里面的规范词条去描述,那么就不会出现很多偏门或者杂乱的描述词汇。GO促进了人们对生物学知识的交流和理解。 GO数据库是基因本体联合会(Gene Onotology Consortium)所建立的数据库,旨在建立一个适用于各种物种的,对基因和蛋白质功能进行限定和描述的,并能随着研究不断深入而更新的语言词汇标准。
GO提供了一系列的词条（terms）,用于描绘基因（基因产物）的特点,这些词条分为3大类：
(1) 细胞学组件（cellular component）,用于描述亚细胞结构、位置和大分子复合物,例如外部封装结构（external encapsulating structure）等。
(2) 分子功能（molecular function）,用于描述基因（基因产物）的功能,比如蛋白质结合转录因子活性（protein binding transcription factor activity）。
(3) 生物学过程（biological process）,指的是分子功能的有序组合以实现更复杂的生物功能,例如树突状细胞的抗原处理和递呈（dendritic cell antigen processing and presentation）。

我们使用emapper^[29]注释工具,将细菌基因组prokka注释结果中的.faa文件的基因蛋白序列作为query查询序列,到eggnog蛋白数据库做比对搜索,获得GO注释信息,用R语言做分类统计。

KEGG注释

KEGG数据库于 1995 年由 Kanehisa Laboratories 推出 0.1 版,目前发展为一个综合性数据库,其中最核心的为 KEGG PATHWAY 和 KEGG ORTHOLOGY 数据库。在 KEGG ORTHOLOGY 数据库中,将行使相同功能的基因聚在一起,称为 Ortholog Groups (KO entries),每个 KO 包含多个基因信息,并在一至多个 pathway 中发挥作用。
在 KEGG PATHWAY 数据库中,将生物代谢通路划分为以下6类：
(1) 细胞过程（Cellular Processes）
(2) 环境信息处理（Environmental Information Processing）
(3) 遗传信息处理（Genetic Information Processing）
(4) 人类疾病（Human Diseases）
(5) 新陈代谢（Metabolism）
(6) 生物体系统（Organismal Systems）

我们使用emapper注释工具,将细菌基因组prokka注释结果中的.faa文件的基因蛋白序列作为query查询序列,到eggnog蛋白数据库做比对搜索,获得KEGG_Pathway_map_id注释信息,用R语言做分类统计。

CAZy注释

碳水化合物活性酶（CAZy）数据库,录入的是能降解、修饰或者生成糖苷键的酶的功能结构域（或称模块）的信息。

CAZy数据库^[31]收录了碳水化合物活性酶的两种常见模块的数据信息:
（1）具有催化活性的模块,分为5类
    GH,Glycoside Hydrolases,糖苷水解酶。
    GT,Glycosyl Transferases,糖基转移酶。
    PL,Polysaccharide Lyases,多糖裂解酶。
    CE,Carbohydrate Esterases,碳水化合物酯酶。
    AA,Auxiliary Activities,辅助活性模块（一般是氧化还原酶,跟其它的碳水化合物活性酶共同发生作用）。
（2）结合在催化活性模块之上的其它模块,现有1类
    CBM,Carbohydrate-Binding Modules,与碳水化合物发生结合作用的模块,一般只起到结合作用,而没有催化作用。

我们使用蛋白序列比对工具diamond,将prokka注释得到的.faa文件中的蛋白序列作为查询序列,到CAZy数据库做比对搜索,获得CAZy数据库注释信息,然后用R语言对注释结果做分类统计。

SWISS-PROT数据库注释

SWISS-PROT是经过注释的蛋白质序列数据库,由欧洲生物信息学研究所(EBI)维护。数据库由蛋白质序列条目构成,每个条目包含蛋白质序列、引用文献信息、分类学信息、注释等,注释中包括蛋白质的功能、转录后修饰、特殊位点和区域、二级结构、四级结构、与其它序列的相似性、序列残缺与疾病的关系、序列变异体和冲突等信息。SWISS-PROT中尽可能减少了冗余序列,并与其它30多个数据建立了交叉引用,其中包括核酸序列库、蛋白质序列库和蛋白质结构库等。

我们使用diamond蛋白序列比对工具,将细菌基因组prokka注释结果中的.faa文件的基因蛋白序列作为query查询序列,到Swiss-Prot蛋白数据库做比对搜索,得到细菌的Swiss-Prot蛋白注释信息,生成表格文件。
推荐用Excel或WPS软件查看*mergeInfo.txt表格文件。

TCDB数据库^[32]注释

TCDB（Transporter Classification DataBase）是对转运蛋白进行分类的一个数据库。类似于对酶进行分类的EC系统,TCDB对于每一个转运蛋白家族,提供了一个TC Nmuber, TC Number 由小数点分隔的5段数字或者字母构成。每一段的数字或字母代表某一个层级的分类,第一级分类包括7个大类。目前TCDB提供了超过800个转运蛋白家族, 包含10000多条唯一的蛋白质序列和10000多篇文献。

我们使用diamond比对工具,将细菌基因组prokka注释结果中的.faa文件的基因蛋白序列作为query查询序列,到TCDB数据库做比对搜索,得到细菌的转运蛋白分类注释信息,然后用R语言做分类统计。

Pfam数据库^[30]注释

Pfam可以理解为是Protein family蛋白质家族的英文单词的缩写。该数据库主要提供蛋白质结构域家族的分类信息,被广泛用于查询蛋白质结构域注释信息及其多序列比对信息。在该数据库中,每个蛋白结构域家族由多序列比对和HMMs（hidden Markovmodels,隐马尔可夫模型）所组成。Pfam-A根据最新的UniProtKB蛋白序列数据库所构建而成,是人工注释和检查的蛋白结构域信息数据库,可信度较高。pfam_scan是Pfam官网提供的工具软件,用来分析蛋白序列具有哪些结构域。

我们使用pfam_scan分析工具,以Pfam-A数据库作为参考数据库,对细菌基因组prokka注释结果中的.faa文件的基因蛋白序列做注释,得到蛋白质结构域家族注释信息,生成表格文件。

---

3.6.2 致病系统分析

VFDB^[26]数据库注释

毒力因子指由细菌,病毒,真菌等产生的带有侵袭力和毒素等毒力性质的分子,主要用于微生物感染宿主时,通过抑制或逃避宿主的免疫反应等出入宿主组织细胞,并从宿主获得营养及自身增殖生长的目的。毒力因子可编码在可移动遗传元件（比如质粒、基因岛、噬菌体等）上并进行水平基因转移（传播）,使无害细菌变成危险的病原菌,所以在鉴定毒力因子时一般会考虑基因岛、分泌蛋白等。VFDB数据库由中国医学科学院研发,收集整理了24个属100多种重要医学病原菌已知毒力因子的组成、结构、功能、致病机理、毒力岛、序列和基因组信息等内容,被广泛应用于毒力因子基因鉴定。

我们使用diamond^[24]蛋白序列比对工具,将细菌基因组prokka注释结果中的.faa文件的基因蛋白序列作为query查询序列,到VFDB细菌毒力因子数据库做比对搜索,得到细菌的毒力因子注释信息。

提示：VFDB官网发布了网页在线分析工具VFanalyzer,推荐用户登录http://www.mgc.ac.cn/cgi-bin/VFs/v5/main.cgi 网站,上传prokka分析结果文件夹中的.gbk注释文件,提交在线分析任务,可以获得关于病原细菌毒力因子更丰富的分析结果。

CARD^[27]数据库注释

微生物抗生素耐药性问题,是在人类为抑制病原微生物生长繁殖,而高频次、大剂量使用抗生素的背景下凸显出来的。微生物通过自身基因突变,或者在环境中由于基因水平转移,得到的这些（突变的）基因,使药物作用的靶位发生变异或药物不能正常地发挥作用,从而获得对特定抗生素的抗性。
CARD(Comprehensive Antibiotic Resistance Database): 综合性抗生素抗性数据库,是目前使用最广泛的抗性基因数据库之一,目前包括约4000个抗性基因分类。

我们使用card_rgi分析工具,以CARD数据库作为参考数据库,对unicycler的组装结果assembly.fasta文件,做细菌的抗生素抗性基因注释,找出基因组上的抗性基因,生成注释结果表格文件*.CARD.txt,推荐使用Excel或WPS软件来查看。

PHI^[25]数据库注释
病原-宿主相互作用,是指微生物在侵染宿主的过程中,跟宿主会有一系列的相互作用,比如宿主需要探测到微生物的某些产物,宿主才能激活正常的免疫反应。那么在这个互作的过程中,微生物的某些基因就发挥着重要的作用。

PHI-base（Pathogen Host Interactions）,病原-宿主相互作用数据库。该数据库收录有,微生物中的、能影响病原-宿主相互作用的基因的信息。
其中,我们重点关注,PHI-base收集整理的,在已发表的研究论文中,人为地对微生物的影响病原-宿主相互作用的基因,做基因突变后,造成的微生物菌株的表型信息。

我们使用diamond蛋白序列比对工具,将细菌基因组prokka注释结果中的.faa文件的基因蛋白序列作为query查询序列,到PHI数据库做比对搜索,得到在病原-宿主相互作用中发挥作用的基因,在人为突变后的菌株表型的相关信息,用R语言做分类统计。

---

3.6.3 代谢系统分析

├── antiSMASH
│   │   ├── svg,    [次级代谢基因簇及相应基因的数量统计]
│   │   ├── *.antismash.tsv     [ 次级代谢基因簇及基因数量统计列表]
│   │   ├── *.cluster.stat.tsv   [ 次级代谢基因簇及基因数量统计]
│   │   ├── *.cluster.stat.png    [次级代谢基因簇及相应基因的数量统计图]
├── Secretory_Protein_dir
│   └── *
│       ├── *.faa
│       ├── *.log
│       ├── *.membrane.membrane.faa
│       ├── *.membrane.secretory.faa
│       ├── *.membrane.sigseq
│       ├── *.region_output.gff3
│       ├── *.signalp.output.gff3
│       ├── *.sigseq.pep
│       ├── *.tmhmm.txt
│       ├── done
│       ├── prediction_results.txt
│       ├── region_output.gff3
│       ├── Statistics.txt
├── T3SS_dir
│   └── *
│       ├── *.isnot_secretion.faa
│       ├── *.is_secretion.faa
│       ├── *.log
│       ├── *.secretion.sum.csv
│       ├── *.t3ss.txt
│       └── done
└── TNSS_dir
    └── *
        ├── *.tnss.sum.txt
        ├── *.tnss.reslut.txt

次级代谢物合成基因簇分析

"anti-biotics and Secondary Metabolite Analysis SHell — antiSMASH" 意为抗生素和次级代谢物的分析流程软件。AntiSMASH是用于寻找次级代谢物合成基因簇的最流行的工具之一。AntiSMASH数据库包含约6,200个完整细菌基因组和18,576个细菌草图基因组的注释。

我们使用antiSMASH^[23]分析工具,以antiSMASH数据库作为参考数据库,对unicycler的组装结果assembly.fasta文件做分析,得到次级代谢物合成基因簇分析结果。

---

蛋白信号肽和蛋白跨膜区预测

信号肽是蛋白质N-末端一段编码长度为5-30的疏水性氨基酸序列,用于引导新合成蛋白质向通路转移的短肽链。信号肽存在于分泌蛋白、跨膜蛋白和真核生物细胞器内的蛋白中。信号肽指引蛋白质转移的方式有两种：（1）常规的分泌（Sec/secretory）通路;（2）双精氨酸转移（Tat/twin-arginine）通路。前者存在于原核生物蛋白质转移到质膜过程中,以及真核生物蛋白质转移到内质网膜的过程中。后者存在于细菌、古菌、叶绿体和线粒体中,信号肽序列较长、疏水性较弱且尾部区含有两个连续精氨酸。相比于前者转运非折叠蛋白质,后者能转运折叠蛋白质跨越双层脂质膜。

信号肽指引蛋白质转运后,将由信号肽酶进行切除。信号肽酶有三种：（1）一型信号肽酶（SPaseI）;（2）二型信号肽酶（SPaseII）;（3）三型信号肽酶（SPaseIII）。大部分信号肽由SPaseI进行移除,SPaseI存在古菌、细菌和真核生物中,且在真核生物的内质网膜上仅存在一型信号肽酶。细菌和古菌脂蛋白的信号肽C端含有一段称为 lipobox 的保守区域,由SPaseII切除其信号肽,且lipobox紧邻切除位点（CS/Cleavage Site）的氨基酸是半胱氨酸,这和锚定到膜的功能是相关的。细菌的四型菌毛蛋白信号肽由SPaseIII进行切除。此外：分泌通路（Sec）相关信号肽能由SPaseI、SPaseII和SPaseIII切除,但是双精氨酸转移（Tat）通路相关信号肽仅由 SPaseI和SPaseII切除。

SignalP^[19]是最为常用的信号肽分析软件,常用于分泌蛋白预测的第一步。SignalP版本更新到了第六版,结合深度神经网络（deep neural network）、条件随机场分类（Conditional random field classification）和迁移学习（transfer learning）方法,能对信号肽进行更准确的预测。

我们使用SignalP6.0分析工具,对细菌基因组prokka注释结果中的.faa文件的基因蛋白序列做分析,得到蛋白的信号肽预测结果。
蛋白跨膜区即蛋白质序列中跨越细胞膜的区域,通常为α-螺旋结构,约20~25个氨基酸残基,构成跨膜区蛋白的氨基酸大部分是疏水性氨基酸。蛋白质含有跨膜区,提示它可能作为膜受体起作用,也可能是定位在定位在膜上的膜锚定蛋白或离子通道蛋白。TMHMM是最常见的跨膜区分析软件之一,它基于马尔可夫模型,综合跨膜区疏水性、电荷偏倚、螺旋长度和膜蛋白拓扑学等性质,可对跨膜区及膜内外取进行整体预测。

我们使用TMHMM^[20]分析工具,对细菌基因组prokka注释结果中的.faa文件的基因蛋白序列做分析,得到蛋白的跨膜区预测结果。

---

TNSS预测

病原菌通过分泌系统TNSS（type N secretion systems）将该类蛋白分泌至胞外或是宿主细胞中，通过控制免疫应答反应以及细胞衰亡引起病理反应，而其中革兰氏阴性菌的T3SS通常用来从分子水平研究病原菌，感染机制，毒力作用等，是研究得比较多的分泌系统。 对于TNSS系统，通过蛋白序列功能数据库注释结果中，提取分泌系统相关蛋白进行注释。

---

T3SS预测

III 型分泌系统（Type III secretion system，T3SS）主要是革兰氏阴性菌的分泌蛋白分泌到细胞外的运输途径，因此 III 型分泌系统效应蛋白（Type III secretion system Effector protein）与革兰氏阴性致病菌致病机理有关。 使用软件 EffectiveT3^[21] 对输入的氨基酸序列进行预测，通过其内部特定的计算模型对每条氨基酸序列进行评分，分值越高，可信度越高，选出评分高于阈值的序列，认为这些序列为 III 型分泌系统效应蛋白。

---

3.6.4 甲基化分析

使用SMRT Link^[35]最终的基因组组装结果进行甲基化位点检测和可能的甲基化转移酶识别的核苷酸基序（motif）的预测。能够预测到的修饰类型包括m6A，m4C和m5C，以及未知类型(modified_base)。 结果统计如下:

---

3.7 基因组可视化

--07.Genome_Drawing/
--*.pdf/    [Circos绘图pdf格式]
--*.svg/    [Circos绘图svg格式]
--*.png/    [Circos绘图结果图]

07.基因组可视化

Circos是由加拿大生物信息学科学家Martin Krzywinski利用Perl语言开发的一款可以用于描述关系型数据和多维数据的基因组圈图可视化软件。2009年Circos发表在Genome research。Circos不仅能将一个物种的整个基因组呈现在一张图片中,还可以给基因组添加丰富的注释信息,如功能注释信息,差异统计信息等等。

利用细菌基因组prokka^[10]蛋白预测信息,功能区注释信息（含正负链、CDS、RNA类型等信息）,CAZyme数据库注释结果（GH、GT、CE、PL、CBM、AA）,以及GC content、GC skew的统计结果绘制Circos圈图。非常规格式文件请用notepade++或记事本打开。

我们使用Circos^[34]软件,分别对细菌基因组的每一条染色体或contig作一个单独的完整圈图,标示基因正负链转录位置信息、RNA位置信息、碳水化合物酶位置信息、基因GC含量信息、基因GC-skew含量等信息。

这里展示的是对contig1绘制的circos圈图。

从外到内每一圈的含义说明如下：

第一圈：
细菌染色体;

第二圈：
红色高亮表示正链中的CDS区;
黑色高亮表示负链中的CDS区;

第三圈：
蓝色三角表示tRNA;
紫色三角表示rRNA;
灰色三角表示tmRNA;

第四圈：
青色方块表示GH;
绿色方块表示GT;
棕色方块表示CE;
粉红方块表示PL;
黄色方块表示CBM;
橙色方块表示AA; （方块较大,相邻区域会有重叠,但是完全可以直观的表现不同CAZy的数量差异）

第五圈;
直方图表示基因的GC content,即(G+C)/(A+G+C+T),高度表示GC content的值的大小,宽度表示基因的长度;
红色直方图,方向朝外表示该基因的GC content的值大于基因的总GC content的值;
蓝色直方图,方向朝内表示该基因的GC content的值小于基因的总GC content的值;

第六圈： 直方图表示基因的GC skew(偏移),GC偏移=(G-C)/(G+C)这个式子用来衡量G和C的相对含量,G>C则GC skew的值为正值,G<C则为负值。在大多数细菌基因组中,前导链(leading strand)和滞后链(lagging strand)在碱基组成上有明显的不同——前导链富含G和T,但是滞后链中的A和C则更多一些。打破A=T和C=G的碱基频率发生的偏移,被称之为“AT(AT-skew)”和“GC(GC-skew)”。由于通常GC偏移比AT偏移发生的更明显,所以习惯上更多地只考虑GC偏移。 橙色直方图表示contig的GC skew>0,高度表示绝对值大小; 绿色的直方图表示contig的GC skew<0,高度表示绝对值大小;

---

3.8 基因组比较分析

---

3.8.1 基因共线性分析

共线性分析

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

使用MUMmer软件对样本基因组和参考基因组进行比对，确定样本基因组和参考基因组之间的比对关系，识别样本基因组和参考基因组之间的共线性区域。样品基因组和参考基因组之间的共线性展示结果如下图

---

3.8.2 SNP分析

SNP和indel分析

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

使用snippy软件对样本基因组和参考基因组进行比对，检测之间的snp位点和indel位点.

全基因组结构变异类型统计如下所示:

---

3.8.3 SV分析

SV分析

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

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

全基因组结构变异类型统计如下所示:

---

3.8.4 泛基因组分析

泛基因组分析

泛基因组是某一物种全部基因的总称。核心基因（core genes）和可变基因(variable genes)是泛基因组研究重点。从结构上来看，核心基因是在所有动植物个体或者菌株中都同时存在的基因。可变基因：是指，在一个或者一个以上的动植物个体或者菌株中存在的基因，它们不是固定的，是多变的。如果某个基因，仅存在某一个动植物个体或者菌株中，该基因还可以细分为品系或者菌株特有基因。 Roary是一款由The Wellcome Trust SangerInstitute主导开发的并且可以快速对大规模原核生物基因组进行泛基因组分析的软件，其优点在于可以在短时间内对数千规模的微生物基因组进行泛基因组分析。

我们使用roary软件对样本的gff文件做泛基因组分析。

四 分析所用软件的版本

软件	版本
bam2fastq	1.3.1
fastp	0.23.1
quast	2.1.5
unicycler	2.2.1
bwa	0.7.17-r1188
Jellyfish	1
GenomeScope	V1
prokka	1.14.6
CompareM	2.2.1
pyani	2.2.1
Tandem Repeat Finder	2.1.5
RepeatMasker	v0.12.6
CRISPRCasFinder	v0.12.6
PhiSpy	v0.12.6
IslandPath-DIMOB	v0.12.6
digIS	v1.2
SignalP	v6.0
tmhmm	v2.0
EffectiveT3	2.0
MacSyFinder	2.0
antiSMASH	v6.0
diamond	v2.0.14
PHI	2022
VFDB	2022
CARD	2020
Diamond	v2.0.14
eggNOGmapper	104.3
pfam	35.0
CAZy	2022
tcdb	2022
emapper	2.0.1
Circos	0.69
SMRT Link	v11

---

五 参考文献

[1] https://github.com/pacificbiosciences/bam2fastx/

[2] Shifu Chen, Yanqing Zhou, Yaru Chen, Jia Gu; fastp: an ultra-fast all-in-one FASTQ preprocessor, Bioinformatics, Volume 34, Issue 17, 1 September 2018, Pages i884–i890, https://doi.org/10.1093/bioinformatics/bty560

[3] Fukasawa, Yoshinori et al. “LongQC: A Quality Control Tool for Third Generation Sequencing Long Read Data.” G3 (Bethesda, Md.) vol. 10,4 1193-1196. 9 Apr. 2020, doi:10.1534/g3.119.400864

[4] Gurevich, Alexey et al. “QUAST: quality assessment tool for genome assemblies.” Bioinformatics (Oxford, England) vol. 29,8 (2013): 1072-5. doi:10.1093/bioinformatics/btt086

[5] Wick, Ryan R et al. “Unicycler: Resolving bacterial genome assemblies from short and long sequencing reads.” PLoS computational biology vol. 13,6 e1005595. 8 Jun. 2017, doi:10.1371/journal.pcbi.1005595

[6] van der Graaf-van Bloois, Linda et al. “RFPlasmid: predicting plasmid sequences from short-read assembly data using machine learning.” Microbial genomics vol. 7,11 (2021): 000683. doi:10.1099/mgen.0.000683

[7] Li, Heng, and Richard Durbin. “Fast and accurate long-read alignment with Burrows-Wheeler transform.” Bioinformatics (Oxford, England) vol. 26,5 (2010): 589-95. doi:10.1093/bioinformatics/btp698

[8] https://github.com/jamesturk/jellyfish

[9] Vurture, Gregory W et al. “GenomeScope: fast reference-free genome profiling from short reads.” Bioinformatics (Oxford, England) vol. 33,14 (2017): 2202-2204. doi:10.1093/bioinformatics/btx153

[10] Seemann, Torsten. “Prokka: rapid prokaryotic genome annotation.” Bioinformatics (Oxford, England) vol. 30,14 (2014): 2068-9. doi:10.1093/bioinformatics/btu153

[11] https://github.com/dparks1134/CompareM

[12] Pritchard, L. et al. “Genomics and taxonomy in diagnostics for food security: soft-rotting enterobacterial plant pathogens.” Analytical Methods 8 (2016): 12-24.

[13] Benson, G. “Tandem repeats finder: a program to analyze DNA sequences.” Nucleic acids research vol. 27,2 (1999): 573-80. doi:10.1093/nar/27.2.573

[14] Tarailo-Graovac, Maja, and Nansheng Chen. “Using RepeatMasker to identify repetitive elements in genomic sequences.” Current protocols in bioinformatics vol. Chapter 4 (2009): Unit 4.10. doi:10.1002/0471250953.bi0410s25

[15] Couvin, David et al. “CRISPRCasFinder, an update of CRISRFinder, includes a portable version, enhanced performance and integrates search for Cas proteins.” Nucleic acids research vol. 46,W1 (2018): W246-W251. doi:10.1093/nar/gky425

[16] PhiSpy: a novel algorithm for finding prophages in bacterial genomes that combines similarity- and composition-based strategies.” Nucleic acids research vol. 40,16 (2012): e126. doi:10.1093/nar/gks406

[17] Bertelli, Claire, and Fiona S L Brinkman. “Improved genomic island predictions with IslandPath-DIMOB.” Bioinformatics (Oxford, England) vol. 34,13 (2018): 2161-2167. doi:10.1093/bioinformatics/bty095

[18] https://github.com/janka2012/digIS

[19] Teufel, Felix et al. “SignalP 6.0 predicts all five types of signal peptides using protein language models.” Nature biotechnology vol. 40,7 (2022): 1023-1025. doi:10.1038/s41587-021-01156-3

[20] Krogh, A., Larsson, B., Von Heijne, G., & Sonnhammer, E. L. (2001). Predicting transmembrane protein topology with a hidden Markov model: application to complete genomes. Journal of molecular biology, 305(3), 567-580. https://www.sciencedirect.com/science/article/pii/S0022283600943158

[21] Arnold, Roland et al. “Sequence-based prediction of type III secreted proteins.” PLoS pathogens vol. 5,4 (2009): e1000376. doi:10.1371/journal.ppat.1000376

[22] Abby, Sophie S et al. “MacSyFinder: a program to mine genomes for molecular systems with an application to CRISPR-Cas systems.” PloS one vol. 9,10 e110726. 17 Oct. 2014, doi:10.1371/journal.pone.0110726

[23] Blin, Kai et al. “antiSMASH 6.0: improving cluster detection and comparison capabilities.” Nucleic acids research vol. 49,W1 (2021): W29-W35. doi:10.1093/nar/gkab335

[24] Buchfink, Benjamin et al. “Fast and sensitive protein alignment using DIAMOND.” Nature methods vol. 12,1 (2015): 59-60. doi:10.1038/nmeth.3176.

[25] Urban, Martin et al. “PHI-base in 2022: a multi-species phenotype database for Pathogen-Host Interactions.” Nucleic acids research vol. 50,D1 (2022): D837-D847. doi:10.1093/nar/gkab1037

[26] Liu, Bo et al. “VFDB 2022: a general classification scheme for bacterial virulence factors.” Nucleic acids research vol. 50,D1 (2022): D912-D917. doi:10.1093/nar/gkab1107

[27] Alcock, Brian P et al. “CARD 2020: antibiotic resistome surveillance with the comprehensive antibiotic resistance database.” Nucleic acids research vol. 48,D1 (2020): D517-D525. doi:10.1093/nar/gkz935

[28] Buchfink, Benjamin et al. “Fast and sensitive protein alignment using DIAMOND.” Nature methods vol. 12,1 (2015): 59-60. doi:10.1038/nmeth.3176.

[29] Cantalapiedra, Carlos P et al. “eggNOG-mapper v2: Functional Annotation, Orthology Assignments, and Domain Prediction at the Metagenomic Scale.” Molecular biology and evolution vol. 38,12 (2021): 5825-5829. doi:10.1093/molbev/msab293

[30] Mistry, Jaina et al. “Pfam: The protein families database in 2021.” Nucleic acids research vol. 49,D1 (2021): D412-D419. doi:10.1093/nar/gkaa913

[31] Elodie Drula, Marie-Line Garron, Suzan Dogan, Vincent Lombard, Bernard Henrissat, Nicolas Terrapon, The carbohydrate-active enzyme database: functions and literature, Nucleic Acids Research, Volume 50, Issue D1, 7 January 2022, Pages D571–D577,

[32] Saier MH, Reddy VS,el al (2021). The Transporter Classification Database (TCDB): 2021 update. Nucleic Acids Res. 49(D1):D461-7

[33] Jaime, H. C. , Kristoffer, F. , Pedro, C. L. , Damian, S. , Juhl, J. L. , & Christian, V. M. , et al. (2016). Fast genome-wide functional annotation through orthology assignment by eggnog-mapper. Molecular Biology & Evolution(8), 2115.

[34] Circos: An information aesthetic for comparative genomics. Genome Research 2009

[35] Reiner J, Pisani L, Qiao W, et al. Cytogenomic identification and long-read single molecule real-time (SMRT) sequencing of aBardet–Biedl Syndrome 9(BBS9) deletion:[J]. Npj Genomic Medicine, 2018, 3(1).