How to create interval list from reference fasta or dict file?

I am using GATK pipeline on WGS data. My BAM files is aligned to GRCh38 from GENCODE. So I want to create interval file for this GRCh38 instead of download from GATKbundle, because some of their contigs have different names. For example "KI270706.1" in GENCODE's GRCh38 is "chr1_KI270706v1_random" in interval list downloaded from GATKbundle. I can't modify interval list downloaded from GATKbundle, because some contigs are different, for example exists in interval list but not GENCODE's GRCh38.

Can anyone tells me how to create my own interval list from GENCODE GRCh38? Thanks.

gatk

• 9,316 views • link

updated 2.7 years ago by Manuel Sokolov Ravasqueira •

• written 5.5 years ago by MatthewP

Hello @MatthewP , were you able to that for a fasta file for example? If you have a fasta file containing a targeted region, how to you generate a intervals.bed file for that fasta?

Thank you!

• 0 views • link

written 2.7 years ago by Manuel Sokolov Ravasqueira •

2 answers

Ok I can answer myself.

# Create dict and faidx index for fasta
gatk CreateSequenceDictionary -R GRCh38.primary_assembly.genome.fa
samtools faidx GRCh38.primary_assembly.genome.fa
# get fasta bed file
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
# remove blacklist
bedtools subtract -a GRCh38.primary_assembly.genome.bed -b ../hg38.blacklist.bed > GRCh38.primary_assembly.genome.interval.bed
# convert to interval list file
gatk BedToIntervalList -I GRCh38.primary_assembly.genome.interval.bed -O \ 
GRCh38.primary_assembly.genome.interval_list -SD GRCh38.primary_assembly.genome.dict

Note, this is picard style interval list, so named pattern is *.interval_list

• 0 views • link

written 5.5 years ago by MatthewP

Hi.

Thanks for your scripts! It is awesome.

Because call of intervals used bed files, I was wondering if gtf could do it as well. So I test gtf2bed command from bedops to see if gencode gtf could be used for this purpose.

It output a hundreds Mb interval list with many repeated information and ensembl_id as well. I did not test whether this test interval gonna work.

Searching for the whole internet, here is the most accurate answer for this question.

Thanks again.

• 0 views • link

written 5.0 years ago by tomas4482 •

Use picard tool

java -jar picard.jar CreateSequenceDictionary \ 
      R=reference.fasta \ 
      O=reference.dict

This will do the job.

• 0 views • link

written 5.5 years ago by tothepoint

Ok, I already created this dict file, but don't know this can serve as interval list, thanks.

• 0 views • link

written 5.5 years ago by MatthewP

Ok, this won't work, parse dict file to -L parameter of gatk command HaplotypeCaller will cause error.

• 0 views • link

written 5.5 years ago by MatthewP