Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-02-18T13:47:41.667Z
Answers.Text: 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.")
Answers.UserDisplayName: aechchiki
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-02-18T13:51:47.300Z
Answers.Text: 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
Answers.UserDisplayName: StupidWolf
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-02-18T13:58:54.303Z
Answers.Text: 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-
Answers.UserDisplayName: aechchiki
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-02-18T14:49:26.807Z
Answers.Text: 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:
Answers.UserDisplayName: aechchiki
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-02-18T15:03:15.387Z
Answers.Text: 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..
Answers.UserDisplayName: StupidWolf
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-02-18T13:25:38.093Z
Answers.LastActivityDate: 2020-02-18T13:25:38.093Z
Answers.UserDisplayName: StupidWolf
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-02-18T14:04:39.807Z
Answers.Text: 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.
Answers.UserDisplayName: aechchiki
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-02-18T14:12:39.593Z
Answers.Text: 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
Answers.UserDisplayName: Emanuel 
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-02-18T14:15:57.917Z
Answers.Text: 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"
Answers.UserDisplayName: aechchiki
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-02-18T13:58:32.173Z
Answers.LastActivityDate: 2020-02-18T13:58:32.173Z
Answers.UserDisplayName: Emanuel 
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-07T02:47:47.947Z
Answers.LastActivityDate: 2020-04-07T02:47:47.947Z
Answers.UserDisplayName: Geraldine_VdAuwera
Body: <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>
  
Tags: <file-formats><bed><gatk>
Title: Get start and end coordinates per chromosome
ThreadText: # Get start and end coordinates per chromosome
  
  (Asked by: username_0 on 2020-02-18)
  
  I will use GATK for SNP calling ([HaplotypeCaller](https://gatk.broadinstitute.org/hc/en-us/articles/360035531412-HaplotypeCaller-in-a-nutshell)). 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](https://gatkforums.broadinstitute.org/gatk/discussion/1319/collected-faqs-about-interval-lists) 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](https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/bam2bed.html))
  * 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...
  
  **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(<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")
  
  ```
  
  > **Comments:**
  > **username_0**: thanks for your input @username_1. 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 @username_1'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)
  
  
  ---
  
  