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，即使加上侧翼长度。

54692160：外显子加上前后50bp
73066288：外显子加上前后100bp
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转换）

awk '{print "chr"$3,$4,$5,$1,0,$2,$4,$5,"255,0,0"}' CCDS.20150512.exon.txt >CCDS.20150512.exon.hg38.bed
~/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数量如下

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
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
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
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
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
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
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
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
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
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
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
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

/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4 Sample3.varscan.snp.vcf > Sample3.annovar
/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4 Sample4.varscan.snp.vcf > Sample4.annovar
/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

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

23515 Sample3.anno.exonic_variant_function
23913 Sample4.anno.exonic_variant_function
24009 Sample5.anno.exonic_variant_function

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

downstream
exonic
exonic;splicing
intergenic
intronic
ncRNA_exonic
ncRNA_intronic
ncRNA_splicing
ncRNA_UTR3
ncRNA_UTR5
splicing
upstream
upstream;downstream
UTR3
UTR5
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

Min coverage: 10
Min reads2: 4
Min var freq: 0.2
Min avg qual: 15
P-value thresh: 0.05
Adj. min reads2: 2
Adj. var freq: 0.05
Adj. p-value: 0.15

Reading input from trio.filter.mpileup

1371416525 bases in pileup file （137M的序列）
83123183 met the coverage requirement of 10 （其中有83M的测序深度大于10X）
145104 variant positions (132268 SNP, 12836 indel) （共发现15.5万的变异位点）
4403 were failed by the strand-filter
139153 variant positions reported (126762 SNP, 12391 indel)
502 de novo mutations reported (376 SNP, 126 indel) （真正属于 de novo mutations只有502个）
1734 initially DeNovo were re-called Germline
12 initially DeNovo were re-called MIE
3 initially DeNovo were re-called MultAlleles
522 initially MIE were re-called Germline
1 initially MIE were re-called MultAlleles
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文件进行简化，只看基因型！

head status.txt （顺序是dad,mom,child）
STATUS=2 0/0 0/1 0/1
STATUS=2 1/1 1/1 1/1
STATUS=2 0/1 0/0 0/1
STATUS=2 1/1 1/1 1/1
STATUS=1 0/1 0/0 0/0
STATUS=1 0/1 0/0 0/0
STATUS=2 1/1 1/1 1/1
STATUS=2 1/1 1/1 1/1
STATUS=2 1/1 1/1 1/1
STATUS=2 0/1 0/1 0/1
#那么总结如下：
26564 STATUS=1 无所以无（0/0 0/1 0/0或者 0/1 0/0 0/0等等）
97764 STATUS=2 有所以有（1/1 1/1 1/1 或者0/1 0/1 1/1等等）
385 STATUS=3 无中生有（0/0 0/0 0/1 或者0/0 0/0 1/1）
1485 STATUS=4 有中生无（1/1 0/1 0/0 等等）

我用annovar注释了一下

/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4 trio.mpileup.output.snp.vcf > trio.snp.annovar
/home/jmzeng/bio-soft/annovar/annotate_variation.pl -buildver hg19 --geneanno --outfile trio.snp.anno trio.snp.annovar /home/jmzeng/bio-soft/annovar/humandb

结果是：

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

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

23794 284345 3123333 trio.snp.anno.exonic_variant_function

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

后记：

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

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

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。