生物信息学

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

2021-01-26  本文已影响0人  BeeBee生信

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 文件直接转换就差不多了。

上一篇下一篇

猜你喜欢

热点阅读