GATK4 —— 获取短变异 (call SNP+indel)

GATK是一款用于基因组数据分析的软件,其强大的处理引擎和高性能计算功能使其能够承担任何规模的项目。

GATK的功能非常强大,这里不详细介绍,大家可以根据自己的要求,从首页进入对应的模块,说明书还是很详细的。

官网:
https://gatk.broadinstitute.org/hc/en-us

GATK 4.2.1.0 命令说明:
https://gatk.broadinstitute.org/hc/en-us/sections/4404575821339-4-2-1-0?page=10#articles

1. 下载安装

1.1下载

下载链接:
https://github.com/broadinstitute/gatk/releases

1.2 安装

解压:

$ unzip gatk-4.2.4.1.zip

查看帮助:

$ cd /your/path/gatk-4.2.4.1
$ ./gatk --help

查看工具列表:

$ ./gatk --list

查看指定工具帮助:

$ ./gatk ToolName --help

2.获取短变异(SNP+indel)

说明书:
https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels-
主要流程如下:

2.1 准备数据

$ samtools index LPF1_R1_MP.sort.dup.bam
  • 参考基因组(.fa)
  • 参考基因组索引(.fai)
$ samtools faidx Mparg_v2.0.fa
$ less Mparg_v2.0.fa.fai
Mpar_chr1       41449654        11      100     101
Mpar_chr2       50537509        41864173        100     101
Mpar_chr3       46570116        92907069        100     101
Mpar_chr4       46692298        139942898       100     101
Mpar_chr5       50691827        187102130       100     101
Mpar_chr6       65498417        238300887       100     101
Mpar_chr7       34026340        304454300       100     101
Mpar_chr8       48943667        338820915       100     101

五列分别对应,序列名、序列长度、第一个碱基的偏移量、行的碱基数、行宽。

  • 参考基因组索引(.dict)
$ gatk CreateSequenceDictionary -R Mparg_v2.0.fa -O Mparg_v2.0.dict

2.2 HaplotypeCaller

HaplotypeCaller是GATK4中的核心工具,可以利用单倍型区域的局部从头组装同时调用种系snv和小indel。详细原理参见说明书。
https://gatk.broadinstitute.org/hc/en-us/articles/360050814612-HaplotypeCaller#--sample-name

$ ./gatk HaplotypeCaller --help

Required Arguments:

--input,-I <GATKPath>         BAM/SAM/CRAM file containing reads  This argument must be specified at least once.
                              Required.

--output,-O <GATKPath>        File to which variants should be written  Required.

--reference,-R <GATKPath>     Reference sequence file  Required.

示例:

$ gatk --java-options "-Xmx100g -XX:ParallelGCThreads=4" HaplotypeCaller \
    -I ~/LPF1_R1_MP.sort.dup.bam \
    -R ~/Mparg_v2.0.fa \
    -ERC GVCF \
    -O ~/vcf/LPF1_R1_MP.g.vcf.gz

报错:A USER ERROR has occurred: Argument emit-ref-confidence has a bad value: Can only be used in single sample mode currently. Use the --sample-name argument to run on a single sample out of a multi-sample BAM file.

解决办法:
picard——修改BAM文件的Read Group - 简书 (jianshu.com)

更改命令:

$ samtools index LPF1_R1_MP.rg.sort.bam
$ gatk --java-options "-Xmx100g -XX:ParallelGCThreads=4" HaplotypeCaller \
    -I ~/LPF1_R1_MP.rg.sort.bam \
    -R ~/Mparg_v2.0.fa \
    -ERC GVCF \
    -O ~/LPF1_R1_MP.g.vcf.gz

gatk好像最多只支持4个线程,如果我我理解有错误的话,还请大佬指正。

这里注意HaplotypeCaller只能处理单样本文件,当有多样本时,官方建议使用HaplotypeCaller对单bam文件分别进行变异检测,生成GVCF文件,GVCF会记录每一个位点到情况,包括有无突变,VCF只记录突变位点情况,之后在下一步对GVCF文件进行合并。

VCF文件和GVCF文件的格式的详细说明参见下链接:
VCF和GVCF格式说明 - 青萍,你好 - 博客园 (cnblogs.com)

3.GVCF文件合并

官方说明书:https://gatk.broadinstitute.org/hc/en-us/articles/4404604801947-CombineGVCFs

gatk --java-options "-Xmx100g -XX:ParallelGCThreads=4" CombineGVCFs \
    -R ~/ref/Mparg_v2.0.fa \
    -V LPF1_R1_MP.g.vcf.gz \
    -V LPF1_R2_MP.g.vcf.gz \
    -V LPF1_R3_MP.g.vcf.gz \
    -O LPF1_MP.g.vcf.gz

参考资料:
https://www.melbournebioinformatics.org.au/tutorials/tutorials/variant_calling_gatk1/variant_calling_gatk1/

https://gencore.bio.nyu.edu/variant-calling-pipeline-gatk4/

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 199,830评论 5 468
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 83,992评论 2 376
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 146,875评论 0 331
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,837评论 1 271
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,734评论 5 360
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,091评论 1 277
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,550评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,217评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,368评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,298评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,350评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,027评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,623评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,706评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,940评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,349评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,936评论 2 341

推荐阅读更多精彩内容