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

Question

[MatthewP]

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.

[Manuel Sokolov Ravasqueira]

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!

Answer 1

[MatthewP]

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

[tomas4482]

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.

Answer 2

[tothepoint]

Use picard tool

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

This will do the job.

[MatthewP]

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

[MatthewP]

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