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

CC BY-SA 4.0
2023-01-12T00:23:07.667Z
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.
rseg
CC BY-SA 4.0
2023-01-11T07:08:40.970Z
2023-01-11T07:08:40.970Z
Steve
<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>

CC BY-SA 4.0
2023-01-12T00:29:59.637Z
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.
rseg
CC BY-SA 4.0
2023-01-13T10:58:48.673Z
@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.
user438383
CC BY-SA 4.0
2023-01-11T14:54:08.437Z
2023-01-11T15:02:26.347Z
2023-01-11T15:02:26.347Z
Ed Morton
<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>

<vcf><awk>
Extract lines from a file based on a range of coordinates listed on a second file
# 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

= " $2 " && $2 <= " $3 }' | awk -f - file1.vcf > output.vcf

I tried to modify the last awk code like this, but it didn't work:
= " $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; i2 <= 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**: @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.  

---

### Answer by username_1 on 2023-01-11. Score: 7

If you have a VCF file with an intact header, you can use 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 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 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 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.  

---

