<p>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):</p>

<pre><code>samtools view -H some.bam | awk '{if(<span class="math-container">$1~/@SQ/){print $2"\t1\t"$</span>3}}' | sed 's/[SL]N://g' &gt; coords.bed
</code></pre>

<p>For the human genome I got this:</p>

<pre><code>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
</code></pre>

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

<p>If you use R, you can do:</p>

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

CC BY-SA 4.0
2020-02-18T13:47:41.667Z
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.")
aechchiki
CC BY-SA 4.0
2020-02-18T13:51:47.300Z
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
StupidWolf
CC BY-SA 4.0
2020-02-18T13:58:54.303Z
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: https://gatk.broadinstitute.org/hc/en-us/community/posts/360057671431-GATK-HaplotypeCaller-errors-with-or-without-the-L-interval-file-
aechchiki
CC BY-SA 4.0
2020-02-18T14:49:26.807Z
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:
aechchiki
CC BY-SA 4.0
2020-02-18T15:03:15.387Z
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..
StupidWolf
CC BY-SA 4.0
2020-02-18T13:25:38.093Z
2020-02-18T13:25:38.093Z
StupidWolf
<p>Not sure if I completely understood your question, but if you want simply an interval of your entire genome  you could do instead:</p>

<pre><code>samtools faidx MyGenome.fasta
</code></pre>

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

<pre><code>000000F 33203223        94      60      61
000001F 28828106        33756799        60      61
000002F 27810542        63065468        60      61
</code></pre>

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

CC BY-SA 4.0
2020-02-18T14:04:39.807Z
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.
aechchiki
CC BY-SA 4.0
2020-02-18T14:12:39.593Z
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
Emanuel 
CC BY-SA 4.0
2020-02-18T14:15:57.917Z
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"
aechchiki
CC BY-SA 4.0
2020-02-18T13:58:32.173Z
2020-02-18T13:58:32.173Z
Emanuel 
<p>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. <strong>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.</strong> 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).</p>

<p>See the GATK doc article for more details:
<a href="https://gatk.broadinstitute.org/hc/en-us/articles/360035531852-Intervals-and-interval-lists" rel="nofollow noreferrer">https://gatk.broadinstitute.org/hc/en-us/articles/360035531852-Intervals-and-interval-lists</a></p>

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

CC BY-SA 4.0
2020-04-07T02:47:47.947Z
2020-04-07T02:47:47.947Z
Geraldine_VdAuwera
<p>I will use GATK for SNP calling (<a href="https://gatk.broadinstitute.org/hc/en-us/articles/360035531412-HaplotypeCaller-in-a-nutshell" rel="nofollow noreferrer">HaplotypeCaller</a>). 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). </p>

<p>The <a href="https://gatkforums.broadinstitute.org/gatk/discussion/1319/collected-faqs-about-interval-lists" rel="nofollow noreferrer">interval file</a> for GATK can be for example a <code>bed</code> file, with 3 columns: <code>chr_name  chr_start  chr_end</code>. I don't have this file, but have the genome and the reads. </p>

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

<ul>
<li>converted the alignment to interval file, using <code>bam2bed</code>    (<a href="https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/bam2bed.html" rel="nofollow noreferrer">BEDOPS</a>)</li>
<li>extract the columns $1,2,3 (using <code>cut</code>)</li>
<li>get the min &amp; max coordinates per chr/scaff with a small R script:</li>
</ul>

<pre><code>bd &lt;- read.table('all.bed', h=FALSE)
bd_min &lt;- aggregate( bd$V3 ~ bd$V1 , bd, function(x) min(x))
bd_max &lt;- aggregate( bd$V2 ~ bd$V1 , bd, function(x) max(x))
bd_coord &lt;- 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)
</code></pre>

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

<p>The ideas for this:</p>

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

<p>The output looks plausible:</p>

<pre><code>$ cat coords.bed | head -3
Bla_chrm1   678 43860826
Bla_chrm10  181 20381540
Bla_chrm11  343 20367560
</code></pre>

<p>My question here:</p>

<ol>
<li>is this a correct way to proceed?</li>
<li>is there a standardized way to perform this?</li>
</ol>

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

<p>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. <code>awk '$3&gt;max[$1]{max[$1]=$3; row[$1]=$0} END{for (i in row) print row[i]}'</code>), but I get different results than with the R approach... </p>

<file-formats><bed><gatk>
Get start and end coordinates per chromosome
# Get start and end coordinates per chromosome

(Asked by: username_0 on 2020-02-18)

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:

is this a correct way to proceed?
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...

**Tags:** file-formats, bed, gatk

---

## 3 answers

### Answer by username_1 on 2020-02-18. Score: 0 | Accepted answer

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()[[1]]$targets
write.table(data.frame(names(CL),start=1,end=CL),
"coords.bed",row.names=FALSE,quote=F,col.names=FALSE,sep="\t")

> **Comments:**
> **username_0**: 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.")  
> **username_1**: 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  
> **username_0**: 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: https://gatk.broadinstitute.org/hc/en-us/community/posts/360057671431-GATK-HaplotypeCaller-errors-with-or-without-the-L-interval-file-  
> **username_0**: 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:  
> **username_1**: 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..  

---

### Answer by username_2 on 2020-02-18. Score: 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.

> **Comments:**
> **username_0**: 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.  
> **username_2**: 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  
> **username_0**: 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"  

---

### Answer by username_3 on 2020-04-07. Score: 0

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)


---

