【直播】我的基因组24：用GATK对SAM格式的文件进行重排

GATK教程我以前写过（:）

这个步骤分成2个小步骤：

首先用RealignerTargetCreator找到需要重新比对的区域，输出文件intervals，然后用输出的 tmp.intervals 做输入文件来进行重新比对，也就是用IndelRealigner在这些区域内进行重新比对

 这个是批量运行的代码，就是用shell脚本写一个循环即可（:）。

sam/bam的染色体顺序不一样，还有可能是染色体的名字不一样，比如>1和>chr1的区别。

下面是对我的基因组bam文件用GATK进行重排，命令分成两步，请注意下面用到的参考基因组文件和软件，需要保证都是按照前面的教程安装，用的的vcf文件下载地址见文末

nohup java -Xmx60g -jar ~/biosoft/GATK/GenomeAnalysisTK.jar -R ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta -T RealignerTargetCreator -I P_jmzeng.filter.rmdup.bam  -o P_jmzeng.filter.rmdup.intervals -known ~/annotation/GATK/1000G_phase1.indels.b37.vcf.gz 1> P_jmzeng.filter.rmdup.RealignerTargetCreator.log

## 请保证你的bam文件比对的时候选择的参考基因组跟你在使用GATK的时候自己指定的参考基因组的一致性。

nohup java -Xmx60g -jar ~/biosoft/GATK/GenomeAnalysisTK.jar  -R ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta  -T IndelRealigner  -I P_jmzeng.filter.rmdup.bam  -targetIntervals P_jmzeng.filter.rmdup.intervals -o P_jmzeng.filter.rmdup.realgn.bam  1>P_jmzeng.filter.rmdup.IndelRealigner.log

输出文件如下：

（nohup的命令在之前讲过，是后台运行的命令，其中1表示运行的日志输入log文件，如果用2的话，表示的是指输出报错的日志）

在探索这个软件的用法的时候可能会报错，根据报错日志搜索解决即可;

第一个错误：

MESSAGE: Fasta dict file /home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.dict for reference /home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta does not exist. Please see http://gatkforums.broadinstitute.org/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference for help creating it.

解决方案，在参考基因组文件所在的目录运行下面代码：

java -jar ~/biosoft/picardtools/picard-tools-1.119/CreateSequenceDictionary.jar R=human_g1k_v37.fasta O=human_g1k_v37.dict

第二个错误：

MESSAGE: An index is required,  K/1000G_phase1.indels.b37.vcf.gz

因为我在broad 网站下载的文件需要gunzip解压。 

文件下载地址是：  

下载代码是：

## ftp://ftp.broadinstitute.org/bundle/b37/

mkdir -p ~/annotation/GATK

cd ~/annotation/variation/GATK

wget ftp://[e-mail]/bundle/b37/1000G_phase1.snps.high_confidence.b37.vcf.gz

wget ftp://[e-mail]/bundle/b37/dbsnp_138.b37.vcf.gz

wget ftp://[e-mail]/bundle/b37/human_g1k_v37.fasta.gz

wget ftp://[e-mail]/bundle/b37/NA12878.HiSeq.WGS.bwa.cleaned.raw.subset.b37.sites.vcf.gz

wget ftp://[e-mail]/bundle/b37/Mills_and_1000G_gold_standard.indels.b37.vcf.gz

wget ftp://[e-mail]/bundle/b37/hapmap_3.3.b37.vcf.gz

wget ftp://[e-mail]/bundle/b37/1000G_phase1.indels.b37.vcf.gz

wget ftp://[e-mail]/bundle/b37/1000G_phase1.indels.b37.vcf.idx.gz

gunzip 1000G_phase1.indels.b37.vcf.idx.gz

gunzip 1000G_phase1.indels.b37.vcf.gz

ref:

文：Jimmy、阿尔的太阳

图文编辑：吃瓜群众