Answers.Body: <p>If you have a VCF file with an intact header, you can use <a href="/questions/tagged/bcftools" class="post-tag" title="show questions tagged &#39;bcftools&#39;" aria-label="show questions tagged &#39;bcftools&#39;" rel="tag" aria-labelledby="bcftools-container">bcftools</a> for this. All you need to do first is to bgzip and index your VCF:</p>
  <pre><code>$ bgzip file.vcf
  $ bcftools index file.vcf.gz
  </code></pre>
  <p>Then, assuming your list of regions are actually 1-based tab-delimited, use <a href="https://samtools.github.io/bcftools/bcftools.html#view" rel="noreferrer"><code>bcftools view</code></a> and the <code>--regions-file</code> option to restrict the search to the regions listed. Use also the <code>--no-header</code> option to suppress the header in the output:</p>
  <pre><code>bcftools view -HR regions.txt file.vcf.gz
  </code></pre>
  <p>From the <a href="https://samtools.github.io/bcftools/bcftools.html#common_options" rel="noreferrer">common options</a> documentation:</p>
  <blockquote>
  <p>-R, --regions-file FILE</p>
  <p>Regions can be specified either on command line or in a VCF, BED, or tab-delimited file (the default). The
  columns of the tab-delimited file can contain either positions
  (two-column format) or intervals (three-column format): CHROM, POS,
  and, optionally, END, where positions are 1-based and inclusive. The
  columns of the tab-delimited BED file are also CHROM, POS and END
  (trailing columns are ignored), but coordinates are 0-based,
  half-open. To indicate that a file be treated as BED rather than the
  1-based tab-delimited file, the file must have the &quot;.bed&quot; or &quot;.bed.gz&quot;
  suffix (case-insensitive). Uncompressed files are stored in memory,
  while bgzip-compressed and tabix-indexed region files are streamed.
  Note that sequence names must match exactly, &quot;chr20&quot; is not the same
  as &quot;20&quot;. Also note that chromosome ordering in FILE will be respected,
  the VCF will be processed in the order in which chromosomes first
  appear in FILE. However, within chromosomes, the VCF will always be
  processed in ascending genomic coordinate order no matter what order
  they appear in FILE. Note that overlapping regions in FILE can result
  in duplicated out of order positions in the output. This option
  requires indexed VCF/BCF files. Note that -R cannot be used in
  combination with -r.</p>
  </blockquote>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-12T00:23:07.667Z
Answers.Text: Hi Steve, I converted my list file to tab delimited and following your suggestions bcftools worked great doing this job, and it's very fast too. It's handy to have a dedicated tool to parse vcf files. Thanks for your help.
Answers.UserDisplayName: rseg
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-11T07:08:40.970Z
Answers.LastActivityDate: 2023-01-11T07:08:40.970Z
Answers.UserDisplayName: Steve
Answers.Body: <p>Using GNU awk for arrays of arrays:</p>
  <pre><code><span class="math-container">$ cat tst.awk
  NR==FNR {
      ranges[$1][$2 FS $3]
      next
  }
  FNR == 1 {
      print
      next
  }
  $</span>1 in ranges {
      for ( range in ranges[<span class="math-container">$1] ) {
          split(range,r)
          if ( (r[1] &lt;= $2) &amp;&amp; ($</span>2 &lt;= r[2]) ) {
              print
              next
          }
      }
  }
  </code></pre>
  <p>
  <pre><code>$ awk -f tst.awk listfile.txt file1.vcf
  #CHROM      POS ID  REF ALT QUAL    FILTER  INFO    FORMAT          COUNTS
  scaffold_1  210 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  340 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  370 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  800 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2      850     .       C       CTA     40      PASS    .       GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  </code></pre>
  <p>Using any awk it'd be:</p>
  <pre><code><span class="math-container">$ cat tst.awk
  NR==FNR {
      ranges[$1] = ranges[$1] $2 FS $3 RS
      next
  }
  FNR == 1 {
      print
      next
  }
  $</span>1 in ranges {
      n = split(ranges[<span class="math-container">$1],list,RS)
      for ( i=1; i&lt;n; i++ ) {
          range = list[i]
          split(range,r)
          if ( (r[1] &lt;= $2) &amp;&amp; ($</span>2 &lt;= r[2]) ) {
              print
              next
          }
      }
  }
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-12T00:29:59.637Z
Answers.Text: Hi Ed, your awk code works very well doing this job. Very clever solution. Thanks again for your timely help. That's right, the question I posted in the usual Unix page was moved to this site by the administrator.
Answers.UserDisplayName: rseg
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-13T10:58:48.673Z
Answers.Text: @rseg whilst I am sure this is a great solution (ed definitely knows awk well!), I would strongly advise to use dedicated tools like bcftools. It is much easier to use, less prone to errors as they will have accounted for any of the peculiarities of vcf files, and is more robust in a pipeline. There's no need to use awk for such a simple example.
Answers.UserDisplayName: user438383
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-11T14:54:08.437Z
Answers.LastActivityDate: 2023-01-11T15:02:26.347Z
Answers.LastEditDate: 2023-01-11T15:02:26.347Z
Answers.UserDisplayName: Ed Morton
Body: <p>I would be highly appreciated to have a script able to perform the following operation:</p>
  <p>Only extract scaffold_1 lines from file1.vcf in which POS number (column 2) falls within column 2 and column 3 numbers of scaffold_1 lines from listfile.txt (inclusive). Similarly, only extract scaffold_2 lines from file1.vcf in which POS number falls within column 2 and column 3 numbers of scaffold_2 lines from listfile.txt. In reality the vcf I have continues until scaffold_12 and has 17,000 lines.</p>
  <p>file1.vcf</p>
  <pre><code>#CHROM      POS ID  REF ALT QUAL    FILTER  INFO    FORMAT          COUNTS
  scaffold_1  160 .   C   A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  210 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  340 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  470 .   T   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  580 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  630 .   CA  C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  40  .   C   A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  100 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  370 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  600 .   T   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  800 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2      850     .       C       CTA     40      PASS    .       GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  1120    .   CA  C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  </code></pre>
  <p>listfile.txt</p>
  <pre><code>scaffold_1  50  100
  scaffold_1  200 350
  scaffold_2  350 450
  scaffold_2  750 850
  </code></pre>
  <p>The desired output file is the following:</p>
  <p>desired_file.vcf</p>
  <pre><code>#CHROM      POS ID  REF ALT QUAL    FILTER  INFO    FORMAT          COUNTS
  scaffold_1  210 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  340 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  370 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  800 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2      850     .       C       CTA     40      PASS    .       GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  </code></pre>
  <p>To do this job I found the below awk code from another similar post but they don't differentiate scaffold_1 from scaffold_2, and I don't how to add this extra variable properly. Thank you.</p>
  <pre><code>awk 'NR==FNR{rng++;start[rng]=$2;end[rng]=$3;next} {for (i=1;i&lt;=rng;i++) if (($2&gt;=start[i])&amp;&amp;($2&lt;=end[i])) {print; next}}' listfile.txt file1.vcf &gt; output.vcf
  </code></pre>
  <pre><code>&lt;listfile.txt awk '{ print &quot;$2 &gt;= &quot; $2 &quot; &amp;&amp; $2 &lt;= &quot; $3 }' | awk -f - file1.vcf &gt; output.vcf
  </code></pre>
  <p>I tried to modify the last awk code like this, but it didn't work:</p>
  <pre><code>&lt;listfile.txt awk '{ print &quot;$1 == &quot; $1 &quot; &amp;&amp; $2 &gt;= &quot; $2 &quot; &amp;&amp; $2 &lt;= &quot; $3 }' | awk -f - file1.vcf &gt; output.vcf
  </code></pre>
  
Tags: <vcf><awk>
Title: Extract lines from a file based on a range of coordinates listed on a second file
ThreadText: # Extract lines from a file based on a range of coordinates listed on a second file
  
  (Asked by: username_0 on 2023-01-10)
  
  I would be highly appreciated to have a script able to perform the following operation:
  
  Only extract scaffold\_1 lines from file1.vcf in which POS number (column 2) falls within column 2 and column 3 numbers of scaffold\_1 lines from listfile.txt (inclusive). Similarly, only extract scaffold\_2 lines from file1.vcf in which POS number falls within column 2 and column 3 numbers of scaffold\_2 lines from listfile.txt. In reality the vcf I have continues until scaffold\_12 and has 17,000 lines.
  
  file1.vcf
  
  ```
  #CHROM      POS ID  REF ALT QUAL    FILTER  INFO    FORMAT          COUNTS
  scaffold_1  160 .   C   A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  210 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  340 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  470 .   T   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  580 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  630 .   CA  C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  40  .   C   A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  100 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  370 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  600 .   T   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  800 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2      850     .       C       CTA     40      PASS    .       GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  1120    .   CA  C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  
  ```
  
  listfile.txt
  
  ```
  scaffold_1  50  100
  scaffold_1  200 350
  scaffold_2  350 450
  scaffold_2  750 850
  
  ```
  
  The desired output file is the following:
  
  desired\_file.vcf
  
  ```
  #CHROM      POS ID  REF ALT QUAL    FILTER  INFO    FORMAT          COUNTS
  scaffold_1  210 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  340 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  370 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  800 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2      850     .       C       CTA     40      PASS    .       GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  
  ```
  
  To do this job I found the below awk code from another similar post but they don't differentiate scaffold\_1 from scaffold\_2, and I don't how to add this extra variable properly. Thank you.
  
  ```
  awk 'NR==FNR{rng++;start[rng]=$2;end[rng]=$3;next} {for (i=1;i<=rng;i++) if (($2>=start[i])&&($2<=end[i])) {print; next}}' listfile.txt file1.vcf > output.vcf
  
  ```
  
  ```
  <listfile.txt awk '{ print "$2 >= " $2 " && $2 <= " $3 }' | awk -f - file1.vcf > output.vcf
  
  ```
  
  I tried to modify the last awk code like this, but it didn't work:
  
  ```
  <listfile.txt awk '{ print "$1 == " $1 " && $2 >= " $2 " && $2 <= " $3 }' | awk -f - file1.vcf > output.vcf
  
  ```
  
  **Tags:** vcf, awk
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_2 on 2023-01-11. Score: 4 | Accepted answer
  
  Using GNU awk for arrays of arrays:
  
  ```
  $ cat tst.awk
  NR==FNR {
      ranges[$1][$2 FS $3]
      next
  }
  FNR == 1 {
      print
      next
  }
  $1 in ranges {
      for ( range in ranges[$1] ) {
          split(range,r)
          if ( (r[1] <= $2) && ($2 <= r[2]) ) {
              print
              next
          }
      }
  }
  
  ```
  
  ```
  $ awk -f tst.awk listfile.txt file1.vcf
  #CHROM      POS ID  REF ALT QUAL    FILTER  INFO    FORMAT          COUNTS
  scaffold_1  210 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_1  340 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  370 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2  800 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  scaffold_2      850     .       C       CTA     40      PASS    .       GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
  
  ```
  
  Using any awk it'd be:
  
  ```
  $ cat tst.awk
  NR==FNR {
      ranges[$1] = ranges[$1] $2 FS $3 RS
      next
  }
  FNR == 1 {
      print
      next
  }
  $1 in ranges {
      n = split(ranges[$1],list,RS)
      for ( i=1; i<n; i++ ) {
          range = list[i]
          split(range,r)
          if ( (r[1] <= $2) && ($2 <= r[2]) ) {
              print
              next
          }
      }
  }
  
  ```
  
  > **Comments:**
  > **username_0**: Hi Ed, your awk code works very well doing this job. Very clever solution. Thanks again for your timely help. That's right, the question I posted in the usual Unix page was moved to this site by the administrator.  
  > **username_3**: @username_0 whilst I am sure this is a great solution (ed definitely knows awk well!), I would strongly advise to use dedicated tools like bcftools. It is much easier to use, less prone to errors as they will have accounted for any of the peculiarities of vcf files, and is more robust in a pipeline. There's no need to use awk for such a simple example.  
  
  ---
  
  ### Answer by username_1 on 2023-01-11. Score: 7
  
  If you have a VCF file with an intact header, you can use [bcftools](/questions/tagged/bcftools "show questions tagged 'bcftools'") for this. All you need to do first is to bgzip and index your VCF:
  
  ```
  $ bgzip file.vcf
  $ bcftools index file.vcf.gz
  
  ```
  
  Then, assuming your list of regions are actually 1-based tab-delimited, use [`bcftools view`](https://samtools.github.io/bcftools/bcftools.html#view) and the `--regions-file` option to restrict the search to the regions listed. Use also the `--no-header` option to suppress the header in the output:
  
  ```
  bcftools view -HR regions.txt file.vcf.gz
  
  ```
  
  From the [common options](https://samtools.github.io/bcftools/bcftools.html#common_options) documentation:
  
  > -R, --regions-file FILE
  >
  > Regions can be specified either on command line or in a VCF, BED, or tab-delimited file (the default). The
  > columns of the tab-delimited file can contain either positions
  > (two-column format) or intervals (three-column format): CHROM, POS,
  > and, optionally, END, where positions are 1-based and inclusive. The
  > columns of the tab-delimited BED file are also CHROM, POS and END
  > (trailing columns are ignored), but coordinates are 0-based,
  > half-open. To indicate that a file be treated as BED rather than the
  > 1-based tab-delimited file, the file must have the ".bed" or ".bed.gz"
  > suffix (case-insensitive). Uncompressed files are stored in memory,
  > while bgzip-compressed and tabix-indexed region files are streamed.
  > Note that sequence names must match exactly, "chr20" is not the same
  > as "20". Also note that chromosome ordering in FILE will be respected,
  > the VCF will be processed in the order in which chromosomes first
  > appear in FILE. However, within chromosomes, the VCF will always be
  > processed in ascending genomic coordinate order no matter what order
  > they appear in FILE. Note that overlapping regions in FILE can result
  > in duplicated out of order positions in the output. This option
  > requires indexed VCF/BCF files. Note that -R cannot be used in
  > combination with -r.
  
  > **Comments:**
  > **username_0**: Hi username_1, I converted my list file to tab delimited and following your suggestions bcftools worked great doing this job, and it's very fast too. It's handy to have a dedicated tool to parse vcf files. Thanks for your help.  
  
  ---
  
  