搜索
查看: 1903|回复: 5

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

[复制链接]

103

主题

133

帖子

814

积分

版主

Rank: 7Rank: 7Rank: 7

积分
814
发表于 2017-2-2 16:42:35 | 显示全部楼层 |阅读模式
【直播】我的基因组24:用GATK对SAM格式的文件进行重排

GATK教程我以前写过([url=]GATK使用注意事项[/url]:[url=]http://www.bio-info-trainee.com/838.html[/url]
这个步骤分成2个小步骤:

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



这个是批量运行的代码,就是用shell脚本写一个循环即可([url=]WES(二)snp-calling[/url]:[url=]http://www.bio-info-trainee.com/1114.html[/url])。

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

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

[AppleScript] 纯文本查看 复制代码
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的时候自己指定的参考基因组的一致性。


[AppleScript] 纯文本查看 复制代码
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的话,表示的是指输出报错的日志)
在探索这个软件的用法的时候可能会报错,根据报错日志搜索解决即可;
第一个错误:
[AppleScript] 纯文本查看 复制代码
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 [url]http://gatkforums.broadinstitute.org/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference[/url] for help creating it.

解决方案,在参考基因组文件所在的目录运行下面代码:
[AppleScript] 纯文本查看 复制代码
java -jar ~/biosoft/picardtools/picard-tools-1.119/CreateSequenceDictionary.jar R=human_g1k_v37.fasta O=human_g1k_v37.dict

第二个错误:
[AppleScript] 纯文本查看 复制代码
MESSAGE: An index is required,  K/1000G_phase1.indels.b37.vcf.gz

因为我在broad 网站下载的文件需要gunzip解压。
文件下载地址是: [url=]ftp://ftp.broadinstitute.org/bundle/[/url]  
下载代码是:
[AppleScript] 纯文本查看 复制代码
## [url]ftp://ftp.broadinstitute.org/bundle/b37/[/url]
mkdir -p ~/annotation/GATK

cd ~/annotation/variation/GATK

wget [url]ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/1000G_phase1.snps.high_confidence.b37.vcf.gz[/url]

wget [url]ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/dbsnp_138.b37.vcf.gz[/url]

wget [url]ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/human_g1k_v37.fasta.gz[/url]

wget [url]ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/NA12878.HiSeq.WGS.bwa.cleaned.raw.subset.b37.sites.vcf.gz[/url]

wget [url]ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/Mills_and_1000G_gold_standard.indels.b37.vcf.gz[/url]

wget [url]ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/hapmap_3.3.b37.vcf.gz[/url]

wget [url]ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/1000G_phase1.indels.b37.vcf.gz[/url]

wget [url]ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/1000G_phase1.indels.b37.vcf.idx.gz[/url]

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

gunzip 1000G_phase1.indels.b37.vcf.gz

ref:
[url=]https://software.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_indels_IndelRealigner.php[/url]
[url=]http://gatkforums.broadinstitute.org/gatk/discussion/1213/what-s-in-the-resource-bundle-and-how-can-i-get-it[/url]
[url=]http://sfg.stanford.edu/SNP.html[/url]
[url=]https://redmine.igm.cumc.columbia.edu/projects/biopipeline/wiki/GATK[/url]

文:Jimmy、阿尔的太阳

图文编辑:吃瓜群众



本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:【直播】我的基因组23:对比对结果文件进行过滤
下一篇:【直播】我的基因组25:用bcftools来call variation
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

0

主题

11

帖子

199

积分

注册会员

Rank: 2

积分
199
发表于 2017-11-8 15:43:12 | 显示全部楼层
疯了,大神,不懂我的基因组做到一般出现了啥问题?~~~
一直提示这个错误:ERROR MESSAGE: Lexicographically sorted human genome sequence detected in reference. Please see https://software.broadinstitute. ... tion/article?id=132
8for more information. Error details: reference contigs = [chr10, chr11, chr11_gl000202_random, chr12, chr13, chr14, chr15, chr16, chr17_ctg5_hap1, chr17, chr17_gl000203_r
andom, chr17_gl000204_random, chr17_gl000205_random, chr17_gl000206_random, chr18, chr18_gl000207_random, chr19, chr19_gl000208_random, chr19_gl000209_random, chr1, chr1_g
l000191_random, chr1_gl000192_random, chr20, chr21, chr21_gl000210_random, chr22, chr2, chr3, chr4_ctg9_hap1, chr4, chr4_gl000193_random, chr4_gl000194_random, chr5, chr6_
apd_hap1, chr6_cox_hap2, chr6_dbb_hap3, chr6, chr6_mann_hap4, chr6_mcf_hap5, chr6_qbl_hap6, chr6_ssto_hap7, chr7, chr7_gl000195_random, chr8, chr8_gl000196_random, chr8_gl
000197_random, chr9, chr9_gl000198_random, chr9_gl000199_random, chr9_gl000200_random, chr9_gl000201_random, chrM, chrUn_gl000211, chrUn_gl000212, chrUn_gl000213, chrUn_gl
000214, chrUn_gl000215, chrUn_gl000216, chrUn_gl000217, chrUn_gl000218, chrUn_gl000219, chrUn_gl000220, chrUn_gl000221, chrUn_gl000222, chrUn_gl000223, chrUn_gl000224, chr
Un_gl000225, chrUn_gl000226, chrUn_gl000227, chrUn_gl000228, chrUn_gl000229, chrUn_gl000230, chrUn_gl000231, chrUn_gl000232, chrUn_gl000233, chrUn_gl000234, chrUn_gl000235
, chrUn_gl000236, chrUn_gl000237, chrUn_gl000238, chrUn_gl000239, chrUn_gl000240, chrUn_gl000241, chrUn_gl000242, chrUn_gl000243, chrUn_gl000244, chrUn_gl000245, chrUn_gl0
00246, chrUn_gl000247, chrUn_gl000248, chrUn_gl000249, chrX, chrY]
;然后我根据网页上的提示,跑了个ReorderSam,结果是这个Exception in thread "main" htsjdk.samtools.SAMFormatException: SAM validation error: ERROR: Record 1490928, Read name B80KJTABXX:1:3:1174:29885#NNNNNNNN, bin field of BAM
record does not equal value computed based on alignment start and end, and length of sequence to which read is aligned
这是啥情况呢~~~~
我用的是hg19的参考基因组~~~~
我前面做过水稻的snp也用的GATK都没这些问题~卡了我好久哇~~~~~求解答!!!!
回复 支持 反对

使用道具 举报

0

主题

11

帖子

199

积分

注册会员

Rank: 2

积分
199
发表于 2017-11-8 15:52:08 | 显示全部楼层
感谢感谢~~~~~
回复

使用道具 举报

1

主题

4

帖子

84

积分

注册会员

Rank: 2

积分
84
发表于 2018-1-11 14:33:54 | 显示全部楼层
393677116@qq.co 发表于 2017-11-8 15:43
疯了,大神,不懂我的基因组做到一般出现了啥问题?~~~
一直提示这个错误:ERROR MESSAGE: Lexicographical ...

你这是sam文件没有按染色体排序?
回复 支持 反对

使用道具 举报

0

主题

11

帖子

199

积分

注册会员

Rank: 2

积分
199
发表于 2018-2-2 10:54:21 | 显示全部楼层
lz13215962853 发表于 2018-1-11 14:33
你这是sam文件没有按染色体排序?

我sort了哦~~
回复 支持 反对

使用道具 举报

1

主题

46

帖子

584

积分

高级会员

Rank: 4

积分
584
发表于 2018-2-2 15:26:31 | 显示全部楼层
请问,为什么要对sam文件进行排序?不排序结果影响很大吗
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2018-5-21 05:42 , Processed in 0.101886 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.