1

I have a .vcf file of 138 first header lines (starts with #) and the others of data (snp in rows(322045) and patients with some infos in columns(the first 10). I used a script bash to calculate for each row the number of cells in that row that are different (in the initial part) from "0|0": this is my script

for j in {139..322045}

    do
     c=0
     awk -v var=$c -v j=$j 'NR==j{for(i=10; i<=NF; i++) {if(substr($i,1,3)!="0|0") var++}} END{ print $1 ":" $2 "\t" var }' file.vcf >> out.txt
    done

This is the input:

> #<info>
> #..
> # . . . 
21  9411245  x  C   A   505  PASS   AC=2   GT:AD:DP:GQ:PL   0|0:11  0|0:12
21  9411246  y  C   T   505  PASS   AC=2   GT:AD:DP:GQ:PL   0|0:11  1|0:13

(the columns are tab separated) Then I print the 1th and 2th columns linked by : and the count; but it does not work completely, if I use a subset containing only 2 row it works perfectly. This is the result

21:48111872 2
21:48111872 1
21:48111872 0
21:48111872 2

It repeats the rows

How can I fix it ? Thanks in advance, if you fix it please write a short explanation of that.

N.B. it takes a lot of time to compute it (also using for {139..160})

1 Answer 1

1

The reason it isn't working is that you are printing $1 and $2 in the END{} block. The END{} is only run once, after the last line of the input file has been read. So $1 and $2 will always be the 1st and 2nd fields of the last line.

In any case, this is an incredibly inefficient way of parsing a text file. You are reading the entire thing for every iteration of the loop. And shell loops are very slow. So you are using a very slow loop and you are needlessly reading thousands of lines in your awk over and over again.

Instead of using a shell loop, just do everything in awk:

$ awk -F"\t" '/^[^#]/{var=0; for(i=10; i<=NF; i++) {if(substr($i,1,3)!="0|0") var++} print $1 ":" $2 "\t" var }' foo.vcf 
21:9411245  0
21:9411246  1

Or, a little less condensed:

awk -F"\t" '/^[^#]/{
        var=0; 
        for(i=10; i<=NF; i++) {
            if(substr($i,1,3)!="0|0"){
                var++
            } 
        }
        print $1 ":" $2 "\t" var 
    }' foo.vcf 

Explanation

  • -F"\t" : set the input field separator to a tab.
  • /^[^#]/{ ... }: do this only for lines that start (/^a/ will match lines starting with a) with a character that isn't a # ([^#]).
  • var=0; : set var back to 0 for every input line.
  • for(i=10; i<=NF; i++) {if(substr($i,1,3)!="0|0") var++} : this is your original code, it counts the number of times it finds a genotype that isn't 0|0.
  • print $1 ":" $2 "\t" var : again, your code, but now outside the END{} block so it is run on every line and not only the end.

That's it. No shell loop needed and it should take just a few seconds.

2
  • Ok thanks, I don't know very well awk. So in this way whit "NF" I can read the all .vcf files for both columns and rows, right ? Because NF is the "cell" ?
    – Will
    Nov 12, 2019 at 23:08
  • @Will NF is the "number of fields". So if a line has 5 fields (cells), then NF will be 5. By the way, for future questions, you might want to check out our sister site: Bioinformatics. People here won't necessarily know what VCF files or SNPs are, and there are sometimes important details about the format that can affect the answers. Asking there is safer for this sort of format.
    – terdon
    Nov 12, 2019 at 23:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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