WES分析七步走

WES——全外显子测序，已经逐步应用到医疗领域，本例的测试数据就是一个自闭症家系的全外显子测序数据结果，我本来雄心勃勃的想帮助他解决病因的问题，后来发现也只是会一点数据分析而已。该系列文章最初写于2015年11月。

目录如下：

  • 测序质量控制

  • snp-calling

  • snp-filter

  • 不同个体的比较

  • 不同软件比较

  • annovar注释

  • de novo变异情况

其实以上还远远不够，剩余的分析在直播我的基因组系列会涉及到，欢迎大家继续围观！

测序质量控制

这一步主要看看这些外显子测序数据的测序质量如何：

首先用fastqc处理，会出一些图表，通常是不会有问题的，毕竟公司不想砸自己的牌子。(同时推荐multiqc这个python程序对多样本的项目的fastqc结果输出动态可交互式html报告)

然后粗略统计下平均测序深度及目标区域覆盖度，这个是重点，不过一般没问题的，因为现在芯片捕获技术非常成熟了，而且实验水平大幅提升，没有以前那么多的问题了。

这个外显子项目的测序文件中mpileup文件是1,371,416,525行，意味着总的测序长度是1.3G，以前我接触的一般是600M左右的。(说明这个项目的测序数据量比较足)因为外显子目标区域并不大，就34,729,283bp，也就是约35M，即使加上侧翼长度。

  1. 54692160：外显子加上前后50bp

  2. 73066288：外显子加上前后100bp

  3. 90362533：外显子加上前后150bp

然后我要根据外显子记录文件对mpileup文件进行计数，统计外显子的coverage还有测序深度，这个脚本其实蛮有难度的。

我前面提到过外显子组的序列仅占全基因组序列的1%左右，而我在NCBI里面拿到 consensus coding sequence (CCDS)记录 CCDS.20150512.txt文件，是基于hg38版本的，需要首先转换成hg19才可以来计算这次测序项目的覆盖度和平均测序深度。(PS:其实可以直接下载hg19版本的)

参考：http://www.bio-info-trainee.com/?p=990 （ liftover基因组版本之间的coordinate转换）

  1. awk '{print "chr"$3,$4,$5,$1,0,$2,$4,$5,"255,0,0"}' CCDS.20150512.exon.txt >CCDS.20150512.exon.hg38.bed

  2. ~/bio-soft/liftover/liftOver CCDS.20150512.exon.hg38.bed ~/bio-soft/liftover/hg38ToHg19.over.chain CCDS.20150512.exon.hg19.bed unmap

下面这个程序就是读取转换好的外显子记录的数据，对一家三口一起统计，然后再读取每个样本的20G左右的mpileup文件进行统计，所以很耗费时间。

统计结果如下表，外显子目标区域平均测序深度接近100X，所以很明显是非常好的捕获效率啦！而全基因组背景深度才3.3，这符合实验原理，即与探针杂交碱基多的片段比少的片段更易被捕获。对非特异杂交的基因组覆盖度非特异的背景 DNA 也进行了测序。

接下来对测序深度进行简单统计，脚本如下，但是这个图没多大意思因为我们的外显子的35M区域平均都接近100X的测序量。

这个其实没必要自己写脚本啦，bedtools的作者给出了WES的QC图的制作方式，简单的谷歌就可以查到，而且非常快速和方便。

snp-calling

准备文件：下载必备的软件和参考基因组数据

软件

ps：还有samtools，freebayes和varscan软件，我以前下载过，这次就没有再弄了，但是下面会用到

参考基因组

参考突变数据

1.下载数据

2.bwa比对

3.sam转为bam，并sort好

4.标记PCR重复，并去除

5.产生需要重排的坐标记录

6.根据重排记录文件把比对结果重新比对

7.把最终的bam文件转为mpileup文件

8.用bcftools 来call snp

9.用freebayes来call snp

10.用gatk来call snp

11.用varscan来call snp

snp-filter

其中freebayes,bcftools,gatk都是把所有的snp细节都call出来了，可以看到下面这些软件的结果有的高达一百多万个snp，而一般文献都说外显子组测序可鉴定约8万个变异。

这样得到突变太多了，所以需要过滤。这里过滤的统一标准都是qual大于20，测序深度大于10。过滤之后的snp数量如下

  1. perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if /INDEL/;/(DP4=.*?);/;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$1"}' Sample3.bcftools.vcf >Sample3.bcftools.vcf.filter

  2. perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if /INDEL/;/(DP4=.*?);/;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$1"}' Sample4.bcftools.vcf >Sample4.bcftools.vcf.filter

  3. perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if /INDEL/;/(DP4=.*?);/;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$1"}' Sample5.bcftools.vcf >Sample5.bcftools.vcf.filter

  4. perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next unless /TYPE=snp/;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]"}'  Sample3.freebayes.vcf > Sample3.freebayes.vcf.filter

  5. perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next unless /TYPE=snp/;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]"}'  Sample4.freebayes.vcf > Sample4.freebayes.vcf.filter

  6. perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next unless /TYPE=snp/;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]"}'  Sample5.freebayes.vcf > Sample5.freebayes.vcf.filter

  7. perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if length($F[3]) >1;next if length($F[4]) >1;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]:$tmp[2]"}'  Sample3.gatk.UG.vcf  >Sample3.gatk.UG.vcf.filter

  8. perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if length($F[3]) >1;next if length($F[4]) >1;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]:$tmp[2]"}'  Sample4.gatk.UG.vcf  >Sample4.gatk.UG.vcf.filter

  9. perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if length($F[3]) >1;next if length($F[4]) >1;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]:$tmp[2]"}'  Sample5.gatk.UG.vcf  >Sample5.gatk.UG.vcf.filter

  10. perl -alne '{@tmp=split/:/,$F[9];next if $tmp[3]<10;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[3]"}' Sample3.varscan.snp.vcf >Sample3.varscan.snp.vcf.filter

  11. perl -alne '{@tmp=split/:/,$F[9];next if $tmp[3]<10;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[3]"}' Sample4.varscan.snp.vcf >Sample4.varscan.snp.vcf.filter

  12. perl -alne '{@tmp=split/:/,$F[9];next if $tmp[3]<10;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[3]"}' Sample5.varscan.snp.vcf >Sample5.varscan.snp.vcf.filter

这样不同工具产生的snp记录数就比较整齐了，我们先比较四种不同工具的call snp的情况，然后再比较三个人的区别。

然后写了一个程序把所有的snp合并起来比较

得到了一个很有趣的表格，我放在excel里面看了看 ，主要是要看生物学意义，但是我的生物学知识好多都忘了，得重新学习了！（其实这里的miss位点并不能确定是野生型，也有可能是该位点没有被芯片捕获，需要修正代码）

不同个体的比较

3-4-5分别就是孩子、父亲、母亲

我对每个个体取他们的四种软件的公共snp来进行分析，并且只分析基因型，看看是否符合孟德尔遗传定律。

结果如下：

粗略看起来好像很少不符合孟德尔遗传定律，然后我写了程序计算。（这里面有一个问题，就是我把VCF没有记录的位点，都当做是野生型，这个其实是错的，所以这样简单粗暴的脚本查找de novo变异，是不可取的）

总共127138个可以计算的位点，共有18063个位点不符合孟德尔遗传定律。我检查了一下不符合的原因，发现我把

chr1 100617887 C T:DP4=0,0,36,3 T:1/1:40 T:1/1:0,40:40 miss T:DP4=0,0,49,9 T:1/1:59 T:1/1:0,58:59 miss T:DP4=0,0,43,8 T:1/1:53 T:1/1:0,53:53 T:1/1:50

计算成了 chr1 100617887 C 0/0 0/0 1/1 所以认为不符合，因为我认为只有四个软件都认为是snp的我才当作是snp的基因型，否则都是0/0。那么我就改写了程序，全部用gatk结果来计算。这次可以计算的snp有个176036，不符合的有20309，而且我看了不符合的snp的染色体分布，Y染色体有点异常。

但是很失败，没什么发现。

不同软件比较

主要是画韦恩图看看，参考：http://www.bio-info-trainee.com/?p=893

对合并而且过滤的高质量snp信息来看看四种不同的snp calling软件的差异

我们用R语言来画韦恩图代码如下：

由下图可以看出不同软件的差异还是蛮大的，所以我只选四个软件的公共snp来进行分析

首先是sample3

然后是sample4

然后是sample5

可以看出，不同的软件差异还是蛮大的，所以我重新比较了一下，这次只比较，它们不同的软件在exon位点上面的snp的差异，毕竟，我们这次是外显子测序，重点应该是外显子snp

然后我们用同样的程序，画韦恩图，这次能明显看出来了，大部分的snp位点都至少有两到三个软件支持

所以，只有测序深度达到一定级别，用什么软件来做snp-calling其实影响并不大。

用annovar注释

使用annovar软件参考自：http://www.bio-info-trainee.com/?p=641

  1. /home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4  Sample3.varscan.snp.vcf > Sample3.annovar

  2. /home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4  Sample4.varscan.snp.vcf > Sample4.annovar

  3. /home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4  Sample5.varscan.snp.vcf > Sample5.annovar

然后用下面这个脚本批量注释

Reading gene annotation from /home/jmzeng/bio-soft/annovar/humandb/hg19_refGene.txt ... Done with 50914 transcripts (including 11516 without coding sequence annotation) for 26271 unique genes

最后查看结果可知，真正在外显子上面的突变并不多

  1. 23515 Sample3.anno.exonic_variant_function

  2. 23913 Sample4.anno.exonic_variant_function

  3. 24009 Sample5.anno.exonic_variant_function

annovar软件就是把我们得到的十万多个snp分类了，看看这些snp分别是基因的哪些位置，是否引起蛋白突变，位置信息如下：

  1. downstream

  2. exonic

  3. exonic;splicing

  4. intergenic

  5. intronic

  6. ncRNA_exonic

  7. ncRNA_intronic

  8. ncRNA_splicing

  9. ncRNA_UTR3

  10. ncRNA_UTR5

  11. splicing

  12. upstream

  13. upstream;downstream

  14. UTR3

  15. UTR5

  16. UTR5;UTR3

de novo变异情况

de novo变异寻找也属于snp-calling的一部分，但是有点不同的就是该软件考虑了一家三口的测序文件，找de novo突变。

功能介绍：http://varscan.sourceforge.net/trio-calling-de-novo-mutations.html

而且还专门有一篇文章讲ASD和autism与de novo变异的关系，但个人感觉文章不清不楚，没什么意思

Trio Calling for de novo Mutations

  1. Min coverage:   10

  2. Min reads2:     4

  3. Min var freq:   0.2

  4. Min avg qual:   15

  5. P-value thresh: 0.05

  6. Adj. min reads2:        2

  7. Adj. var freq:  0.05

  8. Adj. p-value:   0.15

Reading input from trio.filter.mpileup

  1. 1371416525 bases in pileup file （137M的序列）

  2. 83123183 met the coverage requirement of 10 （其中有83M的测序深度大于10X）

  3. 145104 variant positions (132268 SNP, 12836 indel) （共发现15.5万的变异位点）

  4. 4403 were failed by the strand-filter

  5. 139153 variant positions reported (126762 SNP, 12391 indel)

  6. 502 de novo mutations reported (376 SNP, 126 indel) （真正属于 de novo mutations只有502个）

  7. 1734 initially DeNovo were re-called Germline

  8. 12 initially DeNovo were re-called MIE

  9. 3 initially DeNovo were re-called MultAlleles

  10. 522 initially MIE were re-called Germline

  11. 1 initially MIE were re-called MultAlleles

  12. 3851 initially Untransmitted were re-called Germline

然后我看了看输出的文件trio.mpileup.output.snp.vcf

软件是这样解释的

The output of the trio subcommand is a single VCF in which all variants are classified as germline (transmitted or untransmitted), de novo, or MIE.

  • FILTER - mendelError if MIE, otherwise PASS

  • STATUS - 1=untransmitted, 2=transmitted, 3=denovo, 4=MIE

  • DENOVO - if present, indicates a high-confidence de novo mutation call

里面的信息量还是不清楚。

我首先对拿到的trio.de_novo.mutaion.snp.vcf文件进行简化，只看基因型！

  1. head status.txt   （顺序是dad,mom,child）

  2. STATUS=2 0/0 0/1 0/1

  3. STATUS=2 1/1 1/1 1/1

  4. STATUS=2 0/1 0/0 0/1

  5. STATUS=2 1/1 1/1 1/1

  6. STATUS=1 0/1 0/0 0/0

  7. STATUS=1 0/1 0/0 0/0

  8. STATUS=2 1/1 1/1 1/1

  9. STATUS=2 1/1 1/1 1/1

  10. STATUS=2 1/1 1/1 1/1

  11. STATUS=2 0/1 0/1 0/1

  12. #那么总结如下：

  13. 26564 STATUS=1 无所以无 （0/0 0/1 0/0或者 0/1 0/0 0/0等等）

  14. 97764 STATUS=2 有所以有 （1/1 1/1 1/1 或者0/1 0/1 1/1等等）

  15. 385 STATUS=3 无中生有 （0/0 0/0 0/1 或者0/0 0/0 1/1）

  16. 1485 STATUS=4 有中生无 （1/1 0/1 0/0 等等）

我用annovar注释了一下

  1. /home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4  trio.mpileup.output.snp.vcf > trio.snp.annovar

  2. /home/jmzeng/bio-soft/annovar/annotate_variation.pl -buildver hg19  --geneanno --outfile  trio.snp.anno trio.snp.annovar /home/jmzeng/bio-soft/annovar/humandb

结果是：

  1. A total of 132268 locus in VCF file passed QC threshold, representing 132809 SNPs (86633 transitions and 46176 transversions) and 3 indels/substitutions

可以看到最后被注释到外显子上面的突变有两万多个。

  1. 23794  284345 3123333 trio.snp.anno.exonic_variant_function

这个应该是非常有意义的，但是我还没学会后面的分析。只能先做到这里了……

后记：

被软件找到的de novo 变异位点有502个，而文献参考认为WGS角度来看，正常人的繁殖后代新发突变应该是37个左右，如果只测了WES，那么真正的de novo 变异位点应该再少一个数量级，但是NGS测序就是这样，技术的限制，找到的位点非常多，需要用IGV载入bam和vcf去仔细人工检查看看。

因为当时理解有限，现在这些数据又丢失了，所以这个分析就虎头蛇尾了。

赞 (0)

相关推荐

  • NGS数据分析实践：03. 涉及的常用数据格式[5] - vcf格式

    系列文章: 二代测序方法:DNA测序之靶向重测序 NGS数据分析实践:00. 变异识别的基本流程 NGS数据分析实践:01. Conda环境配置及软件安装 NGS数据分析实践:02. 参考基因组及注释 ...

  • 【软件介绍】ANNOVAR注释软件用法

    变异检测得到的结果是检测样本的基因组序列与参考基因组序列之间的差异.本质上是一个将真实的变异从文库准备.样本富集.检测/测序和映射/比对产生的产物中分离出来的过程.想要进一步研究每一个变异的实际意义, ...

  • 【录音回放】专利分析五步走：从专利情报到制胜决策

    本文根据科睿唯安解决方案顾问武彬主讲的<专利情报和专利分析--基于重载轴承技术专利分析案例>网络研讨会内容整理后编制而成. 专利情报是法律信息.技术信息.商业信息的复合体,是所有创新型企业 ...

  • 基层医声：高血压管理七步走

    基层医声：高血压管理七步走

  • 曹锟纪念币鉴定七步走

    曹锟纪念币鉴定七步走

  • 草书入门实战七步走

    草书之难,在于草书将篆隶楷行诸体之笔意融会贯通浑然一体,却又不可得见端倪.对此,刘熙载说得好:"草书之笔划,要无一可移他书,而他书之笔意,草书却要无所不悟". 草书之难,又在于点划 ...

  • 学会草书入门七步走，就能克服新手时期的畏难心理！

    学好草书 3篇原创内容 公众号 欢迎您 查看草书名帖 草书之难,在于草书将篆隶楷行诸体之笔意融会贯通浑然一体,却又不可得见端倪.对此,刘熙载说得好:"草书之笔划,要无一可移他书,而他书之笔意 ...

  • 第十九战：背水一战（4）兵仙的七步走

    韩进率汉军进入了井陉. 在大军行进到距离井陉道东口还有三十里的时候,韩信下令停止进军,让战士们早早休息. 半夜,韩信下令打埋伏专业户灌婴选拔两千轻骑,起营整军,偃旗息鼓,秘密开拔迅速推进,每人发红旗( ...

  • 陈春花：品牌发展之路七步走

    导读:一个成功的品牌一定会获得高的价值认同,致力于构建品牌的企业,也往往会获得持续的发展.品牌发展之路归纳为以下七步:识别力量:价值链管理:始终如一交付价值的经理:清晰沟通价值的员工:可细分的忠诚的顾 ...

  • 七步走，齐步走。

    <永远不要停下前进的脚步> 今天,是你陪着老周的第744天! 有许多问题,出现在脑子里时,感觉特别简单,手到擒来.但出现在生活的具体情境中,一下子会感觉眼高手低,手足无措,说不清道不明. ...

  • 习拳按这七步走，想不出功夫都难！

    太极拳,看起来简单,练起来却奥妙无穷,想要练出太极功除了注意练呼吸.吐纳.意念外,每个阶段都有不同的要术,习练者必须要明了: 基础阶段的技术训练 要领:体松心静.身体中正.安舒自然.动作,技法规范准确 ...