Assembly of diverse Chinese genomes
68个样本:代表了36个少数民族和8个具有语言群体的中国泛基因组联盟核心样本
hifi:28.82X(14.29–60.67×)
9个sample 测了 linked reads;
11个sample测了ONT三代以及hic数据
组装了58个单倍型水平的组装体(图 1a)
同时,作者也整合了5个高深度测序的ONT长读长测序样本,4个中国人群的已经发布的5个组装体,以及2个华中汉族样本(ONT长读长测序)
使用Hifiasm 软件对 68 个 PacBio HiFi 样品组装,初步组装体以及单倍型组装
从初步组装去除了 3 个样品,
从二倍体组装体中去除了 7 个样品
最终保留了 58 个样品或者说是 116 个高质量组件用于进一步分析
Showing here are the geographical locations of all samples included in the CPC
Phase I. See Supplementary Table 1 for details of the CPC sample information.
The China map used in this study was obtained from a standard map service
(http://bzdt.ch.mnr.gov.cn) approved by the Ministry of National Resources of
the People’s Republic of China (GS [2020] 4618).
我们评估了 116 个组装体总基因组长度平均为 3.01Gb的组装,由于性染色体之间的大小差异,范围从 2.88 Gb -3.12 Gb,93.1% 的组装体比无gap的GRCh38 基因组长度2.94 Gb更长
图1b 和补充表 2)
The contigs of T2T-CHM13 and GRCh38 (N-masked) are included for comparison.
116 个组装体平均包含 690.5 个contigs,294 -1568 个(补充表 2)。在 116 个组装体中,contig N50长度中位数从 11.66 Mb 到 87.3 Mb 不等,平均为 35.63 Mb,8.62% 的组装体contig N50 值大于 GRCh38 (57.88 Mb;图1b 和补充表 2)。
我们表明,超过 99.67% 的CCS reads可以mapping到每个组装体中的contigs,116 个组装体的平均质量值为 52.84,在 47.33 和 59.45 之间变化(图 1c)。
The x axis shows the quality value of each assembly, and the y axis shows the rate at which the circular consensus sequencing reads were remapped to the assembly contigs.
所有contig大约24.34 Mb(17.9 -43.01 Mb)未比对到 T2T-CHM13,约 69 Mb(53.32 -89.05 Mb)未比对上GRCh38 ,表明 CPC 样本的基因组未被两个参考基因组中的任何一个完全覆盖(图 1d)。
x轴代表代表比对不上参考基因组(GRCh38 and T2T-CHM13 references)的长度,y轴代表每个组装体比对到参考基于组上的比例!
我们的组装体覆盖了 GRCh38 的 96.54% (92.55-98.03%) ,T2T-CHM13 基因组的 93.59% (89.66-95.77%) (图 1d),表明我们组装体的一些基因组区域是系统地未组装出来的或无法可靠地比对上,这可能是由于高度重复的区域。
此外,我们对contigs的未比对上的区域进行了注释,发现大约 84%(约 39.7-75.1 Mb; 平均58.1 Mb)的未比对上的序列是卫星重复序列(补充表 3)。接下来,我们应用 Inspector 来揭示每个装配中平均有 3,627 个小规模装配错误 (1,889-6,848),包括 327 个扩展 (143-610)、454 个折叠 (154-898) 和 2,846 个替换 (1,287-5,663),这些错误在assembly polishing后减少(图 1e 和补充图3). 每个组装体平均有 38 个大尺度装配误差 (11-78) (补充图 4)。
然而,我们发现了一个有趣的案例,当测序错误和组装错误共存时,评估软件无法确定正确的组装(补充图 5),指出当前评估工具在处理具有挑战性的情况时的局限性。我们最终分别评估了相对于 GRCh38 和 T2T-CHM13 参考的重复率,它们在 58 个样品的 116 个组装体中大致处于同一水平(图1f)
明暗颜色代表每个组装提相对于参考基因组的重复率
Genomic features of the CPC assemblies
为了研究我们 58 个 CPC 样本的基因组组成,我们将 116 个组装体与 T2T-CHM13 对齐,并使用 Phased Assembly Variant caller17 (Methods) 检测每个样品的变体。我们观察到每个 CPC 样本平均有 11.22% (9.35-13.05%) 的小变异和 24.16% (19.72-29.3%) 的结构变异 (SV) 位于高度重复的区域(Extended Data Fig. 2),并且除了 Y 染色体上的 SV 外,超过一半的 SV 是singletons(补充图 6)。
singletons指的是仅在单个CPC样本中观测到的独特结构变异。这意味着这些结构变异在整个CPC样本集中是孤立存在的,即仅出现在一个样本中,其他样本中没有类似的变异。
The highly repetitive sequences in each assembly were identified using dna-brnn. Proportions of (a) small variants and (b) SVs in these regions are shown for each sample. In each plot, the samples are ordered by the total variant proportion within each language family.
使用 dna-brnn 鉴定每个组装中高度重复的序列。展示了每个样本在这些区域中 (a) 小变异和 (b) SV 的比例。在每个图中,样本按每个人群中检测到的变异的数目比例排序。
Different shadows in each bar indicate SVs shared among all CPC samples
that not present in T2T-CHM13 reference. These varians are classified into
major allele variants (AF ≥ 50%, purple), common allele variants (5% ≤ AF <
50%, purple), polymorphisms (≥ 2 haplotypes but AF < 5%), and singletons.
每个柱子中不同阴影代表CPC样本中共享但是在T2T参考基因组中不存在的SV
主要等位基因变异(AF ≥ 50%):群体中常见的主要变异。
常见等位基因变异(5% ≤ AF < 50%)普遍存在但未达到主要等位基因的水平。
多态性变异(AF < 5% 且至少有两个单倍型):罕见但多样化的变异。
单一变异(singletons):只在单个样本中出现的变异。
与 T2T-CHM13 参考相比,我们通过收集大于1kb的insertion插入片段,进一步确定了每个 CPC 样本的平均2802 (2265 –3187)个新序列(Methods and Extended Data Fig. 3a)。这里是2802条序列,不是bp,长度是204Mb。
这些新序列位于基因组的 63,243 个区域,它们在高度重复的序列中特别丰富,例如着丝粒和端粒区域 (Extended Data Fig. 3b) 我们根据插入频率确定了 115 个基因组区域作为插入热点,总长度为 204.8 Mb。
具体是怎么做的,需要看方法部分!!!
a, Number of novel sequences (insertions ≥1 kb) identified by aligning two-phased assemblies of each CPC sample to the T2T-CHM13 reference. b, Chromosome distribution of 115 insertion hotspots of novel sequences.
a,通过将每个CPC样本中的两倍型组装体与T2T参考基因组比对鉴定到的新序列。b,鉴定到的新序列的115个插入位点染色体分布图。
考虑到reference bias通常会阻碍无论长还是短reads数据中的新 SVs 识别,于是使用CPC 高质量和分单倍型序列进一步完成了复杂结构变异 (CSV) 的分析。on average,我们使用新开发的方法 SVision[18]. SVision: a deep learning approach to resolve complex structural variants. Nat. Methods,每个样本鉴定出 70 个 CSV(补充Fig. 7), 最常见的 CSV 类型是 INS+tDUP,平均数量为 44 个(补充Fig. 8)
The type of INS+tDUP is the most abundant CSV across all samples. Others
indicate the CSV types with count smaller than five in each sample.
INS+tDUP是所有样本中最多的CSV,图例中的others代表的是每个样本中CSV类型少于5种!
我们将 CSV 合并为一组 706 个非冗余 CSV(Extended Data Fig. 4a) 将它们分为四类:共享shared(所有样本)、主要major(≥50% 的样本)、多态性polymorphic(<50% 的样本)和singleton(只有一个样本)。我们发现 CSV 长度分布的峰值在 10 kb 处(Extended Data Fig. 4b)。我们从 CPC 汇编集中发现了一个共享的 CSV(chr11:66,245,313–66,246,097 个,INS+INV 类型)、34 个主要 CSV、277 个多态 CSV 和 394 个singleton CSV(Extended Data Fig. 4c)
a, The count of CSVs with different shared sample numbers. b, The length distribution of the final 706 CSVs. We merged CSVs of all 58 samples by considering 80% reciprocal overlap and the same type, and obtained 706 CSVs in total. c, The number of CSVs in each CSV class. We classified the CSVs by the shared sample numbers into Shared (present in all samples), Major (present in ≥50% samples but not all), Polymorphic (present in more than 1 but <50% of all samples), and Singleton (present in only one sample).
a,不同共享样本数时的SCV数;b,706个CSVs的长度分布;合并了58个样本检测到的CSVs,考虑互相重叠的以及相同类型的CSVs,共获得706个CSVs;
c,不同CSVs分类下的CSVs数,我们把所有样本直接中都出现的CSVs定义为Shared,≥50% 样本中都出现的CSVs定义为Major,后续还有多态以及Singleton(只在一个样本中出现的CSVs)
接下来,我们注释了每个组装中相对于 GRCh38 基因组的拷贝数变异 (CNV)。
具体是怎么做的?method
在完整的组装集中有 1,367 个蛋白质编码基因,这些基因在至少一个基因组中的拷贝数有所增加。在每个组装中平均观察到 53 (27-100) 个拷贝数增加的基因(图 2a),13.39% 的 CNV 基因在整个 CPC 组装集中的频率超过 5%,57.86% 的 CNV 基因仅存在于单个单倍型中。我们发现 CPC 组装集中有 1,079 个重复基因,这些基因在 HPRC 组装中没有观察到(图 2b和补充表 4)。这些基因在嗅觉相关功能中显著富集 (例如,GO:0007608,嗅觉的感觉感知,比值比 (OR) = 6.46,Benjamini-Hochberg (BH) 调整的 P = 2.22 × 10-36;和 hsa04740,嗅觉转导,OR = 7.14,BH 调整后的 P = 8.35 × 10-38),并且在皮肤相关疾病中处于边缘显著性水平(例如,DOID:37,皮肤病,OR = 2.33, BH 调整后 P = 0.086;和 DOID:16,外皮系统疾病,OR = 2.23,BH 校正 P = 0.086;补充表 5)。
a, Number of duplicated protein-coding genes per CPC genome assembly compared with the GRCh38 reference.
和hg38参考基因组相比每个CPC组装体和复制相关的基因数
b, Venn diagram showing the number of duplicated genes in CPC, HPRC.EAS
and HPRC.nEAS assemblies.
韦恩图展示三个不同群体的duolicated 基因数
c, The top 20 most common CPC-specific CNV-related genes compared with the HPRC assemblies.
和人类泛基因组联盟相比,中国泛基因组联盟发现的top20最常见CNV相关的基因
d, The five overlapped CNV genes showing a higher frequency (≥5%) in CPC assemblies (blue) than in HPRC assemblies (orange). HPRC.EAS, East Asian in HPRC; HPRC.nEAS, non-East Asian in HPRC.
五个有重叠的疾病相关CNV相关基因在CPC中的频率更高(≥5%)
我们发现,根据全基因组关联研究 (GWAS) 目录,CPC 组装集中的新重复基因中有 562 个 (52.1%) 是性状相关的,其中 207 个 (19.2%) 基因与疾病本体注释的至少一种人类疾病有关,这表明 CPC 组装集在促进疾病和表型关联研究方面具有显着的潜力。 特别是,KCNJ18 和 MTRNR2L6 在 CPC 组件中显示出相当高的频率(116 个中的 115 个),但在 HPRC 组件中不存在(图 D)。KCNJ18 编码内向整流钾通道家族 (Kir2.6) 的一个成员,该家族在骨骼肌中表达,并受甲状腺激素的转录调控。据报道,Kir2.6 的变化与甲状腺毒性周期性麻痹有关,这是亚洲人群中众所周知的甲状腺毒症并发症,但在欧洲人群中很少见。假基因 MTRNR2L6 被预测参与细胞凋亡的调节 (GO:1900118)。值得注意的是,据报道,该基因在地理汉族人群中显示出正选择的特征。许多 CPC 特异性 CNV 基因显示出自然选择的特征,由Tajima’s D 估计值表明(Supplementary Fig. 9 and Supplementary Table 4))。
Tajima’s D was calculated for the 58 CPC samples (116 haploid assemblies)
with a sliding window spanning 20 kb in size, stepping 10 kb. The estimate for
each gene was represented by the top window overlapping that gene. The
horizontal dashed line shows the threshold of P = 0.05, which is obtained by
the two-tailed significance test based on the beta distribution, and is further
adjusted for multiple testing using the FDR method; the vertical dashed line
indicates the CNV counts ≥ 5.
在20kb大小的窗口中以滑动步长10kb计算CPC的58个样本(116个单倍型)Tajima’s D,如果一个基因存在多个窗口的重叠,以最高Tajima’s D得分的窗口为准,基于beta分布的双尾显著性检测获得水平虚线标注的显著线阈值P = 0.05,并且进一步经过FDR校准用于多重测试;垂直虚线代表CPC样本中大于5的CNV
数。
典型例子包括 SPATA31D4 (n单倍型 = 6,Tajima’s D = −3.01, false-discovery rate (FDR)-adjusted P = 2.38 × 10−5)和 PSAPL1 (n单倍型 = 7,Tajima’s D = 4.32, FDR-adjusted P = 0.006)。特别是,SPATA31D4 属于核心duplication家族,该家族被认为对类人猿进化做出了重大贡献,在类人猿中具有强烈的正选择信号。 先前的研究表明,该基因家族具有响应紫外线辐射和 DNA 修复的功能,也可能影响人类寿命。
在 HPRC 组装中也观察到来自 CPC 组装的288 duplicated genes,其中 123 个基因与 4 个 HPRC 东亚样本 (HPRC.EAS)共享,并且 278 个基因与其余 40 个非东亚 HPRC (HPRC.nEAS) 样本共享(Fig. 2b)。在这些重叠的基因中,我们发现 CPC 组装集中的几个基因的频率 (≥5%) 高于 HPRC 组装(Fig. 2d and Supplementary Table 6)。例如,CPC 程序集中 CYP2D6 处的 CNV 频率高于集体 HPRC 程序集和单独的 HPRC 集中的 CNV 频率。EAS 或 HPRC.nEAS (Fig. 2d and Extended Data Fig. 5)。CYP2D6 以其显着的多态性而闻名,已在各种人类群体中进行了系统调查。我们还发现 CFC1 在 CPC 和 HPRC 组装体之间具有高度差异化的 CNV (频率差 = 0.124)。几项调查报告了 CFC1 与心血管疾病之间的关联,尤其是在中国人群中。此外,CTAG1A 和 ZNF658 在 HPRC 组装中鉴定出罕见的 CNV (约 1%),但在 CPC 组装中以常见的 CNV (≥5% 和 <95%) 呈现。CTAG1A 在正常卵巢和睾丸中高水平表达。ZNF658 调节锌匀态相关基因的转录并影响核糖体生物发生,并参与感染过程 (hsa05168,单纯疱疹病毒 1 感染)。
The variation graph of the CPC pangenome
在典型的泛基因组reference中,来自群体的基因组数据可以组织成基于边缘的序列变异图。我们应用 Minigraph-Cactus 管道构建了 CPC 泛基因组的变异图,其中单倍型组装体可以表示为由序列节点组成的不同路径(Fig. 3a)。使用 Minigraph,从参考基因组T2T-CHM13 和 GRCh38 开始创建graph,把 61 个样品的 122 个单倍型组装体添加到graph中,为了将泛基因组扩展到单核苷酸多态性水平,使用 Cactus (一种基因组aligner)将组装体重新比对到graph上。我们使用 dna-brnn 来识别和屏蔽输入组件的着丝粒卫星,然后clip所有长度超过 100 kb 且未与underlying graph比对上的路径,以避免复杂区域对比对alignment的影响。最终我们构建了 Minigraph-Cactus 图,其长度为 3284609818 bp(以所有节点的总和衡量)。此外,我们还测量了从 61 个二倍体基因组中的每一个添加到泛基因组图中的序列的长度。图中总共添加了 194.67 Mb 的非参考序列,其中 62.90 Mb 是singletons。此外,我们发现 4.96 Mb 的非参考序列存在于 ≥95% 的所有单倍型中,代表了我们样本群体中的核心基因组,72.24 Mb 的非参考序列存在于 ≥5% 和 <95% 的单倍型中,代表了common基因组(图 3b)。
a, A variation graph representing CHM13.chr11:34,823,712–34,823,777 of the CPC pangenome. Coloured lines represent five haplotype assemblies. b, Pangenome cumulative growth curves for the CPC pangenome graph. Depth in bars measures how often a non-reference segment occurred in the assembled haplotypes. Non-reference sequences were classified into four categories according to the frequency of presence in the haplotypes: core (≥95%), common (≥5% and <95%), polymorphism (≥2 haplotypes but <5%) and singleton (only one haplotype). Different colours indicate a different ethnic group. c, The number of CPC-specific and common variants between CPC and HPRC assemblies from the joint pangenome graph. d, Number of identified CPC-specific small variants and SVs in different populations. e, The distribution of CPC-specific SVs from the joint pangenome graph on autosomes. The colour scale from white to red represents the density of SVs in 1-Mb regions. The orange triangle marks CPC-specific SVs that were significantly enriched in this 1-Mb window compared with HPRC-specific variants and common variants, on the basis of a one-tailed Fisher’s exact test (FDR-correlated P value < 0.05).
a,代表 CPC 泛基因组的 CHM13.chr11:34,823,712–34,823,777 的variation graph。彩色线条代表五个单倍型组装体。b,CPC 泛基因组图的泛基因组累积生长曲线。Depth in bars 代表非参考片段在组装的单倍型中出现的频率。根据单倍型中存在的频率将非参考序列分为四类:核心 (≥95%)、常见 (≥5% 和 <95%)、多态性 (≥2 个单倍型但 <5%) 和singleton (只有一个单倍型)。不同的颜色表示不同的种族群体。c, 来自联合泛基因组图的 CPC 和 HPRC 组装之间的常见变异以及CPC 特异性的数量。d,不同人群中已鉴定的 CPC 特异性小变异和 SV 的数量。e,来自常染色体上联合泛基因组图的 CPC 特异性 SV 的分布。从白色到红色的色标表示 1 Mb 区域中 SV 的密度。橙色三角形标记了根据单尾 Fisher 精确检验 (FDR 相关 P 值 < 0.05) 与 HPRC 特异性variants和常见variants相比,在该 1 Mb 窗口中显著富集的 CPC 特异性 SV。
由于underlying的122个单倍型组装被编码为图中的路径,因此我们通过图分解(method)表征了组装的变体variants和单倍型。我们进一步检查了最长的等位基因是否出现在每个位点,并从图中确定了 15,916,404 个小变异位点 (<50 bp) 和 78,072 个 SV 位点 (≥50 bp)。每个样本包含 4397706 (s.d. = 35,776) 小变异,每个单倍型包含14557 (s.d. = 224) SV(Supplementary fig.10)
什么叫最长的等位基因? longest allele是否出现在每个位点?什么位点?
这里应该是变异位点,原文用的是variants site, each site, small-variant sites, SV sites等描述
基于 HiFi 数据从每个样本中获得的 SV 数量远远多于基于下一代测序 (NGS) 数据(Supplementary Table 7)。CPC 图基因组中 SVs 的平均长度为 3.2 kb,中位长度为 178 bp,比 NGS 数据中的短读长长。此外,我们在 Alu 的 300 bp 左右和 LINE-1 的 6 kb 处观察到一个峰值(扩展数据图 .6),表明我们的 CPC 泛基因组的可靠性。
The average and median length of SVs are indicated by the two red vertical lines, and both are larger than the read length (usually 150 bp) generated by short-read sequencing. The peaks at 300 bp for Alu insertions and 6 kb for LINE-1 are highlighted.
SVs的长度和均值以及中位数用两条红色垂直线表示,都比短reads的150bp长。
Alu的300bp插入和LINE-1的6kb的变异在图中已经标出!
Alu和LINE 1是研究最多的两类non-LTR,非长末端重复逆转录转座子,后者一般是6kb,前者一般为几百bp
CPC 核心样本在很大程度上覆盖了中国多民族人口的遗传多样性,并显示出比汉族样本更多样化的遗传血统(Extended Data Fig. 7).
We randomly selected 10 unrelated high-quality samples from each CPC population, and inferred the genetic ancestry for each sample using ADMIXTURE assuming 2–12 ancestry components (K). Samples used in the pan-genome construction are labeled with short vertical lines.
我们从每个 CPC 人群中随机选择了 10 个不相关的高质量样本,并使用假设 2-12 个祖先成分 (K) 的 ADMIXTURE 推断每个样本的遗传祖先。泛基因组构建中使用的样本用短垂直线标记。
这个结果的解读和意义
为了评估多民族人群在 CPC 泛基因组中的贡献,我们通过单独使用汉族样本和多民族人群来比较非参考序列的累积长度。当添加第一个单倍型时,汉族样本产生的非参考序列的长度与其他人群产生的非参考序列的长度大致相同(Extended Data Fig. 8)汉族人的平均非参考序列长度为 49.71 Mb,其他人群为 49.16 Mb。然而,通过添加单倍型,我们获得了截然不同的结果。特别是,在多民族人群中检测到的非参考序列的累积长度比在单个人群中检测到的序列(例如汉族人)的累积长度增长得更快(Extended Data Fig. 8)
The cumulative length of non-reference (unaligned to the human reference genome GRCh38) sequences in Han Chinese haplotypes (the red boxes) and that in multi-ethnic haplotypes (the green boxes) are calculated based on the CPC graph genome by pangenome-growth software. We randomly ordered the three Han Chinese samples, and conducted 40 replication analyses with randomly selected samples (three samples for each analysis) from the CPC dataset. The multi-ethnic haplotypes show a relatively larger growth rate of the non-reference sequences than the Han Chinese samples. The lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles). The upper whisker extends from the hinge to the largest value no further than 1.5 * IQR from the hinge (where IQR is the inter-quartile range, or distance between the first and third quartiles). The lower whisker extends from the hinge to the smallest value at most 1.5 * IQR of the hinge. Data beyond the end of the whiskers are called “outlying” points and are plotted individually.
汉族单倍型(红框)和多种族单倍型(绿框)中非参考(与人类参考基因组 GRCh38 不对齐)序列的累积长度是通过泛基因组生长软件(pangenome-growth software)基于 CPC 图基因组计算的。我们随机排序了 3 个汉族样本,并使用从 CPC 数据集中随机选择的样本(每次分析 3 个样本)进行了 40 次重复分析。多种族单倍型显示非参考序列的增长率相对大于汉族样本。下铰链和上铰链对应于第一和第三个四分位数(第 25 个和第 75 个百分位数)。上须从铰链延伸到距铰链不超过 1.5 * IQR 的最大值(其中 IQR 是四分位数间距,或第一和第三四分位数之间的距离)。下须从铰链延伸到铰链的最小值,最大为 1.5 * IQR。超出须线末端的数据称为“外围”点,并单独绘制。
Short-read mapping with the CPC reference
为了评估 Giraffe 比对软件处理不同复杂程度graph genomes时的性能,我们根据节点的路径深度path deoth对 CPC reference进行梯度过滤(Methods)。当过滤更严格时,我们观察到图形参考中的节点和边的数量有所减少(Supplementary Fig. 11a),图复杂度(Supplementary Fig. 11b) 和多样性(Supplementary Fig. 11b)有所减少
(a) The number of nodes and edges in graph. The blue numbers show the filter
criteria (vg clip -d). (b) The complexity of graph genome. It is obtained by
dividing the number of edges by the number of nodes. (c) The length of
sequences in graph, i.e., the sum of the lengths of all nodes.
蓝色数字代表过滤标准,越往左标准越严格!
图的复杂性用边/节点数表示!
图的序列长度代表着图的多样性!
接下来,来自 1000 个基因组计划的东亚人群的 10 个样本通过 vg Giraffe 与这些具有不同复杂性的 CPC 图参考基因组进行比对。结果表明,使用简化版的图表(Supplementary Fig. 12a)时比对率增加并达到峰值,可能是由于当前版本的 Giraffe 比对软件在处理局部复杂区域方面的局限。
(a) Mapping rate. One sample in each figure. The abscissa shows different
filtering criteria. (b) Proportion of perfect alignment is evaluated through vg stats.
横坐标代表不同的过滤标准;完美比对通过vg stats计算得到!
我们还观察到,随着图形的简化,具有完美匹配的reads(perfect matching)比例不断下降(Supplementary Fig. 12b),反映了图形多样性的下降(Supplementary Fig. 13)。因此,比对率和比对质量之间存在权衡,需要折衷来确定图形引用的大小。
The abscissa is the sequence length of graph genome. Larger represents more
alternative sequences.
横坐标代表图基因组的序列长度,越大代表越多的可变序列!
与 HPRC 图基因组相比,CPC 图基因组具有更少的节点、边和多样性,这可能是由于仅包括中国样本,与同时涵盖非洲和欧洲样本的 HPRC 相比(Supplementary Table 8)。然而,在比对东亚基因组时,使用 CPC 图基因组比使用 HPRC 图形参考实现了更好的比对(Supplementary Fig. 14)。相比之下HPRC 图基因组在处理非洲样本方面表现更好(Supplementary Fig. 14)。结果表明,使用人群特异的图基因组提高了短读长的比对质量。
Each figure shows a sample. The red line represents the East Asian sample
and the black line represents the African sample. (a) Use raw CPC graph
reference and raw HPRC graph reference. (b) Use filtered CPC and HPRC
reference. The filtering standard is the minimum frequency of 0.1.
红色代表东亚人口样本,黑线代表非洲样本;
a用的是原始未过滤的图基因组;b用的是最小等位基因为0.1过滤的
为了进行变异检测,我们将图形参考坐标中的 GAM 文件映射到线性参考坐标中的 BAM 文件。结果显示,所有样本的比对率平均下降了0.58%(0.54-0.61%)。我们推测,当使用传统的线性参考来进行call或检测变异时,图参考的优势将会丧失,因为图参考中的novel序列在线性坐标中丢失。
这段到底是怎么做的!!!!映射比对!!!
Comparison with the HPRC pangenome graph
为了调查 CPC 全基因组图中东亚人群贡献的先前未识别的组件,我们构建了一个合并的 Minigraph-Cactus 图,包括 CPC 中的所有 116 个组件和 HPRC 中的 94 个组件(Methods)。我们识别出仅在 CPC 组件中存在的 5850863 (18.4%) 个小变异和 34223 (17.1%) 个的 SV(Fig. 3c),其中每个样本包含170307(sd = 10,904)个小变异,每个单倍型携带543(sd = 39)个SV,并且超过一半的CPC特异性变异是单联体或双联体(singletons or doubletons )(Fig. 3d )
a,代表 CPC 泛基因组的 CHM13.chr11:34,823,712–34,823,777 的variation graph。彩色线条代表五个单倍型组装体。b,CPC 泛基因组图的泛基因组累积生长曲线。Depth in bars 代表非参考片段在组装的单倍型中出现的频率。根据单倍型中存在的频率将非参考序列分为四类:核心 (≥95%)、常见 (≥5% 和 <95%)、多态性 (≥2 个单倍型但 <5%) 和singleton (只有一个单倍型)。不同的颜色表示不同的种族群体。c, 来自联合泛基因组图的 CPC 和 HPRC 组装之间的常见变异以及CPC 特异性的数量。d,不同人群中已鉴定的 CPC 特异性小变异和 SV 的数量。e,来自常染色体上联合泛基因组图的 CPC 特异性 SV 的分布。从白色到红色的色标表示 1 Mb 区域中 SV 的密度。橙色三角形标记了根据单尾 Fisher 精确检验 (FDR 相关 P 值 < 0.05) 与 HPRC 特异性variants和常见variants相比,在该 1 Mb 窗口中显著富集的 CPC 特异性 SV。
在 GIAB 3.0 中定义的 GRCh38 参考基因组的“简单”和“困难”区域 ,大约 39% 的 CPC 特异性小变异无法在 gnomAD v0.1.8 中进行注释(Supplementary Table 9 ),表明用基于长读长的方法鉴定的东亚特有的小变异仍然是对当前基于短读长的遗传资源的有效补充。
Note: To make it comparable to GRCh38-based gnomAD and GIAB, the number of small variants was calculated from the joint MC pangenome vcf file which used CHM13 as the reference (the first assembly) and the GRCh38 (the second assembly) as the coordinate with parameters “--reference CHM13v2 --vcfReference GRCh38” in the “cactus-graphmap-join” step.
为了使得基于hg38注释的gnomAD和GIAB,从joint MC 泛基因组的vcf文件中检测得到小变异的数据,构建时在cactus-graphmap-join这一步中按照参数:--reference CHM13v2 --vcfReference GRCh38,设置CHM13作为参考基因组(第一个组装体),hg38作为坐标(第二个组装体)
coordinate 这个词作为名词时意思就是坐标!!!
我们发现 16898 个(49.4%)CPC 特异性 SV 与 6426 个蛋白质编码基因的附近区域(基因编码区上游和下游各 100 kb)重叠,其中 4344 个基因被SV 跨度超过 1 kb 导致破坏或者打断disrupt,这些基因具有最常见的功能富集,即与免疫功能相关,例如体液免疫反应。 (GO:0002455, OR = 5.11, BH-adjusted P = 8.50 × 10−14; and GO:0006959, OR = 2.91, BH-adjusted P = 1.64 × 10−11; Supplementary Table 10)。 根据疾病GO注释,这些 CPC 特异性 SV 还表现出喉炎相关基因的过表达(DOID:3437 and DOID:786, OR = 16.66, BH-adjusted P = 0.007)。
此外,我们使用基于滑动窗口的常染色体分析估计了 CPC 特异性 SV 的位置分布(Methods)。与HPRC特异性SV和常见的SV类似,大多数CPC特异性SV位于染色体的着丝粒和端粒区域(图3e和补充图15 )。
(a) The quantity difference of SV between CPC and HPRC genomes per 100kb
bin. (b) The number of CPC-specific SV per 100kb bin.
学习如何画100bin范围内的SV数目!!
接下来,我们在不同区域的 HPRC 组件中发现的SV和CPC 特异性 SV 的数量之间应用单尾 Fisher 精确检验,发现了 223 个热点,与其他 SV 相比,CPC 特异性 SV 显着富集(FDR-adjusted P < 0.05),涉及807个蛋白质编码基因(图3e ),过表达的生物学功能如氧运输( GO:0015671, OR = 22.66, BH-adjusted P = 0.008; and GO:0005344, OR = 24.91, BH-adjusted P = 0.001)和血红蛋白结构 (GO:0031838, OR = 28.58, BH-adjusted P = 0.003; GO:0005833, OR = 24.21, BH-adjusted P = 0.003; and GO:0031720, OR = 33.15, BH-adjusted P = 0.002; Supplementary Table 11)。
Note: GeneRatio denotes the ratio of input genes that are annotated in a term, and BgRatio denotes the ratio of all genes that are annotated in a term. The P-values are obtained by one-sided Fisher’s exact test.They are further adjusted for multiple comparisons using the Benjamini-Hochberg procedure (denoted as
the BH-adjusted P-value), and are also adjusted for the false discovery rate (denoted as the q-value).
GeneRatio和BgRatio其实这个主要就是CPC特有的SV和CPC检测到的所有的SV富集到的通路
odds ratio(OR,优势比) 是用来衡量输入基因集中某一功能项富集程度的指标。
若 odds ratio > 1,表示该功能项在输入基因集中富集的程度高于背景基因集,说明该功能在输入基因集中可能具有更重要的作用。
若 odds ratio = 1,表示该功能项在输入基因集和背景基因集中富集程度相同。
若 odds ratio < 1,表示该功能项在输入基因集中富集的程度低于背景基因集。
在表格中,例如 "oxygen transport" 的 odds ratio 为 22.66,这意味着与背景基因相比,该功能项在输入基因集中富集的概率高了 22.66 倍。
长读长测序技术和基于泛基因组图的分析方法使我们能够探索以前难以在NGS数据中定位的大而复杂的SV,从而为这些复杂位点与生理功能或疾病的关联研究提供遗传基础。我们发现上述一些CPC特异性富集的SV与东亚流行的疾病密切相关。一个显着的例子是位于 16 号染色体短臂端粒附近的 α-珠蛋白基因簇,包括 5 个功能基因和 2 个假基因, 5′-zeta–pseudozeta–mu–pseudoalpha-1–alpha-2–alpha-1–theta-3′ (Fig. 4a)。我们根据泛基因组图中的α-珠蛋白基因( HBA1或HBA2 )和ζ-珠蛋白( HBZ或假基因HBZP1 ;Fig.4b)基因的拷贝数变异鉴定了六种主要单倍型 (Supplementary Table 12),主要结构如Fig.4b所示
zeta:ζ-珠蛋白基因:HBZ or HBZP1 (伪基因);
pseudozeta:ζ-珠蛋白基因的伪基因:HBZP1
mu:HBM
pseudoalpha-1:α-珠蛋白基因的第一个伪基因:HBAP1,位于mu基因之后,肺功能性拷贝
alpha-2:HBA2,编码α-珠蛋白
alpha-1:HBA1,编码α-珠蛋白
theta:θ-珠蛋白 ,
血红蛋白亚基 theta-1 (HBQ1) 是一种属于珠蛋白家族的蛋白质。
除了在 CPC 和 HPRC 中发现的涉及 α-珠蛋白拷贝数变化的deletion (Z2A1) 和duplication (Z2A3) 之外,我们还鉴定了两个 CPC 特异性的大 SV:涉及 5 个珠蛋白(Fig.4b中所示的5个珠蛋白,包括伪蛋白)的 20 kb 删除 (Z2A0)基因和涉及 ζ-珠蛋白基因的10-kb 重复(Z3A2 和 Z3A3)(图4c)。两个α-珠蛋白都丢失的长缺失已被广泛报道为东南亚缺失,在我们的报道中是单倍型A0,表示为SEA,主要分布在中国南部和东南亚。正如之前报道的杂合子 SEA 缺失 (A2/A0) 以及 α-珠蛋白基因 (A2/A1) 一个拷贝的丢失在表型上是沉默的。 一个 α-珠蛋白基因 (A1/A1) 的纯合性缺失会导致轻度贫血;丢失三个拷贝 (A1/A0) 会导致血红蛋白 H 病,纯合 SEA 缺失会导致严重的胎儿水肿。 CPC泛基因组图中α-珠蛋白基因簇上复杂SV的精确定位可以为未来贫血相关研究提供潜在参考。
另一个例子是位于7号染色体上的RASA4基因(图4d),与参考基因组的两个拷贝相比,东亚人群中发现了高度多样性的拷贝数(Supplementary Table 13),包括HPRC样本中未发现的六拷贝变异(Fig. 4e)。该基因的 CNV 尚未被描述。由RASA4编码的RAS p21蛋白激活剂4的异常表达已被广泛报道与多种人类癌症的发展密切相关,并且我们观察到人群中剂量频率分布的差异(Supplementary Tables 14 and 15),这可能导致疾病发病率的变化。
a ,α-珠蛋白基因在 CPC 泛基因组子图上的位置。 b ,来自 Minigraph-Cactus 图的 116 个 CPC 单倍体组件和 94 个 HPRC 单倍体组件中所有结构单倍型的等位基因计数和线性结构可视化。图中基因的大小和间距并不代表染色体的实际大小。 c ,不同α-珠蛋白基因单倍型通过联合子图的路径。箭头指示路径的方向。 d ,CPC 子图上RASA4区域中基因的位置。 e ,具有不同RASA4B拷贝数的不同结构单倍型的路径。 “部分”代表RASA4B的 14.9 kb 片段。
接下来我们研究了 CPC 组装体中发现的新型 SV 在多大程度上可以增加我们对疾病遗传学的认识。根据从最新发布的 GWAS 目录中收集的 243465 个表型相关变异(其中在东亚人群中报告或重现了 62393 个变异),我们发现 75.95% 的新型 SV 长度>1 kb(跨越新序列总长度的 83.17%)位于距 GWAS 位点 <50 kb 处,特别是,55.49%(新序列总长度的 72.95%)位于与东亚表型相关的变异周围。我们观察到,在比较已经报道的性状变异时,身高相关变异比其他性状更有可能与大部分(50.7%全球和47.7%东亚人群报告的身高相关变异 )基本上独立的新基因座相关,可能是由于身高显着的多基因性(Extended Data Fig. 9 and Supplementary Table 16)。
Note: Proportion of the novel loci located < 50kb around the GWAS variants was compared with 100 sets of randomly sampled common loci with matched size. Here shows the raw P-values obtained by the onesided Wilcoxon rank-sum test. None of them remains significant after adjusting for multiple comparisons (BH-adjusted P > 0.05).
onesided Wilcoxon rank-sum test. 怎么做?!!!
此外,这些新的 SV 显着富集在与东亚人疾病相关的变异上,包括尿石症、肾结石和甲状腺肿(BH-adjusted P value = 0.043),这些变异在一些亚洲地区非常流行(Extended Data Fig. 9)
将位于 GWAS 变体周围 <50 kb 的近似独立的 CPC 特异性基因座(距离 <50 kb 的相邻新 SV 被合并)的比例与一组随机采样的共同基因座(与 T2T-CHM13 共享)的估计比例进行比较具有相匹配的大小分布。 P值是基于 100 个排列获得的,并且围绕 ( a ) 全球人群中报告的 GWAS 位点和 ( b ) 东亚人报告的 GWAS 位点, CPC 特异性 SV 的富集至少达到边际显着性(BH 调整的P值 < 0.1)。显着的富集在性状标识上用星号 (*) 表示。每个箱线图代表排列检验统计量(灰点)的中位数(粗黑线)、上下四分位数(方框)、1.5×四分位数范围(胡须)。观察到的统计数据用红点表示。
尽管在 CPC 组装体中检测到的所有新 SV 总体上显示出的核苷酸多样性水平(由 Tajima’s D 计算得到)与基因组其他序列(非SV相关序列)相似,但 HPRC 组装中不存在的 CPC 特异性新 SV 表现出明显高于后者的 Tajima’s D (P = 4.65 × 10−9, one-tailed Wilcoxon rank-sum test; Supplementary Fig. 16a)突出的信号包括 2 号染色体中的两个蛋白质编码基因 STARD7 和 ITPRIPL1 (Tajima’s D = −3.06, FDR-adjusted P = 7.54 × 10−12; Supplementary Fig. 16b) 。强调 SPATA31 基因,因为它们不仅可以通过拷贝数变异在类人猿中赋予进化意义(Supplementary Fig. 9) 也与代表性不足的序列变体变异有关 (Tajima’s D = −3.01, FDR-adjusted P = 2.38 × 10−5)。
结合显著的负 Tajima’s D 和校正后的显著性水平,强烈暗示该基因受到了自然选择的作用。所有这些结果都意味着 CPC 组装体中的新型 SV 具有巨大潜力,可以为东亚人类的适应性进化提供新的见解。
Archaic introgression and annotation
古老的渗入和注释
我们应用了 ArchaicSeeker 2.0(并在所有 61个 CPC 样本中鉴定了5338个古老的渗入片段 archaic introgression segments (AISs),总跨度为 703.87 Mb,平均每个样本 84.67 Mb。在这些 AIS 中,2450 个位于 5531 个基因的编码序列中(Supplementary Table 17)。我们发现在至少两个样本中检测到 4126 个 AISs 基因。至少 5 个样品中用 AISs 检测到 2617 个基因,并且显著富集在角化 (GO:0031424, OR = 4.19, BH-adjusted P = 1.23 × 10−5)、I 型干扰素受体结合 (GO:0005132, OR = 130.05, BH-adjusted P = 7.62 × 10−12)、STAT 蛋白肽基丝氨酸磷酸化的正调节 (GO:0033141, OR = 48.22, BH-adjusted P = 1.23 × 10−10),RIG-I 样受体信号通路(hsa04622, OR = 4.04, BH-adjusted P = 2.44 × 10−4)和神经元细胞体(GO:0043025, OR = 1.75, BH-adjusted P = 5.96 × 10−3,Supplementary Fig. 17 and Supplementary Table 18)。在分析所有受 AIS 影响的基因和用高频 AIS 检测到的基因时,我们获得了类似的结果(至少 10 个样本携带的 1510 个带有 AIS 的基因;Supplementary Tables 19 and 20 and Supplementary Figs. 18 and 19)。 极高频 AISs(>40 个样本)影响以下基因(使用 GeneCards Suite在线注释):KRT6C、KRT6A、KRT6B 和 KRT75,它们都是角蛋白基因家族成员; CACNA2D2、CYB561D2、EEF1A2 和 KCNQ2,与发育性和癫痫性脑病相关; GNAT1,在视网膜正常杆状光感受器 (RHO) 介导的光感知中起信号转导作用,与常染色体隐性遗传先天性静止性夜盲症有关; 和 USH2A,参与内耳和视网膜的细胞发育和维持,与中国人群的 Usher 综合征和视网膜色素变性有关。
与 HPRC 数据集中的非洲样本 African samples相比,CPC 组装体富含古老的古人类序列(补充图 20)。
a图是在个体水平上;b图是在群体水平上(对CPC样本中抽取了16100个重复)
我们进一步比较了在 CPC 样本中检测到的 AIS 和在美国样本中检测到的 AIS,它们是 HPRC 数据集中仅次于非洲群体的最大大陆群。我们发现,在个体和种群水平上,阿尔泰尼安德特人样 AISs (Altai Neanderthal-like AISs) 的比例在美国人中(76.44 ± 15.61 Mb)高于 CPC 东亚人(74.39 ± 4.48 Mb),给定可比样本量(美国人为 525.81 Mb,东亚为 434.54 ± 8.42 Mb)。丹尼索瓦人样 AIS (Denisovan-like AIS)比例在东亚基因组中较高(16.92 Mb,覆盖了美国基因组的 0.59%),每个样本为 2.10 Mb ± 0.71 Mb (0.07%);26.04 ± 1.34 Mb,平均覆盖东亚基因组的 0.90%,每个样本为 2.77 Mb ± 0.70 Mb (0.10%)),CPC 东亚群体平均为 26.04 Mb ± 1.34 Mb,覆盖基因组的 0.90%,远高于美洲人群的 16.92 Mb(0.59%),表明东亚基因组中丹尼索瓦人遗传的 AIS 多样性高于美国基因组(补充图 21)。
此外,HPRC 的 CHS 样本很大程度上未能充分体现东亚人中的古人类基因渗入。CPC 组装体中的每个人群平均向现代东亚人的古老序列池中添加了 15.45 Mb 的 AIS(14.16 Mb 的阿尔泰尼安德特人样(Altai Neanderthal-like)序列和 1.29 Mb 的丹尼索瓦样(Denisovan-like)序列)(补充图 21),每个 CPC 基因组贡献了 9.56 Mb 的古样序列。特别是,讲突厥语的人群(例如,维吾尔语、哈萨克斯坦和吉尔吉斯语Uyghur, Kazakh and Kyrgyz)与其他东亚人群的阿尔泰尼安德特人(Altai Neanderthal-like)相似的 AIS 最少,这可能是由于这些人群的欧洲遗传祖先(Extended Data Fig. 10a);一些中国南部语言群体(例如,泰-侗语和南亚语系Tai–Kadai and Austro-Asiatic)最大程度增加了类似丹尼索瓦人的 AIS 多样性(扩展数据图 10b)。
a, Sharing ratio of the Neanderthal-like introgressed sequences. b, Sharing ratio of the Denisovan-like introgressed sequences. The heatmaps were generated using the R package pheatmap 1.0.12. Warm colors in the heatmap indicate higher levels of ancestry sharing, while cold colors indicate lower levels of ancestry sharing. Populations are clustered using the complete-linkage method, and the branches are colored according to language families.
a ,类尼安德特人Neanderthal-like渐渗序列的共享比。 b ,类丹尼索瓦人Denisovan-like渐渗序列的共享比。热图是使用 R 包 pheatmap 1.0.12 生成的。热图中的暖色表示祖先共享水平较高,而冷色表示祖先共享水平较低。使用完全链接方法对群体进行聚类,并根据语系对分支进行着色。
我们进一步研究了受 CPC 特异性 AIS 影响的基因及其潜在功能。我们在 CPC组装体中发现了 1575 个 AIS,大小为 72.41 Mb,而 HPRC 程序集中不存在。
上面应用 ArchaicSeeker 2.0对古基因进行注释
所有 61个 CPC 样本中鉴定了5338个古老的渗入片段 archaic introgression segments (AISs),总跨度为 703.87 Mb,平均每个样本 84.67 Mb。在这些 AIS 中,2450 个位于 5531 个基因的编码序列中。
这些 CPC 特异性 AIS 总共包含 3629 个基因。我们强调了位于编码序列区域的受潜在功能性AIS影响的1211个基因(Supplementary Table 21),这些基因在外源性葡萄糖醛酸化( GO:0052697 ,OR = 77.51,BH调整的P = 1.22 × 10 -6 )、类黄酮代谢( GO:0009812 ,OR = 28.73,BH调整的P = 3.98 × 10 -6 )以及抗坏血酸和醛糖代谢过程中发挥作用。(hsa00053,OR = 8.66,BH调整的P = 8.88 × 10 -4 ;Supplementary Fig. 22 and Supplementary Table 22 )。
根据GeneAnalytics注释,这些基因与多种疾病相关(例如,结直肠癌、乳腺癌、精神分裂症和神经系统疾病)(补充表23 )。我们发现影响BOD1 的CPC 特异性 AIS 由 71 个(61.2%)单倍体组装体携带。据报道,该基因与小脑运动功能障碍有关,,并且对人类认知功能至关重要。另一个受 AIS 影响的基因IL17RA ( n单倍型= 52;44.8%)在许多炎症和自身免疫性疾病中具有致病作用。特别是,该基因的多态性与东亚人群的特应性皮炎、自身免疫性 1 型糖尿病和哮喘有关。。 TWIST2 ( n 单倍型= 41, 35.3%)和CHFR ( n 单倍型= 34; 29.3%)在癌症转移中具有关键作用,并且是各种癌症的常用生物标志物。补充表24中列出了所有大小 >150 kb 的群体特异性 AIS,并且几乎所有 AIS 都仅存在于一个个体中( n单倍型= 1)。在所有研究的民族中,维吾尔族人口在 CPC 特定的古老基因渗入中所占的比例最大(0.51%,14.68 Mb)(补充表25 )。维吾尔族特有的 AIS 影响编码谷氨酰肽基转移酶的QPCT的一个值得注意的例子。该基因与欧洲和中国汉族人群的精神分裂症相关,也可能影响成年女性的骨矿物质密度,导致骨质疏松症的易感性。此外,我们发现了一种广为人知的致癌基因JUN ,它受到壮族人群中古老基因渗入的影响,这可能是造成JUN变异目前差异患病率及其与癌症相关性的原因。
我们发现,CPC 组装体中识别出的 AIS 中有 6.68% 受 SV 影响的基因,17.68% 的 SV 受到古老基因渗入的影响(补充表26 )。此外,在 CPC 组装体的 141 个 CNV 基因中检测到了 0.10% 的 AIS(补充表4 );特别是,在 135 个 CPC 特异性 CNV 基因中检测到了 0.09%。这些结果表明,CPC 数据对于增进我们对亚洲人类进化史的理解具有巨大潜力。
讨论先不看
Methods 方法
我们应用了一种基于主成分分析结果定量评估遗传多样性覆盖度的方法。我们使用统计量D d来选择个体来衡量总体样本的代表性,发现所选样本在前两个主成分图上靠近每个总体的聚类中心。另外,我们在各人群的样本选择上力求保持性别比例的平衡。经过 HiFi 数据的质量控制后,最终保留了 58 个样本用于本研究中的所有下游分析。
DNA sample and library preparation
QIAGEN 基因组试剂盒(目录号 13343,QIAGEN)进行纯化以进行常规测序(补充信息).这些淋巴母细胞样细胞系的生长受到限制,以最大限度地减少在扩展细胞培养过程中可能发生的遗传变化。在 1% 琼脂糖凝胶上监测 DNA 降解和污染。使用 NanoPhotometer 分光光度计 (IMPLEN) 检查 DNA 纯度。使用 Qubit DNA 检测试剂盒在 Qubit 2.0 荧光计 (LifeTechnologies) 中测量 DNA 浓度。
使用 SMRTbell Express 模板制备试剂盒 2.0 (Pacific Biosciences) 构建 15 kb 的文库。构建包括 DNA 剪切、损伤修复、末端修复、发夹接头连接、大小选择和文库纯化(补充信息)。
Hi-C sequencing
为了将杂交支架锚定到染色体上,从血液中提取了 Hi-C 文库的基因组 DNA。我们遵循前面描述的标准协议并进行了一些修改。简而言之,我们构建了 Hi-C 文库并通过 Illumina Novaseq-6000 平台获得了测序数据。将连接的 DNA 剪切成 300-500 bp 的片段,然后进行平末端修复和 A 尾,然后通过生物素-链霉亲和素介导的下拉进行纯化。最后,使用 Illumina Novaseq-6000 平台对 Hi-C 文库进行定量和测序。
NGS
使用每个样品总量为 1.5 μg DNA 作为 DNA 样品制备的起始材料。按照制造商的建议使用 Truseq Nano DNA HT 样品制备试剂盒 (Illumina) 生成测序文库,并将index codes 添加到每个样品的属性序列attribute sequences中。简而言之,通过超声处理将 DNA 样品片段化为 350 bp 的大小,然后对 DNA 片段进行末端抛光、连接A 尾并用全长接头连接,用于 Illumina 测序,并进一步进行 PCR 扩增。最后,纯化 PCR 产物(AMPure XP 系统),使用 Agilent 2100 生物分析仪分析文库的大小分布,并使用实时 PCR 进行定量。上述构建的这些文库通过 Illumina NovaSeq-6000 平台进行测序,并生成 150 bp 的双端读数,插入片段大小约为 350 bp。
HiFi sequencing
根据 PacBio 的标准方案 (Pacific Biosciences) 构建 SMRTbell 靶标大小文库,使用 15 kb 制备方案进行测序。简而言之,每个样品的总量为 15 μg DNA 用于 DNA 文库制备。基因组 DNA 样品根据文库的预期片段大小用 g-TUBE (Covaris) 剪切。然后去除单链突出端,对 DNA 片段进行损伤修复、末端修复和 A 尾。然后用发夹接头连接片段,用于 PacBio 测序。然后用 SMRTbell 酶纯化试剂盒用核酸酶处理文库,并用 AMPure PB 珠子纯化。通过 BluePippin (Sage Science) 筛选目标片段。然后用 AMPure PB 珠子纯化 SMRTbell 文库,并使用 Agilent 2100 生物分析仪 (Agilent Technologies) 检测文库片段的大小。在 PacBio Sequel II 或 Sequel IIe 平台 (Pacific Biosciences) 上使用单个 8M SMRT 芯片cell,在 Grandomics 中使用测序引物 V2 和 Sequel II 结合试剂盒 2.0 对文库进行测序。PacBio SMRT 分析套装 (https://www.pacb.com) 用于原始聚合酶读数的质量控制(补充信息)。
Oxford Nanopore Technologies sequencing
使用每个样品中总量为 3–4 μg DNA 作为 Oxford Nanopore Technologies 文库制备的起始材料。样品定性后,使用 PippinHT 系统 (Sage Science) 对长 DNA 片段进行大小选择(>50,000 和 >100,000)。接下来,修复 DNA 片段的末端,并使用 NEBNext Ultra II 末端修复/dA 加尾试剂盒(目录号 E7546)进行 A 连接反应。SQK-LSK109 (Oxford Nanopore Technologies) 中的接头用于进一步的连接反应,并通过 Qubit 4.0 荧光计 (Invitrogen) 测量 DNA 文库。
Genome assembly and quality control
在 2-5 个 SMRT cell上共对 68 个样品进行了测序。获得的子读数在 PacBio 工具中通过 ccs v6.3.0 转换为 HiFi 读数,使用 --hifi-kinetics --min-passes 3 --min-length 50。在个体水平上,我们结合了来自多个cell的 HiFi 读数,并应用 Hifiasm v0.16.1(参考文献 52)对 68 个 PacBio HiFi 样品进行初级组装和二倍体组装。
我们使用 QUAST v5.2.0 (参考文献 53) 来评估 68 个初级组装和 136 个单倍型组装的组装质量 (补充信息)。经过质量控制,我们去除了 3 个初级组装体 N50 < 20 Mb 或contig数 ≥ 2,000 的样品,以及 7 个具有两个单倍型组装体 N50 < 10 Mb 或contig数≥ 2,000 的样品。最后,保留了 58 个样品和 116 个高质量组件用于后续分析。
Assembly polishing and assessment
我们首先使用 Inspector v1.2(参考文献 54)来评估 116 个组装体的组装错误,这些组装体在 Inspector 的输出目录中记录为 bed 文件。然后,我们使用以下命令根据 Inspector 检测的组装错误评估结果进行组装polish:
inspector-correct.py -i asm.corrected/ --skip_structural -t 64。接下来,我们应用 Inspector 来评估校正前后的错装,以确定组装体polish是否提高了组装质量。使用 Inspector 根据每个组装单倍型中的组装误差计算来估计组装质量值。使用 Inspector 对组装体的组装错误评估命令如下:
inspector.py -c sample.ccs.fastq.gz -o asm.fa -r reference.gff.gz -o $asm/ --large --est-ref-size 31000000000 --no-icarus -t 64.
Variant identification from the phased assembly
使用预设参数 -ax map-hifi 通过 minimap2 将获得的 HiFi 读数与 T2T-CHM13 v2.0 参考基因组比对,然后通过 samtools sort进行排序。通过 samtools stats 评估reads的比对quality。BAM 文件用于识别常规和复杂的结构变异。
我们还使用以下命令应用 minimap2 v2.24(参考文献 55)将每个组装体与 T2T-CHM13 v2.0 参考对齐:
minimap2 -x reference.fa -m 10000 -z 10000, 50 -r 50000 --end-bonus=100 --secondary=no -a -t 20 -eqx -Y -O 5,56 -E 4,1 -B 5.
我们使用 dna-brnn v0.1(参考文献 56)来识别每个 CPC 组装体的高度重复序列。然后,我们使用Phased Assembly Variant (PAV) caller v1.2 (参考文献 17),它可以修剪prune由 minimap2 引起的重复比对,根据 minimap2 的输出 CIGAR 序列检测每个单倍型中的small variants和 SV。我们指定了参考基因组 (T2T-CHM13),每个样本的分型组装(two-phased assemblies 二倍体组装体)和配置文件作为启动 PAV 的输入。在我们的 PAV 流程中,合并了 200 bp 以内的 SV,每个样品的两个组装体之间有 50% 的重叠。
为了识别相对于T2T-CHM13 参考的CPC 样本的推定(putative)新序列,我们首先应用 SV-pop v3.0(参考文献 17)使用样本集合并和带有nr::szro(0.8,4):match(0.8) 参数的大小-倒数-重叠模型(size-reciprocal-overlap mode)来合并CPC核心样本的SV。在这种情况下,当插入之间的距离小于插入长度的四倍时,具有至少 80% 重叠的插入被合并。我们将插入定义为每个 SV 检出结果中相对于 T2T-CHM13 大于 1 kb 的新序列。我们应用 primatR v0.1.0(https://github.com/daewoooo/primatR) 使用函数 hotspotter 检测新序列的热点hotspots ,并设置参数 bw = 20,000 和 num.trial = 1,000。
Identification of gene duplication
我们使用 liftoff v1.6.3(参考文献 57)来检测每个组装中相对于 GRCh38 参考的重复基因。GENCODE v38 annotation注释文件用于构建liftoff的注释数据库。由 liftoff 注释到的额外基因拷贝数的蛋白质编码基因被鉴定为每个组装体的重复基因。我们使用参数 -sc 0.9以至少 90% 的置信度运行liftoff,命令如下:
liftoff -p 64 -sc 0.9 -copies -db GENCODE_V38.db -u $asm.unmapped -o asm.gff3 -polish asm.fa GRCh38.fa
我们使用 liftoff 输出 gff3 文件中的extra_copy_number 信息鉴定基因复制事件。
Liftoff 是一种基因注释转移工具,用于在参考基因组和目标基因组之间对基因注释进行映射。
Pangenome graph construction
我们应用 Minigraph-Cactus 管道构建了两个版本的 CPC 泛基因组变异图:基于 GRCh38 和基于 CHM13。我们首先使用 Minigraph (v0.19)60 构建了仅具有大于 50 bp 的 SV 的图形。从参考组装 GRCh38 或 CHM13 开始,作为初始图,剩余的一个参考基因组和 122 个组装依次mapping到图上。
我们使用 Minigraph-Cactus 管道将 SNP 级别的变异添加到图中。首先,鉴定着丝粒和端粒区域并用 dna-brnn软掩蔽,以排除高度重复序列对后续分析的影响。将组装体重新比对回 Minigraph,以在输入程序集的重叠群和 Minigraph 节点序列之间产生精确的对齐。