1

I have some .vcf files and I want to filter some variants out. This is just small part of my file: there are some header lines at the beginning of the file (starting with ##) and then variants (one row per variant).

##fileformat=VCFv4.2 ##source=combiSV-v2.2 ##fileDate=Mon May 8 11:32:53 2023 ##contig=<ID=chrM,length=16571> ##contig=<ID=chr1,length=249250621> ##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record"> ##INFO=<ID=SVCALLERS,Number=.,Type=String,Description="SV callers that support this SV"> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##FORMAT=<ID=DR,Number=1,Type=Integer,Description="# High-quality reference reads"> ##FORMAT=<ID=DV,Number=1,Type=Integer,Description="# High-quality variant reads"> #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample 1 10862 id.1 N <INS> . PASS SVTYPE=INS;SVLEN=101;END=10862;SVCALLERS=cutesv,SVIM GT:DR:DV 1/1:0:26 1 90258 id.2 N <INS> . PASS SVTYPE=INS;SVLEN=118;END=90258;SVCALLERS=SVIM,NanoSV GT:DR:DV 1/1:0:9 1 90259 id.3 N <INS> . PASS SVTYPE=INS;SVLEN=36;END=90259;SVCALLERS=Sniffles GT:DR:DV 0/1:44:7 1 185824 id.4 N <DEL> . PASS SVTYPE=DEL;SVLEN=80;END=186660;SVCALLERS=Sniffles,cutesv GT:DR:DV 1/1:0:15 1 186241 id.5 N <DEL> . PASS SVTYPE=DEL;SVLEN=418;END=186662;SVCALLERS=SVIM,NanoSV GT:DR:DV 1/1:2:12 1 526111 id.6 N <DEL> . PASS SVTYPE=DEL;SVLEN=624;END=526735;SVCALLERS=Sniffles,cutesv GT:DR:DV 0/1:8 2 91926078 id.3958 N <BND> . PASS SVTYPE=BND;SVLEN=.;END=;SVCALLERS=Sniffles,NanoSV GT:DR:DV 0/1:60:15

While keeping the header lines, I want to remove rows with SVLEN < 100 and those with only one SVCALLERS included. These are two criteria that both must be met, in other words I want to keep only rows with SVLEN > 100 and at least two SVCALLERS).

In addition there are some rows where ALT is BND and the file does not provide any SVLEN for this type of variant, so if the row contains BND, I just want to keep it if it is supported by two callers.

Examples: I want to drop this variant because SVLEN is less than 100 and only one SVCALLERS detected it:

SVTYPE=INS;SVLEN=36;END=90259;SVCALLERS=Sniffles GT:DR:DV 0/1:44:7 1 185824 id.4 N <DEL> . PASS

Or this row as well, although there are two callers but SVLEN is less than 100:

SVTYPE=DEL;SVLEN=80;END=186660;SVCALLERS=Sniffles,cutesv GT:DR:DV 1/1:0:15 1 186241 id.5 N <DEL> . PASS

Is there an easy way to do it? Thanks

My final file should look like this:

##fileformat=VCFv4.2 ##source=combiSV-v2.2 ##fileDate=Mon May 8 11:32:53 2023 ##contig=<ID=chrM,length=16571> ##contig=<ID=chr1,length=249250621> ##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record"> ##INFO=<ID=SVCALLERS,Number=.,Type=String,Description="SV callers that support this SV"> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##FORMAT=<ID=DR,Number=1,Type=Integer,Description="# High-quality reference reads"> ##FORMAT=<ID=DV,Number=1,Type=Integer,Description="# High-quality variant reads"> #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample 1 10862 id.1 N <INS> . PASS SVTYPE=INS;SVLEN=101;END=10862;SVCALLERS=cutesv,SVIM GT:DR:DV 1/1:0:26 1 90258 id.2 N <INS> . PASS SVTYPE=INS;SVLEN=118;END=90258;SVCALLERS=SVIM,NanoSV GT:DR:DV 1/1:0:9 1 186241 id.5 N <DEL> . PASS SVTYPE=DEL;SVLEN=418;END=186662;SVCALLERS=SVIM,NanoSV GT:DR:DV 1/1:2:12 1 526111 id.6 N <DEL> . PASS SVTYPE=DEL;SVLEN=624;END=526735;SVCALLERS=Sniffles,cutesv GT:DR:DV 0/1:8 2 91926078 id.3958 N <BND> . PASS SVTYPE=BND;SVLEN=.;END=;SVCALLERS=Sniffles,NanoSV GT:DR:DV 0/1:60:15
Improve this question
3
  • hey, since you're already settled on awk: Do you have an existing approach that you've already worked on and are stuck at? Commented May 18, 2023 at 17:32
  • Hey, not necessarily, I am just trying to find a way and I thought awk can be a way to achieve it Commented May 18, 2023 at 17:35
  • That's a LOT of requirements/checks to be asking for help with. Why not ask for help with a couple of checks and then try adding on to it to implement the rest yourself? Right now it's like you're asking someone to develop your whole project for free. Commented May 18, 2023 at 21:00

2 Answers 2

Reset to default
7

Here's a perl way:

$ perl -F'\t' -lane ' if(/^#/){ print; next }; $F[7] =~ /\bSVLEN=(\d+)/; $svlen=$1; $F[7] =~ /\bSVCALLERS=([^;]+)/; @callers=split(/,/,$1); print if $svlen > 100 && scalar(@callers) > 1' file.vcf ##fileformat=VCFv4.2 ##source=combiSV-v2.2 ##fileDate=Mon May 8 11:32:53 2023 ##contig=<ID=chrM,length=16571> ##contig=<ID=chr1,length=249250621> ##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record"> ##INFO=<ID=SVCALLERS,Number=.,Type=String,Description="SV callers that support this SV"> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##FORMAT=<ID=DR,Number=1,Type=Integer,Description="# High-quality reference reads"> ##FORMAT=<ID=DV,Number=1,Type=Integer,Description="# High-quality variant reads"> #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample 1 10862 id.1 N <INS> . PASS SVTYPE=INS;SVLEN=101;END=10862;SVCALLERS=cutesv,SVIM GT:DR:DV 1/1:0:26 1 90258 id.2 N <INS> . PASS SVTYPE=INS;SVLEN=118;END=90258;SVCALLERS=SVIM,NanoSV GT:DR:DV 1/1:0:9 1 186241 id.5 N <DEL> . PASS SVTYPE=DEL;SVLEN=418;END=186662;SVCALLERS=SVIM,NanoSV GT:DR:DV 1/1:2:12 1 526111 id.6 N <DEL> . PASS SVTYPE=DEL;SVLEN=624;END=526735;SVCALLERS=Sniffles,cutesv GT:DR:DV 0/1:8

Explanation

  • perl -F'\t' -lane: the -a makes perl work kinda like awk in that it automatically splits each input line on the character given by -F (whitespace, by default, like awk) and saves it in the array @F. So, since arrays start counting at 0, the first field is $F[0[, the second $F[1] and so on. Next, the -l removes trailing newlines from each input line, and also adds a \n to each print call, the -n means "read the input file line by line and apply the script given by -e to each line.
  • if(/^#/){ print; next }; : if this is a header line, just print it and move on to the next line.
  • $F[7] =~ /\bSVLEN=(\d+)/; $svlen=$1;: match the longest string of numbers after SVLEN= in the 8th field and save it as $svlen. This will be the length. The \b ensures we only match at word boundaries so this won't fail if you have something like NOTSVLEN= in your file.
  • $F[7] =~ /\bSVCALLERS=([^;]+)/; @callers=split(/,/,$1);: now search the 8th field for the string SVCALLERS=, take the longest stretch of non-; characters after it, and then split that on , into the array @callers. This now has the list of CNV callers used for this CNV.
  • print if $svlen > 100 && scalar(@callers) > 1: we now print the line if the length was more than 100 and the number of callers (scalar(@array) gives the number of elements in an array) is more than 1.

And here's the same basic thing in a more concise and less clear way if you prefer your commands mildly golfed:

perl -F'\t' -lane '$F[7]=~/\bSVLEN=(\d+)/;$s=$1;$F[7]=~/\bSVCALLERS=([^;]+)/; /^#/ || ($s>100&&scalar(split(/,/,$1)) > 1) || next; print' file.vcf

If you also want to keep lines with no SVLEN as long as they have at least two callers, use this:

$ perl -F'\t' -lane 'if(/^#/){ print; next }; $F[7] =~ /\bSVLEN=([.\d]+)/; $svlen=$1; $F[7] =~ /\bSVCALLERS=([^;]+)/; next unless ($svlen > 100 || $svlen == ".") && scalar(split(/,/,$1)) > 1; print' file.vcf ##fileformat=VCFv4.2 ##source=combiSV-v2.2 ##fileDate=Mon May 8 11:32:53 2023 ##contig=<ID=chrM,length=16571> ##contig=<ID=chr1,length=249250621> ##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record"> ##INFO=<ID=SVCALLERS,Number=.,Type=String,Description="SV callers that support this SV"> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##FORMAT=<ID=DR,Number=1,Type=Integer,Description="# High-quality reference reads"> ##FORMAT=<ID=DV,Number=1,Type=Integer,Description="# High-quality variant reads"> #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample 1 10862 id.1 N <INS> . PASS SVTYPE=INS;SVLEN=101;END=10862;SVCALLERS=cutesv,SVIM GT:DR:DV 1/1:0:26 1 90258 id.2 N <INS> . PASS SVTYPE=INS;SVLEN=118;END=90258;SVCALLERS=SVIM,NanoSV GT:DR:DV 1/1:0:9 1 186241 id.5 N <DEL> . PASS SVTYPE=DEL;SVLEN=418;END=186662;SVCALLERS=SVIM,NanoSV GT:DR:DV 1/1:2:12 1 526111 id.6 N <DEL> . PASS SVTYPE=DEL;SVLEN=624;END=526735;SVCALLERS=Sniffles,cutesv GT:DR:DV0/1:8 2 91926078 id.3958 N <BND> . PASS SVTYPE=BND;SVLEN=.;END=;SVCALLERS=Sniffles,NanoSV GT:DR:DV0/1:60:15
Improve this answer
2
  • Hey@terdon, nice this indeed does what I want. However I noticed I have some rows which are translocations and the vcf file does not provide ant SVLEN. (they are shown as: SVLEN=.). So if I face with such a variant I just need to meet the second criteria; variant confirmed by two caller. Is there any way to adjust it in your code? I just updated my question based on this. Sorry about that, totally forgot about it Commented May 18, 2023 at 19:14
  • @Anna1364 see update. Note that structural variant VCF files can also have <CNV> in the ALT field, and some other values, I think, they can be imprecise in which case they will have no length, and various other issues. Commented May 18, 2023 at 21:57
1

This ignores the header lines, just handles the data lines.

start cmd:> awk -F '=|;| +' '$11<100 || $15 !~ "," { next; }; { print $0; }' input 1 10862 id.1 N <INS> . PASS SVTYPE=INS;SVLEN=101;END=10862;SVCALLERS=cutesv,SVIM GT:DR:DV 1/1:0:26 1 90258 id.2 N <INS> . PASS SVTYPE=INS;SVLEN=118;END=90258;SVCALLERS=SVIM,NanoSV GT:DR:DV 1/1:0:9 1 186241 id.5 N <DEL> . PASS SVTYPE=DEL;SVLEN=418;END=186662;SVCALLERS=SVIM,NanoSV GT:DR:DV 1/1:2:12 1 526111 id.6 N <DEL> . PASS SVTYPE=DEL;SVLEN=624;END=526735;SVCALLERS=Sniffles,cutesv GT:DR:DV 0/1:8

Update 1

For the adjusted question

start cmd:> awk -F '=|;| +' '/^#/ { print; next; }; $5 != "<BND>" && $11<100 || $15 !~ "," { next; }; $5 == "<BND>" && $15 !~ "," { next; }; { print $0; }' input ##fileformat=VCFv4.2 ##source=combiSV-v2.2 ##fileDate=Mon May 8 11:32:53 2023 ##contig=<ID=chrM,length=16571> ##contig=<ID=chr1,length=249250621> ##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record"> ##INFO=<ID=SVCALLERS,Number=.,Type=String,Description="SV callers that support this SV"> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##FORMAT=<ID=DR,Number=1,Type=Integer,Description="# High-quality reference reads"> ##FORMAT=<ID=DV,Number=1,Type=Integer,Description="# High-quality variant reads"> #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample 1 10862 id.1 N <INS> . PASS SVTYPE=INS;SVLEN=101;END=10862;SVCALLERS=cutesv,SVIM GT:DR:DV 1/1:0:26 1 90258 id.2 N <INS> . PASS SVTYPE=INS;SVLEN=118;END=90258;SVCALLERS=SVIM,NanoSV GT:DR:DV 1/1:0:9 1 186241 id.5 N <DEL> . PASS SVTYPE=DEL;SVLEN=418;END=186662;SVCALLERS=SVIM,NanoSV GT:DR:DV 1/1:2:12 1 526111 id.6 N <DEL> . PASS SVTYPE=DEL;SVLEN=624;END=526735;SVCALLERS=Sniffles,cutesv GT:DR:DV 0/1:8 2 91926078 id.3958 N <BND> . PASS SVTYPE=BND;SVLEN=.;END=;SVCALLERS=Sniffles,NanoSV GT:DR:DV 0/1:60:15
Improve this answer
3
  • I slightly updated my question, there was one missing point. Also I need to keep the headers. Commented May 18, 2023 at 19:24
  • @Anna1364 See the update Commented May 18, 2023 at 19:41
  • Thank so much for proving your answer as well, that also worked. Commented May 18, 2023 at 22:11

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.