全基因组数据分析完整流程+外显子测序数据分析例子

但如果没条件的话，可以尝试利用现有的数据（比如：千人基因组项目，GIAB  http://jimb.stanford.edu/giab/等）复现它们的成果，甚至只是构建一个分析流程也行，这样子学起来才会比较高效，同时也有利于夯实所学的知识。

找例子：数据库：ATGG

流程图：

流程图


常用软件：bwa，samtools，picard，GATK，bedtools，bcftools，vcftools，FastQC，MultiQC，VEP

例子：http://rosalind.info/problems/locations/    包括：RNA，DNA，蛋白

GWAS

全基因组测序分布-China

软件：比对软件-BWA；

     数据转换-samtools

            - Picard

基因数据变异检测工具-GATK

[if !supportLists]1.   [endif]质控：

     软件FASTQC质控，


1.   [endif]read各个位置的碱基质量值分布

2.   [endif]碱基的总体质量值分布

3.   [endif]read各个位置上碱基分布比例，目的是为了分析碱基的分离程度

4.   [endif]GC含量分布

5.   [endif]read各位置的N含量

6.   [endif]read是否还包含测序的接头序列

7.   [endif]read重复率，这个是实验的扩增过程所引入的

2． 比对

BWA比对


 通常把同一个样本的所有比对结果变成一个：因为有的测序深度深或者需要不同的lane来测序； $ samtools merge [ ...]

gatk \

-T RealignerTargetCreator \    ##目的是定位出所有需要进行序列重比对的目标区域

 -R /path/to/human.fasta\    

 -Isample_name.sorted.markdup.bam \

 -known/path/to/gatk/bundle/1000G_phase1.indels.b37.vcf \

 -known/path/to/gatk/bundle/Mills_and_1000G_gold_standard.indels.b37.vcf \

 -osample_name.IndelRealigner.intervals


Gatk \

-T IndelRealigner \     ##  对所有在第一步中找到的目标区域运用算法进行序列重比对，最后得到捋顺了的新结果。

-R /path/to/human.fasta \

-I sample_name.sorted.markdup.bam \

-known /path/to/gatk/bundle/1000G_phase1.indels.b37.vcf\

-known /path/to/gatk/bundle/Mills_and_1000G_gold_standard.indels.b37.vcf \

-o sample_name.sorted.markdup.realign.bam \

--targetIntervals sample_name.IndelRealigner.intervals


注意：1000G_phase1.indels.b37.vcf和Mills_and_1000G_gold_standard.indels.b37.vcf。这两个文件来自于千人基因组和Mills项目，里面记录了那些项目中检测到的人群Indel区域。

重新校正碱基质量值（BQSR）Base Quality Score Recalibration

主要是通过机器学习的方法构建测序碱基的错误率模型，然后对这些碱基的质量值进行相应的调整。

流程：

gatk \

-T BaseRecalibrator \

-R /path/to/human.fasta \

-I sample_name.sorted.markdup.realign.bam \

--knownSites /path/to/gatk/bundle/1000G_phase1.indels.b37.vcf \

--knownSites /path/to/gatk/bundle/Mills_and_1000G_gold_standard.indels.b37.vcf \

--knownSites /path/to/gatk/bundle/dbsnp_138.b37.vcf \

-o sample_name.recal_data.table

###BaseRecalibrator，这里计算出了所有需要进行重校正的read和特征值，然后把这些信息输出为一份校准表文件（sample_name.recal_data.table）

gatk

-T PrintReads \

-R /path/to/human.fasta \

-I sample_name.sorted.markdup.realign.bam \

--BQSR sample_name.recal_data.table \

-o sample_name.sorted.markdup.realign.BQSR.bam

##PrintReads，这一步利用第一步得到的校准表文件（sample_name.recal_data.table）重新调整原来BAM文件中的碱基质量值，并使用这个新的质量值重新输出一份新的BAM文件。


变异检测   SNP、Indel，CNV和SV等

  方法：使用GATK

HaplotypeCaller模块对样本中的变异进行检测，它也是目前最适合用于对二倍体基因组进行变异（SNP+Indel）检测的算法。

一般来说，在实际的WGS流程中对HaplotypeCaller的应用有两种做法，差别只在于要不要在中间生成一个gVCF：

   （1）直接进行HaplotypeCaller，这适合于单样本，或者那种固定样本数量的情况，也就是执行一次HaplotypeCaller之后就老死不相往来了。否则你会碰到仅仅只是增加一个样本就得重新运行这个HaplotypeCaller的坑爹情况（即，N+1难题），而这个时候算法需要重新去读取所有人的BAM文件，这将会是一个很费时间的痛苦过程；

gatk \

-T HaplotypeCaller \

-R /path/to/human.fasta \

-I sample_name.sorted.markdup.realign.BQSR.bam \

-D /path/to/gatk/bundle/dbsnp_138.b37.vcf \

-stand_call_conf 50 \

 -A QualByDepth \

 -A RMSMappingQuality \

 -A MappingQualityRankSumTest \

 -A ReadPosRankSumTest \

 -A FisherStrand \

 -A StrandOddsRatio \

 -A Coverage \

-o sample_name.HC.vcf

这里我特别提一下-D参数输入的dbSNP同样可以再GATK bundle目录中找到，这份文件汇集的是目前几乎所有的公开人群变异数据集。另外，由于我们的例子只有一个样本因此只输入一个BAM文件就可以了，如果有多个样本那么可以继续用-I参数输入：


（2）每个样本先各自生成gVCF，然后再进行群体joint-genotype。这其实就是GATK团队为了解决（1）中的N+1难题而设计出来的模式。gVCF全称是genome VCF，是每个样本用于变异检测的中间文件，格式类似于VCF，它把joint-genotype过程中所需的所有信息都记录在这里面，文件无论是大小还是数据量都远远小于原来的BAM文件。这样一旦新增加样本也不需要再重新去读取所有人的BAM文件了，只需为新样本生成一份gVCF，然后重新执行这个joint-genotype就行了。

  gatk \

 -T HaplotypeCaller \

    -R reference.fasta \

    -I sample1.bam [-I sample2.bam ...] \

    ...一样了

然而，基因组上各个不同的染色体之间其实是可以理解为相互独立的（结构性变异除外），也就是说，为了提高效率我们可以按照染色体一条条来独立执行这个步骤，最后再把结果合并起来就好了，这样的话就能够节省很多的时间。下面我给出一个按照染色体区分的例子：

gatk \

 -T HaplotypeCaller \

 -R /path/to/human.fasta \

 -I sample_name.sorted.markdup.realign.BQSR.bam\

 -D /path/to/gatk/bundle/dbsnp_138.b37.vcf \

 -L 1 \

 -stand_call_conf 50 \

 -A QualByDepth \

 -A RMSMappingQuality \

 -A MappingQualityRankSumTest \

 -A ReadPosRankSumTest \

 -A FisherStrand \

 -A StrandOddsRatio \

 -A Coverage \

 -o sample_name.HC.1.vcf


注意到了吗？其它参数都没任何改变，就只增加了一个 -L 参数，通过这个参数我们可以指定特定的染色体（或者基因组区域）！我们这里指定的是 1 号染色体，有些地方会写成chr1，具体看human.fasta中如何命名，与其保持一致即可。其他染色体的做法也是如此，就不再举例了。最后合并：

   gatk \

 -T CombineVariants \

 -R /path/to/human.fasta \

 --genotypemergeoption UNSORTED \

 --variant sample_name.HC.1.vcf \

 --variant sample_name.HC.2.vcf \

 ...

 --variant sample_name.HC.MT.vcf \

 -o sample_name.HC.vcf


第二种，先产生gVCF，最后再joint-genotype的做法：


java -jar /path/to/GenomeAnalysisTK.jar \

 -THaplotypeCaller \

 -R/path/to/human.fasta \

 -Isample_name.sorted.markdup.realign.BQSR.bam \

 --emitRefConfidence GVCF \

 -osample_name.g.vcf


#调用GenotypeGVCFs完成变异calling

java -jar /path/to/GenomeAnalysisTK.jar \

 -TGenotypeGVCFs \

 -R/path/to/human.fasta \

 --variant sample_name.g.vcf \

 -osample_name.HC.vcf


其实，就是加了–emitRefConfidence GVCF的参数。而且，假如嫌慢，同样可以按照染色体或者区域去产生一个样本的gVCF，然后在GenotypeGVCFs中把它们全部作为输入文件完成变异calling。也许你会担心同个样本被分成多份gVCF之后，是否会被当作不同的多个样本？回答是不会！因为生成gVCF文件的过程中，GATK会根据@RG信息中的SM（也就是sample name）来判断这些gVCF是否来自同一个样本，如果名字相同，那么就会被认为是同一个样本，不会产生多样本问题。

变异检测质控和过滤（VQSR）

这是我们这个流程中最后的一步了。在获得了原始的变异检测结果之后，我们还需要做的就是质控和过滤。这一步或多或少都有着一些个性化的要求，我暂时就不做太多解释吧（一旦解释恐怕同样是一篇万字长文）。只用一句话来概括，VQSR是通过构建GMM模型对好和坏的变异进行区分，从而实现对变异的质控，具体的原理暂时不展开了。


下面就直接给出例子吧：


## SNP Recalibrator

java -jar /path/to/GenomeAnalysisTK.jar \

   -TVariantRecalibrator \

   -Rreference.fasta \

  -input sample_name.HC.vcf \

  -resource:hapmap,known=false,training=true,truth=true,prior=15.0/path/to/gatk/bundle/hapmap_3.3.b37.vcf \

  -resource:omini,known=false,training=true,truth=false,prior=12.0/path/to/gatk/bundle/1000G_omni2.5.b37.vcf \

  -resource:1000G,known=false,training=true,truth=false,prior=10.0/path/to/gatk/bundle/1000G_phase1.snps.high_confidence.b37.vcf \

  -resource:dbsnp,known=true,training=false,truth=false,prior=6.0/path/to/gatk/bundle/dbsnp_138.b37.vcf \

  -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP \

  -mode SNP \

  -recalFile sample_name.HC.snps.recal \

  -tranchesFile sample_name.HC.snps.tranches \

  -rscriptFile sample_name.HC.snps.plots.R

java -jar /path/to/GenomeAnalysisTK.jar -TApplyRecalibration \

  -R  human_g1k_v37.fasta \

  -input sample_name.HC.vcf \

  --ts_filter_level 99.5 \

  -tranchesFile sample_name.HC.snps.tranches \

  -recalFile sample_name.HC.snps.recal \

  -mode SNP \

   -osample_name.HC.snps.VQSR.vcf

## Indel Recalibrator

java -jar /path/to/GenomeAnalysisTK.jar -TVariantRecalibrator \

  -R  human_g1k_v37.fasta \

  -input sample_name.HC.snps.VQSR.vcf \

  -resource:mills,known=true,training=true,truth=true,prior=12.0/path/to/gatk/bundle/Mills_and_1000G_gold_standard.indels.b37.vcf \

  -an QD -an DP -an FS -an SOR -an ReadPosRankSum -an MQRankSum \

  -mode INDEL \

  -recalFile sample_name.HC.snps.indels.recal \

  -tranchesFile sample_name.HC.snps.indels.tranches \

  -rscriptFile sample_name.HC.snps.indels.plots.R

java -jar /path/to/GenomeAnalysisTK.jar -TApplyRecalibration \

   -Rhuman_g1k_v37.fasta\

  -input sample_name.HC.snps.VQSR.vcf \

  --ts_filter_level 99.0 \

  -tranchesFile sample_name.HC.snps.indels.tranches \

  -recalFile sample_name.HC.snps.indels.recal \

  -mode INDEL \

   -osample_name.HC.snps.indels.VQSR.vcf

最后，sample_name.HC.snps.indels.VQSR.vcf便是我们最终的变异检测结果。对于人类而言，一般来说，每个人最后检测到的变异数据大概在400万左右（包括SNP和Indel）。

http://www.huangshujia.me/2017/09/19/2017-09-19-Begining-WGS-Data-Analysis-The-pipeline.html

例子：http://starsyi.github.io/2016/05/24/BWA-%E5%91%BD%E4%BB%A4%E8%AF%A6%E8%A7%A3/

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 211,817评论 6赞 492
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,329评论 3赞 385
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 157,354评论 0赞 348
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,498评论 1赞 284
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 65,600评论 6赞 386
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,829评论 1赞 290
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,979评论 3赞 408
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,722评论 0赞 266
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,189评论 1赞 303
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,519评论 2赞 327
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,654评论 1赞 340
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,329评论 4赞 330
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,940评论 3赞 313
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,762评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,993评论 1赞 266
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,382评论 2赞 360
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,543评论 2赞 349

推荐阅读更多精彩内容