搜索
查看: 4849|回复: 2

【直播】我的基因组(七):从整体理解全基因组测序数据...

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-1-28 10:28:52 | 显示全部楼层 |阅读模式
本帖最后由 zckoo007 于 2017-1-28 10:31 编辑

【直播】我的基因组(七):从整体理解全基因组测序数据的变异位点
目次
1.    变异与突变的区别
2.    变异分析流程
3.    变异分析代码

1    变异与突变的区别
首先记住一个很重要的知识点,变异是相对的
简单说一下什么是找变异,变异跟突变有什么区别呢?举个栗子:有国际组织规定了人类的参考基因组(如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,但总之,我们得先找到自己的测序数据跟国际规定的参考基因组有什么区别(变异)吧!


2    变异分析流程
变异分成4种,即snv、indel、cnv、sv,大部分情况下只能分析到SNV,另外3个要么不准确,要么有点难度!
bwa软件的作者,大名鼎鼎的 Heng Li给出的流程如下: http://www.htslib.org/workflow/
根据Heng Li的博客自己也完成过几十个外显子数据的找变异分析,其中还包括一个自闭症家系的分析,通过与参考基因组比较找到变异并不难,但是如何给找到的几万到几百万个变异一个合理的解释才是问题所在。
第一步,下载数据
第二步,bwa比对
第三步,sam转为bam,并sort好
第四步,标记PCR重复,并去除
第五步,产生需要重排的坐标记录
第六步,根据重排记录文件把比对结果重新比对
第七步,把最终的bam文件转为mpileup文件
第八步,用bcftools来call snp
第九步,用freebayes来call snp
第十步,用gatk来call snp
第十一步,用varscan来call snp


本次处理全基因组数据我也准备走同样的流程,因为找到变异并不是重点,即使中间有什么不妥,我们也可以随时回过头来看看问题出在哪里!
其中需要安装的软件及参考基因组及注释文件在我之前的文章里都提到了。



3    变异分析代码
大家可以简单用下面的代码处理一下KPGP0001这个个体的全基因组测序数据,如下:

[Perl] 纯文本查看 复制代码
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



d26fcc006d1fca3b4cfa4fc21fa73ec8.jpg
有学者处理了Korean Personal Genomes Project (KPGP)中的35 Koreangenomes里面的WGS数据,文章中用了两套SNVcalling流程来处理:http://bmcbioinformatics.biomedc ... 1471-2105-15-S11-S6 流程如下,大家可以进行一下参考。
c10e0d0b59be8067ed96ba1f834c92dc.jpg




上一篇:生信编程直播第六题:批量根据基因list来提取信息(shell)
下一篇:【直播】我的基因组(八):原始测序数据质量报告
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

0

主题

7

帖子

512

积分

高级会员

Rank: 4

积分
512
发表于 2017-2-21 16:35:33 | 显示全部楼层
感谢感谢
回复

使用道具 举报

0

主题

1

帖子

31

积分

新手上路

Rank: 1

积分
31
发表于 2017-3-1 17:41:27 | 显示全部楼层
楼祝前辈您好,谢谢您的分享,我看您 bwa 使用的参数 是 mem,请问是否可以用 aln 呢? 这两种方法,我在进行外显子组测序分析和基因组重测序分析时,改如何选择呢?
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2020-6-1 13:38 , Processed in 0.026459 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.