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