2
$\begingroup$

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...

$\endgroup$
0
0
$\begingroup$

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")
$\endgroup$
5
  • $\begingroup$ thanks for your input @StupidWolf. yes I'm using the same reference, to which i also aligned all the reads. it helps the question, on how to get chromosome coordiantes - although not solving the GATK error, still getting "Input files reference and reads have incompatible contigs: No overlapping contigs found.") $\endgroup$ – aechchiki Feb 18 '20 at 13:47
  • $\begingroup$ Weird... are you providing the fasta file with the same header? grep ">" <your.fasta> and samtools view -H <bamfile> | grep "@SQ" ... you should see the same chromosome names $\endgroup$ – StupidWolf Feb 18 '20 at 13:51
  • $\begingroup$ yes, weird, agreed. and yes, i verified that they contain the same fasta IDs, which I also see in the bam header... I opened an issue on GATK forum, as self-comment to my question: gatk.broadinstitute.org/hc/en-us/community/posts/… $\endgroup$ – aechchiki Feb 18 '20 at 13:58
  • $\begingroup$ I upvoted your answer since it's what I was asking for in the original question. I'll deal with the GATK issue on my own. :up: $\endgroup$ – aechchiki Feb 18 '20 at 14:49
  • $\begingroup$ Thanks dude.. Can you make a small bam file and provide the fasta via ftp or something? Most likely the GATK people might ask for it. I cannot use GATK right now, but I can try with bcftools.. $\endgroup$ – StupidWolf Feb 18 '20 at 15:03
1
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ yes, I do have the genome, and have indexed it. as reported in the question itself, I could obtain (what I understood being) an interval file (just: chrID, start, end). and with @StupidWolf's answer I just tried to expand this interval to start from pos1 for each chr/scaf feed it in GATK. but I still get the same GATK errors (full link in my self-comment to my question) I got for the other interval files. weird, I agree. 100% sure I used one and only ref to map the reads to. I only have that one. $\endgroup$ – aechchiki Feb 18 '20 at 14:04
  • $\begingroup$ I think the reported error from GATK is clear, the annotation which should look as following: Bla_scaf1190 830 763 Is wrong as the start is higher then the stop. The question is why you would make an interval file which covers only mapped regions. Isnt the point of intervals to reduce the mapped space ? Otherwise speed should be identical for full space and your intervals $\endgroup$ – Emanuel Feb 18 '20 at 14:12
  • $\begingroup$ that's one test I tried, and that one error is clear. in the two other options (1.no intervals at all, technically taking the whole genome; 2.intervals for mapped regions of the chr; and I could add what I just tried, so 3.intervals for the whole chr length), I have the other error still persisting "Input files reference and reads have incompatible contigs: No overlapping contigs found" $\endgroup$ – aechchiki Feb 18 '20 at 14:15
0
$\begingroup$

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)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.