【直播】我的基因组（七）:从整体理解全基因组测序数据的变异位点
简单说一下什么是找变异，变异跟突变有什么区别呢？举个栗子：有国际组织规定了人类的参考基因组（如UCSC,ENSEMBL,NCBI等，前面帖子都有讲)，就是 AAAAA(这里简化一下，就5个碱基，其实人类基因组多达30亿个) 。现在通过给自己测序得知，我与之对应的是AGCAA，那么我相比国际基因组来说，就是2个变异位点，位于基因组的坐标2和3，但是它们还不能说就是突变。
如第二位碱基，虽然我的是G，参考基因组是A，但是全球已经测序了几百万人，而我查看了他们的测序结果，其中99万人都是G，这说明是参考基因组出现了问题，可能是国际组织当年恰好选择了一个人是A，所以就规定第二个碱基是A。所以虽然我用软件找到了我的这个位点相对于参考基因组是来说，是一个变异，但是这恰好是好事，完全不用担心，我们也不需要用突变这个单词来描述它！
那么接下来看第3位碱基，同样，国际组织规定了是A，而我却测了个C，但是全球已经公布的一百万人里面99.999万人都跟参考一样，就是A。有一个人和参考基因组对应的碱基不一样，不一样的那个人是个有病的患者，这个时候，你就惨了，这个变异，就是突变了！
很多变异其实只是造成人种多样性的原因，是构成人独特性的基础，而那些跟疾病相关的变异，我们通常就会叫做是突变！因我我只举了2个极端的例子，所以大家可能会误以为，跟大多数人一样，就没事了！其实也并不是这样，一般来说，在正常人的数据库里面出现了5%的变异就可以认为没什么大的危害，而且变异还可以分成germline、somatic、de novo等情况，如果是特定性的针对某种疾病还可以找driver的mutation，但总之，我们得先找到自己的测序数据跟国际规定的参考基因组有什么区别（变异）吧！
变异分成4种，即snv、indel、cnv、sv，大部分情况下只能分析到SNV，另外3个要么不准确，要么有点难度！
bwa软件的作者，大名鼎鼎的 Heng Li给出的流程如下：
根据Heng Li的博客自己也完成过几十个外显子数据的找变异分析，其中还包括一个自闭症家系的分析，通过与参考基因组比较找到变异并不难，但是如何给找到的几万到几百万个变异一个合理的解释才是问题所在。
我当初的流程如下:()
第一步，下载数据
第二步，bwa比对
第三步，sam转为bam，并sort好
第四步，标记PCR重复，并去除
第五步，产生需要重排的坐标记录
第六步，根据重排记录文件把比对结果重新比对
第七步，把最终的bam文件转为mpileup文件
第八步，用bcftools 来call snp
第九步，用freebayes来call snp
第十步，用gatk来call snp
第十一步，用varscan来call snp本次处理全基因组数据我也准备走同样的流程，因为找到变异并不是重点，即使中间有什么不妥，我们也可以随时回过头来看看问题出在哪里！
其中需要安装的软件及参考基因组及注释文件在我之前的文章里都提到了。
大家可以简单用下面的代码处理一下KPGP0001这个个体的全基因组测序数据，如下：
ls *gz |xargs ~/biosoft/fastqc/FastQC/fastqc -t 10
for i in $(seq 1 6) ;do (nohup ~/biosoft/bwa/bwa-0.7.15/bwa mem -t 5 -M ~/reference/index/bwa/hg19 KPGP-00001_L${i}_R1.fq.gz KPGP-00001_L${i}_R2.fq.gz 1>KPGP-00001_L${i}.sam 2>KPGP-00001_L${i}.bwa.align.log &);done
for i in $(seq 1 6) ;do (nohup samtools sort -@ 5 -o KPGP-00001_L${i}.sorted.bam KPGP-00001_L${i}.sam &);done
for i in $(seq 1 6) ;do (nohup samtools index KPGP-00001_L${i}.sorted.bam &);done
samtools merge KPGP-00001.merge.bam *.sorted.bam
samtools sort -@ 50 -O bam -o KPGP-00001.sorted.merge.bam KPGP-00001.merge.bam
samtools index KPGP-00001.sorted.merge.bam
for i in $(seq 1 6) ;do ( samtools flagstat KPGP-00001_L${i}.sorted.bam >KPGP-00001_L${i}.flagstat.txt );done
有学者处理了Korean Personal Genomes Project (KPGP)中的 35 Korean genomes里面的WGS数据，文章中用了两套SNV calling流程来处理： 流程如下，大家可以进行一下参考。
文：Jimmy、吃瓜群众
图文编辑：吃瓜群众