《A complete telomere-to-telomere assembly of the maize genome》方法总结

研究背景

完整的T2T基因组组装一直是基因组研究的长期追求。

研究方法

通过生成高深度覆盖的超长 Oxford Nanopore Technology (ONT) 和 PacBio HiFi 测序数据，报道了玉米的完整基因组组装。
每条染色体均以单一contig的形式完整覆盖。

结果概述

基因组特征：
- Mo17基因组的大小为 2,178.6 Mb，基础准确率超过 99.99%。
- 揭示了基因组所有重复区域的结构特征。
简单序列重复：
- 发现多个超长简单序列重复阵列，其中连续的 TAG 三核苷酸重复长度可达 235 kb。
核仁组织区域：
- 完整组装了 26.8 Mb 的核仁组织区域，包含 2,974 个 45S rDNA 拷贝，揭示了rDNA重复和转座子插入的复杂模式。
中心粒分析：
- 完成了十个中心粒的完整组装，精确剖析了CentC丰富和CentC贫乏中心粒的重复组成。

研究意义

Mo17基因组的完整组装标志着在理解高等植物基因组中难以处理的重复区域复杂性方面的重要进展。

1.植物材料与基因组测序

植物材料:
- 采用玉米（Zea mays）自交系Mo17，在温室中生长。
- 生长条件：光照下30°C 16小时，黑暗中25°C 8小时。
- 在生长14天后，采集新鲜幼叶组织，立即用液氮冷冻以进行DNA提取。
DNA提取:
- 使用CTAB法提取高分子量基因组DNA，随后使用Qiagen基因组试剂盒（Qiagen, 13343）进行纯化。
测序文库构建:
- 利用高分子量DNA构建PacBio HiFi和ONT（Oxford Nanopore Technologies）测序文库。
- 还采用核提取法构建了超长测序文库。
测序平台:
- 在Nanopore PromethION测序仪上进行ONT测序。
- 在PacBio Sequel II平台上进行PacBio测序。

2.全基因组组装

ONT数据组装:
- 使用428.7 Gb的高质量过滤超长ONT数据。
- 选择最长的332.1 Gb读取，N50为100.7 kb。
- 使用NextDenovo（v2.2-beta.0）进行组装，使用了特定参数。
PacBio HiFi数据组装:
- 利用151.1 Gb的PacBio HiFi数据。
- 同时使用Hifiasm（v.0.7）和Canu（v.2.0）进行组装。
- Hifiasm组装使用默认参数；Canu组装使用了定制参数，包括内存和计算资源选项。

3.全基因组抛光

抛光过程:
- 使用ONT数据、PacBio HiFi数据和Illumina无PCR数据对初始ONT组装的567个contig进行抛光，采用NextPolish（v1.1.0）。
- 先使用ONT数据抛光三轮，再用PacBio HiFi数据抛光三轮，最后用Illumina数据抛光四轮。
- 由于PacBio HiFi数据质量高，因此未对使用PacBio数据组装的contig进行进一步抛光。
- 手动提取闭合间隙的序列及其上下游区域，并进行迭代抛光。

4.伪分子的构建

伪分子构建:
- 基于使用玉米泛基因组基因型测序（GBS）标签的高密度遗传图构建伪分子。
- 将初始ONT组装的567个contig与基因组图进行比对，采用bwa mem工具。
- 从19个非冗余contig中构建伪分子，以覆盖其余的548个contig。
- 由于组装错误，手动将一个contig拆分为两个，根据GBS标签的位置确定contig的顺序和方向。

5.本地覆盖异常区域的识别

覆盖分析:
- 将417.9 Gb质量合格的超长ONT读取（长度超过10 kb）与基础Mo17组装进行比对，使用Minimap2。
- 根据特定FLAG值和查询覆盖率选择主要和补充比对。
- 使用BEDTools（v2.29.2）计算每个碱基的读取深度。
- 根据1 kb的区间计算平均深度，定义低覆盖区域（LCR）和高覆盖区域，并根据需要合并相邻区间。

6.基础Mo17组装与PacBio组装的整合

使用Minimap2（参考文献78）将PacBio Hifiasm和Canu组装的contig比对到基础Mo17组装上，参数设置为‘-x asm5 -f 0.02’。根据比对结果，将两类PacBio contig整合到基础Mo17组装中，包括能够填补空缺的contig和用于修正基础Mo17组装错误的contig。

7.TAG重复阵列和45S rDNA阵列相关缺口的闭合

根据超长的ONT读取手动闭合了5个与TAG重复阵列相关的缺口，并基于PacBio HiFi读取闭合了与45S rDNA阵列相关的缺口，具体方法详见补充材料。

8.FISH实验

收集新鲜的Mo17花序，并用Carnoy溶液（乙醇：冰醋酸，3:1，体积比）固定。按照参考文献80中的方法进行FISH实验。FISH探针（TAG）15和（TTTAGGG）7（端粒重复序列）分别标记为digoxigenin和Cy5。通过与Rhodamin（Roche）结合的抗-digoxigenin抗体检测digoxigenin标记的探针。将Mo17的花药压碎，并通过相差显微镜进行分期。选择有早期细胞期的切片进行FISH实验。使用配备有CCD相机的Olympus BX61荧光显微镜（QImaging; RETGA-SRV FAST 1394）捕捉细胞学图像。使用Image-Pro Plus 6.0软件（Media Cybernetics）分析数字图像。

9.rDNA拷贝数估计

通过基于BLAST的方法使用ONT超长数据和PacBio HiFi数据，采用k-mer基于的方法使用Illumina无PCR数据，估计基因组中5S和45S rDNA的拷贝数。此外，使用Naica Crystal Digital PCR System（Stilla Technologies）估计基因组中45S rDNA的拷贝数，具体方法详见补充材料。

10.使用ONT数据验证T2T Mo17组装

将417.9 Gb质量合格的超长ONT读取（长度超过10 kb）使用Minimap2（参考文献78）比对到T2T Mo17组装上，参数设置为‘-x map-ont -r 10000 -N 50’。根据比对结果，识别出合适的比对读取，标准为主比对和辅助比对（SAM格式文件中的FLAG值必须为0、16、2048或2064），且查询序列覆盖率最低为0.85。识别融合读取的标准是该读取比对到多个基因组区域，且所有比对的查询覆盖率合计达到0.85。识别对称读取的标准是该读取可大致分为两部分，序列互为反向互补，并且两部分对齐到同一区域，查询覆盖率合计达到0.85。仅使用合适的ONT读取进行读取覆盖分析。使用BEDTools（v2.29.2）计算T2T组装每个碱基的读取深度，参数为‘genomecov -bga -split’。然后计算基因组中所有1 kb bin的平均深度。将深度低于100或高于250的bin标识为局部覆盖异常。相邻的bin被合并。除了合适的比对读取、融合读取和对称读取外，剩余的ONT读取进一步被判断为源自细胞器基因组或嵌合读取（补充材料）。

11.T2T Mo17组装的验证

使用PacBio和Illumina数据

精度估算：
- 使用Merqury（v1.1）从映射的k-mer估计T2T Mo17组装的准确性。
- 使用251 Gb的Illumina PCR-free数据生成k-mer数据库（K=21）。
- 在Mo17组装中，每个k-mer的存在性与Illumina数据库进行比较，缺失的k-mer被视为基于的“错误”。
- 发现2178604120个k-mer中有2676840个仅在组装中检测到，计算得出质量值(QV)为42.3252。
- 通过以下公式得出Mo17组装的准确性：
  
  100−(10(42.3252/−10))×100=99.9941100−(10(42.3252/−10))×100=99.9941
覆盖异常区域识别：
- 使用151.1 Gb的PacBio HiFi读取对T2T Mo17组装进行比对，采用Minimap2（参数：‘-x map-pb -r 1000 -N 50’）。
- 根据比对结果，识别出映射正确的读取，标准为主对齐和补充对齐，最小查询序列覆盖率为0.85。
- 使用BEDTools计算每个基因组位置的读取深度，并识别深度低于20或高于105的区域作为局部覆盖异常。

12.重复元素注释

转座子元素库：
- 从玉米TE联盟手动整理转座子元素库（maizeTE02052020），用于识别Mo17基因组中新转座子的注释。
- 使用EDTA（v1.7.0）进行新转座子的识别，参数为‘-species maize -curatedlib maizeTE02052020’。
重复元素的发现与标注：
- 使用RepeatMasker（v4.1.1）发现和识别Mo17基因组中的重复元素。
- 使用Tandem Repeats Finder（TRF，v4.09.1）识别微卫星、迷你卫星和卫星，并对结果进行冗余剔除。
图形表示：
- 不同颜色被分配给注释中不同类型的重复元素，以便在基因组区域或ONT读取的图形表示中进行可视化。

13.转座子阵列的识别

将所有基因间区域划分为100kb的窗口，并计算每个窗口的转座子比例，比例超过0.95的被视为转座子阵列，重叠的阵列合并。

14.生成ISO-seq数据

为了帮助基因组注释，生成了来自五种不同组织的混合RNA的ISO-seq数据，包括根、丝、穗和苞片（收集于抽穗阶段），以及14天的幼苗。
提取RNA的纯度、浓度和完整性，以确保使用高质量的RNA进行文库构建。
使用ab initio预测和基于证据的预测方法对Mo17基因组中的蛋白编码基因进行预测，方法包括RNA测序、ISO-seq、蛋白质同源搜索和MAKER预测。

16.同源基因分析

使用OrthoFinder（v2.5.2）进行Mo17基因组和NAM创始系基因组的同源性分析。
通过BLASTP对Mo17基因组中的着丝粒基因进行同源基因识别，E值阈值设定为小于1×10^-3。

17.基因重复的识别

使用 BLASTP 计算每个基因最长转录本编码的蛋白质之间的成对相似性，设定E值阈值为 < 1.0 × 10⁻²⁰。
然后使用 MCscanX 进行分类，采用默认参数。

18.几乎相同的片段重复的识别

将 T2T Mo17基因组 进行自身比对，使用 Mummer 工具。
从Mummer的结果中，识别出相似度超过0.99且长度超过1,000 bp的区域对作为几乎相同的片段重复。
将重叠和相邻的几乎相同片段重复合并。

19.卫星数组的识别

使用 BLAST 识别与五个已报道的卫星重复（knob180、TR-1、CentC、Cent4和tRNASAT_ZM）以及三个新识别的卫星重复（sat268、sat261和sat112）在Mo17基因组中的匹配。
knob180、TR-1和CentC的序列来源于 MTEC，Cent4和tRNASAT_ZM的序列来源于 Repbase。
将基因组分成100 kb的窗口，步长为10 kb，计算每个窗口中八种卫星的比例。
对于每种卫星，窗口中序列超过10%的定义为数组。将相同类型的重叠数组合并，并修剪合并后的数组边界，以去除末端非卫星序列。

20.高阶重复分析

为knob180、TR-1和CentC重复生成位置概率矩阵（PPM），并根据先前报道的方法计算每个卫星重复副本与PPM的变异距离。
所有分析的卫星完整拷贝使用 MAFFT 进行多序列比对，并计算每个比对位置的核苷酸频率以生成PPM。
对于每个单体，通过将每个位置的不同核苷酸概率求和（即1减去相应的核苷酸频率）计算“变异距离”。
每种卫星的单体与其他拷贝进行比较，变异得分为5或更少的单体被聚类为同一组，称为高阶重复组。一个单体可以同时属于不同的高阶重复组。

21.5S和45S rDNA的基因型分析

从最终的T2T Mo17基因组组装中提取每个完整的5S和45S rDNA拷贝序列。
分别使用 MAFFT 进行多序列比对。
对于5S rDNA，选择支持超过10% 5S rDNA拷贝的SNPs和插入缺失（indels）用于基因型分析。
对于45S rDNA，仅选择支持超过10% 45S rDNA拷贝、且大于5 bp的插入缺失用于基因型分析，所有选定的插入缺失位于间隔区。

22.中心粒的识别

使用兔多克隆抗体抗CENH3进行抗体染色质免疫沉淀测序（ChIP-seq）库的构建，抗体针对的肽序列为RPGTVALREIRKYQKSSTSATPERAAGTGGR。
收集约10 g的新鲜Mo17叶片进行细胞核提取，提取的细胞核用微球菌核酸酶消化，然后使用抗CENH3进行ChIP。
每25 µg的染色质使用5 µl的抗体（浓度为0.83 mg/ml）。设置两个生物重复。
ChIP-seq库在Illumina HiSeq平台上进行测序，生成150-nucleotide成对读取。
原始读取进行接头修剪和过滤，使用 fastp 进行处理，随后将结果映射到Mo17基因组，使用 Bowtie2 进行映射。
提取唯一映射的读取，用于进一步分析。使用 bamCompare 计算CENH3的富集水平，保留富集水平大于5的bins，并合并距离小于1 Mb的bins，最终通过可视检查确定Mo17基因组的功能性中心粒区域。