华大cell文献学习笔记

快，关注这个公众号，一起涨姿势～

近日，华大基因在Cell上发表中国人基因组大数据研究成果引起行业内的关注。

cell官网截图

我昨天也阅读了这篇网红文献，按照自己平时读文献的惯例，整理了文献学习笔记。在内部分享时，有同事建议我将这篇笔记放到公众号上。

由于个人的水平有限，有对原文理解不正确的地方，还请各位老师指出。欢迎留言讨论，我们大家一起学习。

文/徐大梨

一、主要工作及亮点

该文章的工作内容是利用14万中国孕妇的NIPT数据完成大规模的基因组学研究，包括群体遗传学，多表型GWAS，血浆病毒谱三个主要方面。

工作亮点：

1. 利用超低深度（0.06x-0.1x）NIPT数据，开拓了一个全新的组学大数据研究思路，证明了低深度测序数据对于大规模的组学研究来说是可行的。与那些经费昂贵的大规模基因组学计划相比，数据挖掘的性价比高。

2.揭示中国人群的遗传结构和特点。（1）群体变异检测，得到了904万个高质量SNP（包含20万个新发现的SNP），建立中国人群等位基因频率数据库CMDB。（2）揭示人群遗传差异，包括：a. 汉族与36个少数民族的遗传距离；b. 31个省人群与欧洲、南亚以及东亚人群的基因交流程度；c. 找到中国南北方人群存在遗传差异的6大基因。

3.用NIPT数据作GWAS研究。针对4个表型：身高、BMI、怀孕年龄、双胎怀孕，发现了若干新关联位点，填补了相关研究的空白。（1）发现并且验证了48个与身高以及13个与BMI显著相关的基因位点，可解释48%的身高遗传率以及10%的BMI遗传率；有7个与身高以及3个与BMI相关的基因位点属于首次发现。（2）发现了2个与怀孕年龄显著相关的基因位点。（3）发现了1个和双胞胎妊娠显著相关的突变位点。

4.血浆病毒谱分析。（1）展现了全国各省人群病毒感染发生率和病毒在血浆中丰度的分布。（2）发现中国人血浆的病毒组与欧洲人相比差异较大，乙肝病毒是第1位。（3）将病毒感染作为表型，GWAS找到了1个与HHV-6易感性极显著相关的基因突变。

二、数据和方法

我是做生信的，所以重点关注方法。原文方法部分的内容较多，只选取了个人感兴趣的部分做了笔记。

1.样本和数据信息

入组样本：所有受试者均是2012-2013期间在BGI抽取外周血，接受NIPT检测(Zhang et al., 2015)。

一些人口学信息的获取：

（1）地理位置信息：以匿名方式，通过身份证前6位（示受试者的出生地）得到。点评：这个方式值得我们在项目中借鉴。

（2）出生年份和年龄信息：从知情同意书上得到，28和35是两个峰值。

（3）身高和BMI：身高和体重是在医院取血时候进行测量的；BMI=体重（kg）/身高^2（m^2）。

（4）胎儿浓度：通过男胎样本的chrY的信息推断。总体胎儿浓度在3.5%-30%之间，均数8%。

出组样本：剔除测序错误率大于1%的样本（3278）；剔除检测到染色体非整倍性异常或者父亲报告核型异常的样本（502）。

纳入分析141, 431人，其中包括测序为35bp read的118, 576人；测序为49bp read的22, 855人。

2. 测序数据质控

（1）用SOAPnuke过滤低质量的reads，过滤标准为含有30%的低质量碱基（Q≤2）或者N。

通常1个样本含有5-10M clean reads，覆盖度为0.06x-0.1x。

3. 数据比对和选取accessible region

比对的步骤：

（1）将clean read用bwa比对到hg19；

（2）用samtools rmdup去除可能由PCR带来的dup

（3）用GATK中的modules重新比对indel区域的reads和校正碱基的质量值。

（4）最终的比对文件存储为标准的CRAM格式。

比对文件QC的步骤：

（1）用samtools的stats功能来统计样本的测序错误率。中位数为0.3%，剔除测序错误率大于1%的样本。

（2）用samtools depth评估比对质量值大于30的read和碱基质量值大于20的碱基的覆盖度。比对质量值低和碱基质量值低的不纳入后续分析。

考虑read太短，比对可能存在的问题，在hg19上定义了accessible region，只在这个region范围里做后续的call 变异，群体分析，关联分析。

accessible region的选取：

（1）剔除Heng Li给的35-kmer problematic alignment bed文件中的区域，这些是35bp kmer难以unique map的区域。剔除897,319,085 bp。

（2）保留UCSC的wgEncodeDukeMapabilityUniqueness35bp.bedGraph文件中mappability uniqueness score为1的区域。剔除6,459,019 bp。

（3）保留测序深度在3000到30,000之间的区域。再剔除4,340,936 bp。

最终hg19上的accessible region包含22条常染色体和chrX，范围是2,128,184,806 bp.。

点评：QC做得很细致，这个accessible region的选取方法很值得我们借鉴。

这个accessible region更针对35bp的样本。我很好奇，他们对49bp的样本也用同样的accessible region吗？文章没提到。

4.群体变异检测和估计等位基因频率

（1）使用了一个最大似然的方法（S. Liu, unpublished data）来识别多态性位点和推断等位基因频率。文章QUANTIFICATION AND STATISTICAL ANALYSIS部分有该方法的数学细节。

（2）这里值得一提的是用到了一个single-read sampling的策略，即在某个多态性位点上，同一个样本中如果有多条read都覆盖该位点，只能选取其中的1条。这样能够保证覆盖这个多态性位点的全部read都来于自不同的人，目的是消除样本胎儿浓度的不同造成的影响。点评：思路好棒！

5. 40位受试者的中等深度测序

（1）对其中40位受试者用Hiseq X10做15X测序，用bwa-mem比对到hg19，再用GATK multi-sample best practice对这40个样本call SNP并得到genotype。

目的：用这些结果作为超低深度NIPT数据call SNP和impute genotype的基准参照。

6. Imputation genotype

（1）用STITCH（v1.127）对14万样本来impute genotype probabilities。

相关参数和细节：

最关键参数祖先单倍型的数目K选取为10，这个值是基于在chr3中5M区域的测试以及和STITCH的作者讨论得到。

在EM算法的迭代中，AF初始值来自于1KGP中中国人群(CHB+CHS+CDX, n = 301)的AF。

Imputed Loci的选取来自于1KGP中东亚人群的22条常染色体和chrX上AF>0.01的位点。

点评：虽然无法得到低深度数据的准确SNP，但文章用imputed genotype probabilities作为GWAS中的基因型，一下子解决了！

（2）为了保证impute的准确率，对于每个impute出的位点，计算IMPUTE2-style info score和Hardy Weinberg equilibrium的p值。

保留的标准：IMPUTE2-style info score>0.4且Hardy Weinberg equilibrium的p值大于10的负6次方。

点评：从40个中等测序样本的genotype和这里impute的genotype的相关性的R方来看，impute的准确性还是靠谱的。

8. 主成分分析

（1）计算任意两个个体间的covariance

选取的位点是已知的，且AF>0.05，且这里依旧用了single-read sampling策略。

点评：用于做主成分分析的变异位点的选取上，我要给个差评！

一般合理的做法是先thin位点，保证位点彼此间隔至少为0.001cM，然后删除掉主要组织相容性复合体（MHC）区域，乳糖酶区域（2q21），倒位区域8p23 和17q21.31 上的位点。

文章并没有提到这些，个人猜测是不是accessible region里已经去除了这些区域？

（2）用R里的Spectra包对人群的covariance matrix做decomposition，就能得到主成分。

发现前三个主成分中，第一个主成分对应read长度，第三个主成分对应测序错误率。考虑这一情况，只选取了35bp测序和测序错误率低于0.00325的96,880 个样本重新做进一步的主成分分析。

对汉族做主成分分析的时候，还剔除了民族信息缺失的样本，只用了45,387个样本。

（3）利用前几个主成分能够进一步计算人群的遗传距离。

9. PAF的分析

Private allele指的是在一个人群中是多态性的位点，但是在其它人群中的固定的位点。在1KGP的数据中，CEU和ITU人群分别带有3, 485, 371 和 4, 324,376 private alleles。用MAF>0.05来筛选，分别剩余66,700 and 45,536个 private alleles.

对于每一个省的数据，可以计算其与1KGP中人群K的PAF，公式如下

这里，M为各省数据中带有人群K的private alleles的位点数目。Ns是位点s上所有等位基因出现的次数，ns是位点s上K人群私有的等位基因出现次数。

PAF越高，两个人群的基因交流程度越高。

10.Clinvar的致病性位点上，等位基因频率的差异

（1）根据纬度划分的三个地理区域的人群，分别统计等位基因频率。

（2）选取了3,238个 bi-allelic 可能致病的位点，证据等级是5。

（3）对于南部，北部和中部人群，分别计算risk allele和non risk allele B的数目。点评：在这一步给个赞，并没有简单对两种allele进行计数，给出的公式里还考虑了测序质量的权重。

（4）对于每个位点，都能得到一个2X3的表，2代表2个等位基因，3代表3个区域。对该表做fisher检验，得到p value。

11. GWAS分析

对imputed genotype probabilities和表型的关联信号进行检测。

在回归分析的协变量（covariate）选取方面：

（1）对于身高表型，用到的协变量：前五个主成分，孕妇年龄，胎儿性别。

（2）对于BMI表型，协变量里还增加了孕周。

（3）对于双胎和怀孕年龄的表型，所用的协变量与身高表型一样。

点评1：此前，华大基因公众号推文的官方解读是“在NRG1基因中发现了一个和双胞胎妊娠显著相关的突变位点，也就说携带有NRG1基因特定突变的孕妇，有更高的几率怀上双胞胎。”

个人对此保留质疑。

因为双胎在很大程度上与试管婴儿（IVF）有关。个人曾处理过NIPT数据，回顾某阶段的数据，共1663个IVF样本，其中有545个为双胎，也就是说IVF样本中双胎比例大于30%，远远高于普通人群。因此IVF信息很有必要在GWAS分析中被列为协变量。

由于文章的协变量部分并未提及包含IVF的信息，我很怀疑双胎妊娠的关联信号是false positive。甚至故事有可能是相反的，NRG1上该位点的突变，导致甲状腺亢进，进而导致怀孕困难，导致患者选取去做IVF，因而增加了双胎的概率。

当然这些都是我的猜测，并没有证据。NRG1基因是否是双胎基因还需要后续的验证。

点评2：双胎妊娠的表型还可以进一步细分，因为同卵和异卵的生物学机制不同。譬如可以根据常染色体胎儿浓度定量新方法，与chrY得到的胎儿浓度相比较，找出其中的龙凤胎表型（龙凤胎一定是异卵），进一步做GWAS。

虽然提出一些质疑，但是瑕不掩瑜。我认为这篇文章真的特别棒，让我从中学习了很多，得到不少启发。