3
$\begingroup$

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
Improve this question
$\endgroup$
0

2 Answers 2

Reset to default
4
$\begingroup$

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
        }
    }
}
Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$
    – rseg
    Jan 12 at 0:29
  • $\begingroup$ @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. $\endgroup$
    – user438383
    Jan 13 at 10:58
7
$\begingroup$

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

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ 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. $\endgroup$
    – rseg
    Jan 12 at 0:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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