华大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。

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

  • 发表于:
  • 原文链接https://kuaibao.qq.com/s/20181012G24XVC00?refer=cp_1026
  • 腾讯「云+社区」是腾讯内容开放平台帐号(企鹅号)传播渠道之一,根据《腾讯内容开放平台服务协议》转载发布内容。
  • 如有侵权,请联系 yunjia_community@tencent.com 删除。

同媒体快讯

扫码关注云+社区

领取腾讯云代金券