从 fasta 创建用于 GATK 流程的 interval list

GATKbundle 提供了用于 WGS 的 interval list 文件,但如果你用的参考基因组跟他提供的 interval list 不一致,比方说 Contig 命名不同,就容易造成问题。像我喜欢用 GENCODE 的参考基因组,那么需要自己制作 interval list 文件。
下面的代码是需要的步骤和注释。

# 用 fasta 创建 faidx 和 dict 索引
gatk CreateSequenceDictionary -R GRCh38.primary_assembly.genome.fa
samtools faidx GRCh38.primary_assembly.genome.fa

# 用 fadix 索引文件生成参考基因组的 bed 文件
awk -v FS="\t" -v OFS="\t" '{print $1 FS "0" FS ($2-1)}' GRCh38.primary_assembly.genome.fa.fai > GRCh38.primary_assembly.genome.bed
# 移除 blacklist,如果有其他区域需要移除,同样操作
bedtools subtract -a GRCh38.primary_assembly.genome.bed -b ../hg38.blacklist.bed > GRCh38.primary_assembly.genome.interval.bed
# 从 bed 格式转换到 interval list 格式
gatk BedToIntervalList -I GRCh38.primary_assembly.genome.interval.bed -O \ 
GRCh38.primary_assembly.genome.interval_list -SD GRCh38.primary_assembly.genome.dict

要注意 picard style 的 interval list 命名模式为 *.interval_list. GATK style 的为 *.list.

如果是 WES 那么从测序厂商拿到合适的 bed 文件直接转换就差不多了。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容

  • 00 写在前面 仅针对人类WGS或WES数据,供参考。时间管理某一点:能自动化的工作尽量自动化,不要时间用在毫无意...
    fatlady阅读 13,029评论 23 54
  • 上次我们整理到bwa比对后得到bam文件,下一步我们要通过GATK流程从bam文件中call variant。 一...
    耕读者阅读 2,059评论 0 4
  • GATK(全称The Genome Analysis Toolkit)是Broad Institute开发的用于二...
    Greatji阅读 34,240评论 0 18
  • 上次我们整理到bwa比对后得到bam文件,下一步我们要通过GATK流程从bam文件中call variant。 一...
    耕读者阅读 2,296评论 0 2
  •   GATK现在最新的稳定版已经到了3.8,测试版是4.0。3.8版和之前的版本还是有比较大的不同的,但核心算法与...
    生信杂谈阅读 27,459评论 1 51