全文11,594字
，17 图，阅读 35 分钟。封面图源：
china.org.cn
2024 年 1 月 31 日，农历冬天的最后一个月，我们广州出生队列基因组学研究的文章——“
The Born in Guangzhou Cohort Study enables generational genetic discoveries
”[1]——正式在
Nature
主刊上发表了，这是广州出生队列第一个里程碑意义的成果，也是我第一篇作为第一作者发表到
Nature
主刊上的文章。这篇文章上线时 Nature 发了两个篇独立的报道[2,3]，对我们的成果做了评论，感兴趣的朋友可以查看文末的『原文链接』。
Nature 文章截图(第一作者：黄树嘉，刘斯洋，黄铭曦和何健荣；通讯作者：邱琇和夏慧敏)
这里我不打算重新对新闻报道的内容做赘述，我们队列的公众号也有
报道
。这里，我想做更专业、深入和完整的解读，讲讲这项研究的具体发现，以及项目过程中我的一些思考和体会。希望这可以更符合关注了我这个公众号的朋友们的口味，文章有些长，本想拆开成两篇的，想想还是算了，就做一篇吧，这样才连贯。
我喜欢从事生物信息学和生命科学研究，它既不枯燥又充满挑战。这个项目的进展并非一帆风顺，整个过程解决了很多问题，从我们讨论确定大致的研究路线到实际执行和最终发表成果，用了将近三年时间，其中审稿就花了一年半，来来回回审了 5 次（包括编辑的最后一次校审）！这是我截至目前为止历经过最漫长的审稿时间，中间我不止一次担心是不是要杯具了。
从整个文章的结构来说，归结起来，我们主要做了如下 4 件事：
-
完整构建了广州出生队列（BIGCS）自然人群的基础遗传图谱；
-
从群体遗传学的角度，把南方中国人（主要是粤语人群）的遗传历史讲得更清楚了；
-
母婴性状 GWAS 研究以及年龄特异性遗传效应的发现；
-
实现母婴遗传效应分离的方法并建立跨代孟德尔随机化分析方法。
研究背景
广州出生队列
广州出生队列（Born in Guangzhou Cohort Study, 简称 BIGCS）从 2012 年启动至今已经整整做了 12 年——应该是我们国内较早的一批自然人群出生队列，取得了很多好成绩，这些已不需我来重述。但我们一直希望能够在基因组学研究上有所突破，流行病学研究并非总能有效回答在疾病和性状上所发生的现象，我们希望可以从更加根本的遗传特征上去寻找关键信息。但这得一步一步来，BIGCS 第一步是先从入组人群中随机抽取了 2000 个健康的母婴对(其中原计划是包含 400 个家系)，总共是 4,215 个人（我必须说这个样本并不算大），考虑到测序成本的原因，只做了 7x 的全基因组测序（WGS）。2020 年 6 月我加入了广妇儿，来到了出生队列研究室的团队里，并从那个时候开始正式投入到这个项目中。在这个过程中，我首先搭建了整个基因组学分析框架和流程以及项目的初步研究方案（与技术细节），并和团队共同讨论最终确定了整个项目的故事线和方向。
这个项目采用的是大样本+低测序深度的策略，对自然人群的组学研究来说，我会认为样本量对于科学发现而言是更加重要且更具有决定性的因素。所以在同样的件下，应该检测尽可能多的样本，虽然如此我们的第一阶段也只能做到 4,215 例。再加上每个样本只测了7x——最终原始数据在经过我质控之后平均深度就只剩下 6.63x 了——这个深度应该是参考了千人基因组计划的策略，我同时做了一些数学上的推导分析了这个深度的合理性以及它的潜在挑战，如下图：
图 1. 测序深度与基因组覆盖率
这里的 图1.a 是不同测序深度条件下，对基因组的覆盖率，可以看到 6.6x 的深度理论上其实已经能够覆盖整个基因组 99.86% 的区域了，10x 以上时更是无限接近 100%，所以从数据的覆盖度来说这个深度是 OK 的。但光看这个是不够的，因为我们是要进行变异检测的，最终分析的落脚点是在突变数据集上。
对于突变数据来说，我们要准确全面地检测出来，就得确保在两个同源染色体（人是二倍体，有来自父母的两个同源染色体）上能测到足够的 DNA 序列，否则就可能漏掉信息导致变异检测不准确——特别是对于杂合突变来说更是如此。所以，我分析了在不同的测序深度条件下，基因组上任意一个位点都能被
测自两个同源染色体
的 read 覆盖 N 次的概率（图1.b）。我觉得每一个位点起码要被两个同源染色体覆盖分别两次才可以在变异检测的时候得到较准确的结果。在图里我们可以看到，对于 6.6x 的深度来说，理论上只有 71% 的位点符合这个要求（红色曲线），应该说对于基因组的纯合突变来说问题不大，但对于杂合突变的检测而言，就可能有 ~30% 的杂合突变位点会出现
杂合丢失
的假阴/假阳性结果（当然我的假设条件较严苛，实际比例应该不至于如此高），这个也是我后面必须想办法解决的地方。
当然理论分析也仅仅只是数学计算，最终情况还是要看真实数据的表现的。但理论分析的好处就是，可以让我们心里有一个衡量标准，而不是两眼一抹黑，然后瞎搞。
我们都做了哪些研究
那么，接下来我将逐一解释我们主要都做了哪些事情，以及为什么要这样去做。
1. 构造分析流程，完成变异检测，构建 BIGCS 自然人群的基础遗传图谱
第一步自然是要完成所有样本的全基因组变异检测，获得准确的突变数据集，这是实现所有下游分析的关键。
刚开始时我们没有任何关于基因组数据分析流程，团队除了我之外其他人也没有足够的经验。如果只是为了完成一两个人的全基因组数据分析，其实也不麻烦，参考 GATK 的文档也能做下来，但我们是四千多个样本，中间分析步骤繁多，产生的数据文件更是成千上万，手动检查是不可能的，所以除了核心步骤之外，要完成这个事情就还得具备对任务完成情况进行监测的能力。为此，我从头设计并实现了一个适合于大规模 WGS 数据分析的半自动化流程（如下图 2）——ilus，我在公众号里也分享过（
ilus：一个轻量级全基因组（WGS）和全外显子（WES）最佳实践分析流程生成器
)
。有了这个之后，即便只是一个人也可以轻松搞定这个分析。
图 2. 主分析流程（https://github.com/ShujiaHuang/ilus）
*目前 ilus 已经到了版本 2.0，除了 GATK 最佳实践之外还兼容了 Sentieon，并且可以进行 WES 或者其他捕获型测序策略的数据分析流程的构建。如果可能的话，我特别推荐你用 Sentieon。
过程虽然繁琐但就都比较常规，我就不再赘述，Nature 文章的方法部分对此也做了很详细的描述。但在对突变数据集进行深入质控方面，
我想特别指出的是一定不要只做一个 VQSR 就完事了——我的这个判断应该可以适合大多数 WGS 项目，对于低深度 WGS 更是如此
。仅仅只是 VQSR 得到的结果和芯片数据比起来质量是不行的——我们的数据显示 SNP 还不足 95% 的一致性。
原因其实我也知道，主要还是来源于上面理论分析时提到的杂合丢失，如何解决这个问题呢？以前在用
NIPT 超低深度数据进行 14 万人的遗传学研究时
，我就已经知道，这个时候进行 Imputation 将是一个合理的选择。但是，这里不应借助 reference panel，否则数据会被局限在 panel 之下，这对我们是不但得不偿失，而且不正确。所以，我自然地想到了通过 reference panel free 的方式来进行——也就是不用 panel 了，仅仅依靠自己数据的情况来进行推断。但要考虑到深度低，每个样本此刻的基因型会具有一定程度的不确定性，因此在重新填补和推断的时候应该用基因型似然概率（genotype likelihood）而不是直接用 genotype，这也就是我最终的解决路线，这个步骤我们称之为
LD-Based refinement
，因为这其实是通过连锁不平衡来实现内部基因型填补和校正的。那么最后效果如何呢？可以说非常不错，如下图 3：
图 3. 变异质控
可以看到 refinement 之后，基因型的整体一致性就达到了 99% 以上，这就符合预期了。当然，这个过程中还有其他细节，文章里也有详细说明，但比较细枝末节，我就不在这里重述。
最后，我再提一嘴，我们原本的样本一共是 4,215 例，做完初步变异检测之后（还未质控），我们对样本的测序质量、性别一致性、亲缘关系一致性等等各个方面都做了层层过滤，最终就只剩下了正文里说到的 4,053 例（含 332 个三人家系，1,406 个亲子对和 245 个无亲缘关系个体）。
其实，上面定下的突变数据过滤条件十分严格，被过滤掉的数据很多，大约 33% 的原始突变数据被过掉了(最后剩下 56M 个群体突变位点)——是有些可惜的，但我也不想得到假的东西。有了这个干净的突变数据集之后，我们还用它构造了一个新的特别适合于中国人群的参考图谱（reference panel）可以为研究人员提供更好的基因填补（Imputation），并将它做成了在线的数据库，名字是：Genome Database of Born In Guangzhou Cohort study (简称：GDBIG，
http://gdbig.bigcs.com.cn/
)，对于有需要的朋友直接到这个网站注册就可以使用了。这里我吸收了之前我在华大构建 Chinese Millionome Database（CMDB）的不足，把这个 GDBIG 做得更好了，你如果用了应该可以体会到，功能也更多（包括，突变频率查询、基因型填补和特定母婴表型的 GWAS meta-analysis 等），还可以通过 API 在命令行下查询。
GDBIG 主页
2. 群体遗传学分析进一步揭开南方中国人的遗传历史
对于广州出生队列的亲子对基因组数据，原本我们并不需要将群体遗传学作为重点内容来凸显我们的特色。我们之所以还是要这么做有两个原因：
-
第一，项目刚开始时，我也想将焦点集中在 GWAS 上，但我做了功效分析（Power analysis）之后，很快就意识到风险太大了，样本量不足，恐怕只能发现一些别人发现过的信息（而且还不一定可以看全）。初看我们有四千多个样本，似乎不少，但实际上如果要做 GWAS，母子需要分组，这样一来，母亲组和小孩组分别就都不到两千人了，结果会如何我没有信心。所以，最好是能有其他的突破口，这个时候人群的群体遗传研究就成了一个重要的选择，而且这恰好是南方中国人群，此前研究并没有很好覆盖到，那么做得好的话是有可能有新发现的，填补一些空白也将是有意义的。
-
第二，作为潮汕人，我其实非常好奇南方中国人的群体遗传特征，特别是广东省三个主要方言人群之间是否有什么有趣的特点。我们的随访数据也有很好的方言信息，刚好适合做这个事情。
所以，群体遗传学分析自然就成为我们重点要探索的一个方面。
我们也的确有了一些新的发现，首先是在 PCA 分析中，我们发现遗传信息与方言的匹配度要显著高于籍贯和民族，这其实并不让人感到意外，也许你不做这个分析也会觉得应该如此，
但我们用数据给出了直接的证明
。而且，这应该是第一回用人群全基因组数据把中国人中 10 种主要的方言与他们的遗传信息建立了联系（图 4）。眼尖的朋友应该还隐隐可以看出，这个图谱和方言人群在地理上的分布十分相近。我想这对于中国传统方言的研究将是有积极意义的。
有了上面的分析之后，我们就想往历史深处看得更远一些。我们十分好奇，是否还可以从遗传学角度推断出更多和这些方言有关的人口统计学历史（demographic history）？看看我们的方言祖先们在历史的长河中是否发生过什么有趣的故事。那么，说干就干，考虑到样本中方言人群的样本量问题，我们就以五个样本量大于 50 的方言群体（普通话、湘语、粤语、客家话和闽南语）作为代表开展研究。
经过一番分析之后，
我们清楚地看到这五个方言群体之间的人口动态变化是相似的，都大约在 15 万年前进入人口瓶颈期，有效人口数不断减少，并大约在 2.5 万年前达到最低点，随后才开始逐步恢复（图 5a）
。为什么会这样？我查找了该时段的地球地质历史，发现该时期恰好是地球的末次冰期，并且在大约 2.65 万至 2 万年前达到冰期顶点——即
末次盛冰期
。在这一时期中，地球上的冰川规模达到最大，极端严寒和干燥的气候导致生存环境变得极其恶劣，有效人口的大幅减少就发生在这一时期（从全球其他人群的基因组上同样可以观察到类似情形），之后才开始逐渐恢复。
另外，我们还计算了这五种方言群体之间的遗传分化时间，发现他/她们的分化始于新石器时期（史前一万年），并在大约 3千-4 千年前逐渐稳定下来，而这个时间恰好与中国古代农业社会体系的形成和发展相吻合（图 5b）。
图 5. 通过 BIGCS 基因组数据推断中国人的人口统计学历史
最后，既然粤语占优，我们当然希望再看看这个群体的其他特征。首先是结合古人 DNA 做了 ADMIXTURE 分析，我们可以非常直观地看到，
生活于南方的粤语人群主要是中国北方人（包括黄河流域一带的北方古人，图中橙金色)与南方壮侗语系人群（Tai-Kadai，图中浅蓝色）的混合体，并且有更明显的北方古人痕迹，而南方古人的成分则很低（图中深蓝色）——也就是说粤语人的祖先应该是北方古人（图 6）
，当然做出这个判断并非只靠 ADMIXTURE 分析，而是还进行了其他的分析，特别是 f4 statistics 结果的支持（这里要特别感谢付巧妹老师团队在这一块的分析和解读）。
图 6. ADMIXTURE
既然我们看到粤语群体的祖先大概率来自北方，那么一个顺理成章的问题就是，他们是什么时候南下的？我们通过计算人群开始发生遗传混合的时间点来推断南下时间（图 7）。
结果我们发现，这个时间大约是 1,464 年前。但要注意这个时间要从出生时算起，按照 BIGCS 的年龄结构，要往前推大约 30 年，考虑到样本入组和测序时间，我们从 1990 年算起应该比较合适。所以，粤语人群和南方人发生混合的时间大概是公元 1990 - 1464 = 526 年（南北朝时期）。但是，这已经是和当地人发生基因混合的时间点，人群迁入时间则还要再往前推一些，可能是始于五胡乱华、中原汉人被迫“永嘉南渡”之时（图 8），这是一场前无古人的大迁徙，这个时期西晋朝廷渡过长江，建立东晋政权，此时同宗同族的人们抱团取暖，社会强调血统与身份，家族传承形成书面谱系。东晋朝廷依据父系家谱，判定官员的品阶甚至发展出“谱学”——“上品无寒门，下品无势族”（出自《晋书·卷四十五·刘毅传》）。
