Get start and end coordinates per chromosome

Question

I will use GATK for SNP calling (HaplotypeCaller). I need to feed the interval file in the command, otherwise I get errors (even though I want to use the whole genome, not a subset - btw it's not the point of the question, so will not elaborate).

The interval file for GATK can be for example a bed file, with 3 columns: chr_name chr_start chr_end. I don't have this file, but have the genome and the reads.

To obtain the intervals per chr/scaf, I proceeded this way:

• converted the alignment to interval file, using bam2bed (BEDOPS)
• extract the columns $1,2,3 (using cut)
• get the min & max coordinates per chr/scaff with a small R script:

bd <- read.table('all.bed', h=FALSE)
bd_min <- aggregate( bd$V3 ~ bd$V1 , bd, function(x) min(x))
bd_max <- aggregate( bd$V2 ~ bd$V1 , bd, function(x) max(x))
bd_coord <- merge(bd_min, bd_max, by='bd$V1')
write.table(bd_coord, file='coords.bed', sep='\t', quote=FALSE, col.names=FALSE, row.names=FALSE)

These should thus represent the mapped intervals, conceptually same as not feeding the intervals at all (= considering all the sequences).

The ideas for this:

• in the bed file: $1 is be the feature, $2 the start position, $3 the end position
• so, I need to get the min of $2 for each unique item in $1(for the start coord) and the max of $3 for each item in $1 (for the end coord)

The output looks plausible:

$ cat coords.bed | head -3
Bla_chrm1   678 43860826
Bla_chrm10  181 20381540
Bla_chrm11  343 20367560

My question here:

1. is this a correct way to proceed?
2. is there a standardized way to perform this?

My main concern is that I will have to re-build this "coordinate" files for every GATK run, because I think the coordinates per chr/scaf will shift, even slightly, at each GATK run with different datasets.

My other minor concern, is that this approach is a bit slow (given the 160M lines in the current bed file), so an unix-tools solution will also be accepted. I tried to compute myself the max/min using awk (e.g. awk '$3>max[$1]{max[$1]=$3; row[$1]=$0} END{for (i in row) print row[i]}'), but I get different results than with the R approach...

Answer 1

You can try this, since you have the bam file, all of them are aligned to the same reference genome (I hope)(and apologies for the awful awk.. been a while since I used it):

samtools view -H some.bam | awk '{if($1~/@SQ/){print $2"\t1\t"$3}}' | sed 's/[SL]N://g' > coords.bed

For the human genome I got this:

chr1    1   195471971
chr2    1   182113224
chr3    1   160039680
chr4    1   156508116
chr5    1   151834684
chr6    1   149736546
chr7    1   145441459
chr8    1   129401213
chr9    1   124595110

What I did is to get the chromosome length, and make every line start at 1.

If you use R, you can do:

library(Rsamtools)
CL = scanBamHeader(<bam file>)[[1]]$targets
write.table(data.frame(names(CL),start=1,end=CL),
"coords.bed",row.names=FALSE,quote=F,col.names=FALSE,sep="\t")

Answer 2

Not sure if I completely understood your question, but if you want simply an interval of your entire genome you could do instead:

samtools faidx MyGenome.fasta

Obviously that means that you need access to that, not sure if you have. The resulting fai file contains a format similar to that:

000000F 33203223        94      60      61
000001F 28828106        33756799        60      61
000002F 27810542        63065468        60      61

Where the first column is your chromosome and the second one your length. I guess you could then just do a awk '{OFS=FS="\t"}{print $1,"0",$2-1}' to get your file.

Answer 3

This probably comes too late, but for the record GATK supports multiple kinds of interval file format. BED is the most common, and there's the Picard format (.interval_list) which is BED-like with a few tweaks. Both require coordinates. However, GATK also supports a simpler format (.list or .intervals) that allows you to just give the contig names (eg chr1) without any coordinates. Sounds like that would satisfy your need. Incidentally, GATK also supports using VCF files as interval file if you want to re-genotype specific regions based on variants located there (ideally paired with the interval padding arguments).

See the GATK doc article for more details: https://gatk.broadinstitute.org/hc/en-us/articles/360035531852-Intervals-and-interval-lists

(The doc has a few formatting issues but I just notified the team so they'll fix it)