内生镰刀菌基因组染色体级别组装和注释
Chromosome-Scale Genome Assembly and Annotation Method of Endophyte Fusarium    

材料与试剂

1. 内生镰刀菌Fusarium culmorum Class2-1B、Fusarium pseudograminearum Class2-1C，分离自沿海滩涂植物滨麦Leymus mollis，与植物共生可以促进植物生长和提高植物耐盐性 (Rodriguez et al.，2008; Redman et al.，2011; Pan et al.，2018)

仪器设备

1. 三代测序仪 (Pacific Biosciences PacBio RS II)
2. 二代测序仪 (Illumina HiSeq 2500)
   

软件和数据库

1. MECAT2 (https://github.com/xiaochuanle/MECAT2)
2. BUSCO v2.0 (https://busco.ezlab.org/)
3. tRNAscan-SE (http://lowelab.ucsc.edu/tRNAscan-SE)
4. RepeatModeler: http://www.repeatmasker.org/RepeatModeler
5. RepeatMasker: http://repeatmasker.org
6. NR (https://www.ncbi.nlm.nih.gov/refseq/about/nonredundantproteins)
7. Swiss-Prot (https://www.uniprot.org/statistics/Swiss-Prot)
8. KEGG databases (https://www.genome.jp/kegg/kegg1.html)
9. Repbase database: https://www.girinst.org/server/RepBase
10. Fungi odb10 dataset: https://busco.ezlab.org/frames/fungi.htm
11. TRF (Tandem repeats finder) http://tandem.bu.edu/trf/trf.unix.help.html
12. LTR_FINDER http://tlife.fudan.edu.cn/tlife/ltr_finder
13. Augustus http://bioinf.uni-greifswald.de/augustus/
14. GlimmerHMM http://ccb.jhu.edu/software/glimmerhmm/
15. Piler http://www.drive5.com/piler
16. RepeatScout https://github.com/mmcco/RepeatScout
17. TrEMBL https://www.uniprot.org/statistics/TrEMBL
18. Interpro https://www.ebi.ac.uk/interpro/
19. Maker http://www.yandell-lab.org/software/maker.html
    
20. Fusarium culmorum strain PV, whole genome shotgun sequencing project https://www.ncbi.nlm.nih.gov/nuccore/PVEM00000000
21. Fusarium pseudograminearum CS3096, whole genome shotgun sequencing project https://www.ncbi.nlm.nih.gov/nuccore/AFNW00000000
    

实验步骤

一、测序

1. 使用太平洋生物科学公司开发的单分子实时 (SMRT) 测序和Illumina HiSeq 2500测序技术来组装完整的基因组。测序在北京诺禾致源生物信息技术有限公司进行。
2. 取单孢分离后培养15天的内生镰刀菌Fusarium culmorum、Fusarium pseudograminearum PDA平板，使用Omega真菌DNA提取试剂盒提取DNA，DNA浓度大于100 ng/μl，DNA纯度 (OD_260/280 在1.8-2.0 之间；OD_260/230 在2.0-2.2 之间) ，使用50 mg DNA构建PacBio和Illumina测序文库。
3. 对PacBio文库，构建每个菌株的20 kb插入片段大小的标准SMRTbell文库，用PacBio Sequel II系统对PacBio长读序列进行测序。
4. 为了完善基于PacBio long-read的基因组组装，在Illumina HiSeq 2500上对插入大小为500 bp的双端Illumina DNA文库进行了测序。
5. 基于Illumina Short Reads的数据，分析了两株内生镰刀菌基因组的K-mer分布，并估计了两株内生镰刀菌基因组的大小。
   
   
   图 1. Illumina和PacBio测序流程图
   
   图1展示了第二代测序Illumina和第三代测序PacBio技术的测序流程，结合二代和三代测序数据进行了高质量的基因组组装。
   
   
   图 2. 真菌染色体结构模式图
   
   图1模式图展示了真菌的染色体两端具有端粒结构，在基因组组装中，染色体端粒到端粒的组装代表染色体的完整性，也是高质量基因组组装结果的标志。
   

二、基因组组装和注释

1. 获得了16.7 GB的long-read数据F. culmorum Class2-1B，其中Scaffold N50的长度为9.63 M；而F. pseudograminearum Class2-1C，获得了19.7 GB的long-read数据，Scaffold N50的长度为9.15 M。
2. 基于PacBio测序数据，我们使用MECAT2来进行基因组的组装.MECAT2是一个超快速的准确比对，误差校正、组装工具。利用MECAT2进行了基因组组装和纠错，读取的文件格式为FASTA，基因组大小设置为40 Mb。MECAT2首先利用reads相互比对进行纠错，生成consensus序列。然后再把纠错后的最长的30X reads与低质量的reads进行比对，利用overlap对低质量序列进行纠正。最后，根据reads间的高质量重叠关系进行组装。然后使用Pilon (v1.22) 进行两轮纠错，接下来，我们对以上的组装结果再进行Polish，消除其中的Indel错误。因为三代Pacbio 数据有很高的错误率，所以在使用三代 Pacbio 数据完成组装之后，依然存在少量的Indel，需要将二代测序clean数据用bwa比对到mecat的组装结果上，结合比对结果文件，用pilon纠错。最后我们将二代小片段数据再比对回polish后的序列，用pilon再次对组装结果进行纠错。
3. MECAT2的软件参数，在config配置文件中，如下
   PROJECT=F. culmorum #项目名称
   RAWREADS=m54220.fasta #处理的原始reads（FASTA格式）
   GENOME_SIZE=40000000 #基因组大小为40000000 bp。
   THREADS=20 #使用20个CPU线程
   MIN_READ_LENGTH=2000 #用于纠错和修建的reads的最低长度为2000 bp。
   CNS_OVLP_OPTIONS="" #在纠错阶段是检测候选overlap的参数, 会传给mecat2pw
   CNS_OPTIONS="-r 0.6 -a 1000 -c 4 -l 2000" #原始reads纠错参数，传递给mecat2cns
   TRIM_OVLP_OPTIONS="-B" #在trim阶段，用于检测重叠的参数,会传给v2asmpm
   ASM_OVLP_OPTIONS="-n 100 -z 10 -b 2000 -e 0.5 -j 1 -u 0 -a 400" #在组装阶段，用于检测overlap的参数，传给v2asmpm.sh
   FSA_OL_FILTER_OPTIONS="--max_overhang=-1 --min_identity=-1" #过滤overlap的参数，传递给fsa_ol_filter
   FSA_ASSEMBLE_OPTIONS="" #组装trimm reads的参数, 传给v2asmpm
   USE_GRID=false #无多个计算节点
   CLEANUP=0 #运行结束后不删除MECAT2的中间文件。
   CNS_OUTPUT_COVERAGE=30 #选择30倍的最长纠错覆盖后，将reads进行trim和组装。
   GRID_NODE=0 #当USE_GRID=1时，设置用到的计算节点数，单节点服务器不需要设置。
   
   Pilon (v1.22) 主要参数：
   --changes 指定生成一个列出output.fasta中更改的文件
   --vcf #指定生成vcf文件
   --tracks #vcf文件将包含QE(质量加权)数据
   --diploid #样本来自二倍体有机体；最终会影响杂合子SNPs的检测。
   --fix bases #尝试解决的问题类别用逗号分隔。
4. F. culmorum Class2-1B和F. pseudograminearum Class2-1C都分别得到了6和7个Scaffold，参照两株镰刀菌的参考基因组：F. culmorum PV，F. pseudograminearum CS3096 (Schmidt et al., 2018; Gardiner et al., 2012)。端粒是真核生物染色体末端的DNA重复序列，作用是保持染色体的完整性和控制细胞分裂周期。将F. culmorum组装成的四条染色体，其中两条是端粒对端粒，而另外两条只在一端有一个已识别的端粒。将F. pseudograminearum组装成的四条染色体，其中三条两端都有端粒结构，而另外一条只在一端有一个已识别的端粒。在Scaffold末端发现了TTAGGG的串联重复序列 (或互补DNA链序列，AATCCC)。F. culmorum和F. pseudograminearum的Scaffold至少有一端含有端粒结构，每条Scaffold都接近完整染色体的长度 (Aksenova and Mirkin, 2019)。如上所示， 两株内生镰刀菌基因组的染色体都含有图2中的端粒结构，说明我们得到了两株高质量组装的内生镰刀菌基因组。通过BLAST搜索鉴定了F. culmorum中的2个短Scaffold为线粒体基因组，总GC含量为31.2%。同样，通过BLAST搜索鉴定了F. pseudograminearum中的3个短Scaffold为线粒体基因组，总GC含量为34.6%，进一步分别比较它们的同种镰刀菌线粒体基因组时，发现这两株内生镰刀菌线粒体基因组都显示出大于98%的序列同源性 (Kulik et al., 2020)。
   
5. 通过结合de novo注释和基于同源的预测结果进行蛋白质编码基因的结构注释 (Rigden, 2017) 。然后在MAKER软件的帮助下，将上诉两种方法预测得到的基因集整合成一个非冗余的、更加完整的基因集，同时通过手动整合得到最终的可靠基因集。然后借助于外源蛋白数据库（SwissProt、TrEMBL、KEGG、InterPro等）对基因集中的蛋白进行功能注释。 
   Maker（v2.31.9）基因预测使用的主要参数：
   genome=F. culmorum.FASTA #基因组序列((FASTA文件或嵌入GFF3文件的FASTA))
   organism_type=eukaryotic #真核生物
   snaphmm=#SNAP HMM 文件
   gmhmm= #GeneMark HMM 文件
   augustus_species= #Augustus基因预测物种模型
   fgenesh_par_file= #FGENESH参数文件
   pred_gff= #GFF3的从头计算预测文件
   est2genome=0 #禁止从EST进行基因预测
   （Maker利用现有的软件工具（其中一些是基因预测）并对其输出进行集成，以根据证据比对得出给定位置的Maker认为是最佳可能的基因模型，上面用到的软件参数均为默认参数。）
   
6. 使用BUSCO (Benchmark Universal Single-Copy Orologs) Fungi odb10数据库 (v.4.0.6) 对基因注释和基因组组装质量进行评估，结果显示Class2-1B和Class2-1C的基因注释和基因组组装质量分别为98.8%和99.1% (总共搜索了758个保守核心蛋白) ，这表明俩个基因组的组装质量是非常高的 (Simão et al., 2015)。使用Busco软件评估时使用的为默认参数，仅依据是蛋白质的评估/基因组评估调整“-mode =prot/genome #蛋白质或者基因组”。使用的busco库为fungi odb10。
7. 重复序列是基因组的重要组成部分，主要包括两大类：分别为串联重复序列（Tandem repeat）和散在重复序列(Interpersed repeat)。其中串联重复序列包括有微卫星序列，小卫星序列等等；散在重复序列又称转座子元件，包括以DNA-DNA方式转座的DNA转座子和反转录转座子(retrotransposon)。对于转座子 (TEs) 注释，RepeatMasker (v.4.07) 用于Repbase数据库 (v.23.06) (Bao et al.，2015) 来识别已知的TEs。同时，还使用RepeatModeler (v1.0.11) 和LTR finder (Jurka et al，2005) 进行从头检测。在Class2-1B和Class2-1C中分别鉴定出约1.55Mb和2.04Mb的TEs (占总基因组的4.13%和5.37%) 。
   RepeatMasker主要参数：
   -nolow #不会掩盖低复杂度的DNA或简单的重复序列
   -no_is #跳过细菌插入元素检验
   –norna #不掩盖小RNA（伪）基因
   -parallel 1 #并行使用1个处理器（仅适用于批处理超过50 kb的文件或序列）
   
   RepeatModeler主要参数：默认。
   
   LTR finder主要参数:
   -w 2 #输出格式为表格
   -s LTR_FINDER.x86_64-1.0.5/tRNAdb/Sc-tRNAs.fa#使用指定RNA数据库(名称)
   
   

结果分析

一、组装结果

1. 基因组评估
   在基因组组装前，为了用测序所得的read信息估计基因组特征，我们使用二代测序数据基于K-mer的分析方法来估计基因组大小和杂合率等基本信息。
   在本项目的分析当中，我们取K为17来进行分析，两个物种分析结果如下：
   
   
   图3 Fusarium culmorum Kmer分布图
   
   
   图4 Fusarium pseudograminearum Kmer分布图
   
   注：图中横坐标为深度（depth），纵坐标为各深度下的K-mer种类占所有K-mer种类的比例。
   
   从以上的结果可以看到两个物种的Kmer分布图均只有一个峰，说明这两个物种的杂合度都很低。根据kemer分析结果，预估的基因组大小: F. culmorum为40,053,289 bp，F. pseudograminearum为42,895,120 bp。
   
2. 基因组特征 
   两株内生镰刀菌（F. culmorum和F.pseudograminearum）的基因组测序后分别得到6和7个Scaffold，组装的总基因组大小为40.05 Mb和42.90 Mb，其中黄色镰刀菌被组装成4条染色体，2条为端粒对端粒，2条为有一端为端粒结构；假禾谷镰刀菌也被组装成4条染色体，3条为端粒对端粒，1条为有一端为端粒结构；N50 Scaffold长度分别为9.63 Mb和9.15 Mb；GC含量分别为47.4和47%；转座子大小分别为1.5Mb和2.04 Mb；黄色镰刀菌中的2个短Scaffold为线粒体基因组，大小为136406 bp；假禾谷镰刀菌中的3个短Scaffold为线粒体基因组，大小为136045 bp。
   
   表1. 黄色镰刀菌和假禾谷镰刀菌的基因组特点和预测特征
   
   
   

二、注释结果

通过结合de novo注释和基于同源的预测结果进行蛋白质编码基因的结构注释。使用Maker (v.2.31.9) 分别在F. culmorum和F. pseudograminearum中预测了11,450和11,221个完整的蛋白编码基因模型。发现F. culmorum和F. pseudograminearum中分别有97.06%和96.93%的基因可以在InterProScan、Gene Ontology、KEGG以及NR数据库被注释（见表2）。在F. culmorum和F. pseudograminearum中分别鉴定出总的重复序列为1572977 bp和2060381 bp，占总基因组的4.18%和5.43%。其中鉴定出TEs约为1.55 Mb和2.04 Mb (占总基因组的4.13%和5.38%)；串联重复分别鉴定出1.95 Mb和2.33 Mb，占总基因组的5%和6%（见表3和表4）。




表2 F. culmorum和F. pseudograminearum基因组功能注释统计结果


表3 F. culmorum重复序列分类结果统计


表4 F. pseudograminearum重复序列分类结果统计




致谢

本protocol的研究工作得到课题“内生镰刀菌促进树木生长和耐盐性的分子调控机制研究”资助经费，课题编号为76B2018001。

参考文献

1. Rodriguez, R. J., Henson, J., Van Volkenburgh, E., Hoy, M., Wright, L., Beckwith, F., Kim, Y. O. and Redman, R. S. (2008). Stress tolerance in plants via habitat-adapted symbiosis. ISME J 2: 404-416.
2. Redman, R. S., Kim, Y. O., Woodward, C. J. D. A., Greer, C., Espino, L., Doty, S. L. and Rodriguez, R. J. (2011). Increased fitness of rice plants to abiotic stress via habitat adapted symbiosis: A strategy for mitigating impacts of climate change. PLoS One 6: 1-10.
3. Pan, X. Y., Sun, H. J. and Yuan, Z. L. (2018). Toxin accumulation of three Leymus mollis-associated endophytic Fusarium Isolates and their effects 200 on growth and salt tolerance of Liquidambar styraciflua seedlings. For Res 31: 64–73.
4. Schmidt, R., Durling, M. B., de Jager, V., Menezes, R. C., Nordkvist, E., Svatoš, A., Dubey, M., Lauterbach, L., Dickschat, J. S., Karlsson, M. and Garbeva, P. (2018). Deciphering the genome and secondary metabolome of the plant pathogen fusarium culmorum. FEMS microbiology ecology 94(6): fiy078.
5. Gardiner, D. M., McDonald, M. C., Covarelli, L., Solomon, P. S., Rusu, A. G., Marshall, M., Kazan, K., Chakraborty, S., McDonald, B. A. and Manners, J. M. (2012). Comparative pathogenomics reveals horizontally acquired novel virulence genes in fungi infecting cereal hosts. PLoS Pathog 8(9): e1002952.
6. Aksenova, A. Y. and Mirkin, S. M. (2019). At the beginning of the end and in the middle of the beginning: structure and maintenance of telomeric dna repeats and interstitial telomeric sequences. Genes (Basel) 10: 118.
7. Kulik, T., Brankovics, B., van Diepeningen, A. D., Bilska, K., Żelechowski, M., Myszczyński, K., Molcan, T., Stakheev, A., Stenglein, S, Beyer, M., Pasquali, M., Sawicki, J., Wyrȩbek, J. and Baturo-Cieśniewska, A. (2020).Diversity of mobile genetic elements in the mitogenomes of closely related Fusarium culmorum and F. graminearum sensu stricto strains and its implication for diagnostic purposes. Front Microbiol 11: 1002.
8. Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. and Zdobnov, E. M. (2015). BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31: 3210–3212.
9. Bao, W., Kojima, K. K. and Kohany, O. (2015). Repbase Update, a database of repetitive elements in eukaryotic genomes. Mob DNA 6: 4–9.
10. Rigden, D. J. (2017). From protein structure to function with bioinformatics: second edition.
11. Jurka, J., Kapitonov, V. V., Pavlicek, A., Klonowski, P., Kohany, O., Walichiewicz, J. (2005). Repbase Update, a database of eukaryotic repetitive elements. Cytogentic and Genome Research 110: 462-467.
12. Benson, G. (1999). Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acids Res 27: 573-580.
13. Price, A. L., Jones, N. C. and Pevzner, P. A. (2005). De novo identification of repeat families in large genomes. Bioinformatics 21: i351-i358.
14. Edgar, R. C. and Myers, E. W. Piler: (2005). Identification and Classification of genomic repeats. Bioinformatics 21: i152-158.
15. Xu, Z. and Wang, H. Ltr. (2007). Finder: an efficient tool for the prediction of full-length ltr retrotransposons. Nucleic Acids Res 35: W265-268.
16. Kent, W. J. (2002). BLAT-the BLAST-like alignment tool. Genome Res 12: 656–664.
17. Guy, S. and Ewan, B. (2005). Automated generation of heuristics for biological sequence comparison. BMC Bioinformatics 6: 31
18. Stanke, M., Keller, O., Gunduz, I., Hayes, A., Waack, S. and Morgenstern, B. (2006). AUGUSTUS: ab initio prediction of alternative transcripts. Nucleic Acids Res 34: W435-W439.
19. Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., Salzberg, S. L., Wold, B. J. and Pachter, L. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol 28: 511-515.
20. Majoros, W. H., Pertea, M. and Salzberg, S. L. (2004). TigrScan and GlimmerHMM: two open. Bioinformatics 20(16): 2878-9.
21. Carson, H. and Mark, Y. (2011). MAKER2: an annotation pipeline and genome-database management tool for second-generation genome projects. BMC Bioinformatics 12: 491.
22. Bairoch, A. and Apweiler, R. (2000). The SWISS-PROT protein sequence database and its supplement TrEMBL in 2000. Nucleic Acids Res 28: 45-48.
23. Zdobnov, E. M. and Apweiler, R. (2001). InterProScan - an integration platform for the signature-recognition methods in InterPro. Bioinformatics 17: 847-848.
24. Ashburner, M. Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., Davis, A. P., Dolinski, K., Dwight, S. S., Eppig, J. T. Harris, M. A., Hill, D. P., Issel, Tarver, L., Kasarskis, A., Lewis, S., Matese, J. C., Richardson, J. E., Ringwald, M., Rubin, G. M. and Sherlock, G (2000). Gene Ontology: tool for the unification of biology. Nat Genet 25(1): 25-29.
25. Kanehisa, M. and Goto, S. (2000). KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res 28: 27-30.
26. Griffiths-Jones, S., Moxon, S., Marshall, M., Khanna, A., Eddy, S. R. and Bateman, A. (2005). Rfam: annotating non-coding RNAs in complete genomes. Nucleic Acids Res 33: D121-4.
27. Todd M. Lowe and Sean R. Eddy. (1997). tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 25(5): 966-64.