CODEWITHSUNDEEP
bashawk
The Nightman
The Nightman

Reputation: 5759

Identifying columns from text file with awk

I am analyzing some genomic sequencing data, and I am trying to use my script to look for lines starting 'chr' and then grab the values in columns 1,2,4, and 5 then search for the numbers corresponding to RO and AO and calculate the frequency of allele variant by AO/(AO+RO).

My script works much of the time when I go back and manually check it by calculating AO/(AO+RO). But sometimes like for the following read it fails.

For the following read, my script calculates that AO/(AO+RO) is 0.214286, but I calculate it to be 3/(3+7757)=0.000387. If someone can see why this is I would appreciate it.

chr4    110541280   .   G   A   0   .   AB=0;ABP=0;AC=0;AF=0;AN=2;AO=3;CIGAR=1X;DP=7760;DPB=7760;DPRA=0;EPP=9.52472;EPPR=16838.4;GTI=0;LEN=1;MEANALT=1;MQM=60;MQMR=60;NS=1;NUMALT=1;ODDS=10694.5;PAIRED=0;PAIREDR=0;PAO=0;PQA=0;PQR=0;PRO=0;QA=111;QR=280123;RO=7757;RPP=9.52472;RPPR=16838.4;RUN=1;SAF=3;SAP=9.52472;SAR=0;SRF=7757;SRP=16847.1;SRR=0;TYPE=snp   GT:DP:RO:QR:AO:QA:GL    0/0:7760:7757:280123:3:111:0,-10,-10

My awk statement:

awk 'v {
        ro=gensub(/^.*;RO=([0-9]*).*$/, "\\1", "1");
        printf "%s %f\n", f, (ao/(ao + ro)); v=0
    }
    /^chr/ {ao=gensub(/^.*;AO=([0-9]*).*$/,"\\1", "1");
        v=1;
        f=$1 FS $2 FS $4 FS $5
        }' $VCFDIR/$currentfile.vcf >> $VCFDIR/totalvariants.txt #output all to single file

awk '{print $5}' $VCFDIR/totalvariants.txt > $VCFDIR/allelefreqs.txt

For Ref this is the entire Script I'm running using Karakfa's answer.

currentfile=finalOutput.bam

VCFDIR=/dir

awk '/^chr/{match($0,/;AO=([^;]*);.*;RO=([^;]*);/,v); 
            print $1,$2,$4,$5,v[1]/(v[1]+v[2])}' $VCFDIR/$currentfile.vcf >> $VCFDIR/totalvariants.txt #output all to single file

awk '{print $5}' $VCFDIR/totalvariants.txt > $VCFDIR/allelefreqs.txt

Upvotes: 1

Views: 100

Answers (3)

Ed Morton
Ed Morton

Reputation: 203209

$ cat tst.awk
BEGIN { FS="[ =;]+" }
/^chr/ {
    for (i=1;i<=NF;i++)
        a[$i] = $(i+1)
    print  a["AO"]/(a["AO"]+a["RO"])
}

$ awk -f tst.awk file
0.000386598

Upvotes: 1

karakfa
karakfa

Reputation: 67467

awk to the rescue!

awk '/^chr/{match($0,/;AO=([^;]*);.*;RO=([^;]*);/,v); 
            print $1,$2,$4,$5,v[1]/(v[1]+v[2])}'

expects AO to appear before RO

Using the test data dna

chr4    110541280   .   G   A   0   .   AB=0;ABP=0;AC=0;AF=0;AN=2;AO=3;CIGAR=1X;DP=7760;DPB=7760;DPRA=0;EPP=9.52472;EPPR=16838.4;GTI=0;LEN=1;MEANALT=1;MQM=60;MQMR=60;NS=1;NUMALT=1;ODDS=10694.5;PAIRED=0;PAIREDR=0;PAO=0;PQA=0;PQR=0;PRO=0;QA=111;QR=280123;RO=7757;RPP=9.52472;RPPR=16838.4;RUN=1;SAF=3;SAP=9.52472;SAR=0;SRF=7757;SRP=16847.1;SRR=0;TYPE=snp   GT:DP:RO:QR:AO:QA:GL    0/0:7760:7757:280123:3:111:0,-10,-10
chr17 17474586 . C G 197.853 . AB=0;ABP=0;AC=2;AF=1;AN=2;AO=8;CIGAR=1X;DP=8;DPB=8;DPRA=0;EPP=20.3821;EPPR=0;GTI..=0;LEN=1;MEANALT=1;MQM=60;MQMR=0;NS=1;NUMALT=1;ODDS=15.6955;PAIRED=0;PAIREDR=0;PA..O=0;PQA=0;PQR=0;PRO=0;QA=296;QR=0;RO=0;RPP=20.3821;RPPR=0;RUN=1;SAF=0;SAP=20.3821..;SAR=8;SRF=0;SRP=0;SRR=0;TYPE=snp GT:DP:RO:QR:AO:QA:GL 1/1:8:0:0:8:296:-10,-2.40824,0
chr1 58910197 . T G 52.2409 . AB=0;ABP=0;AC=2;AF=1;AN=2;AO=2;CIGAR=1X;DP=2;DPB=2;DPRA=0;EPP=7.35324;EPPR=0;GTI..=0;LEN=1;MEANALT=1;MQM=60;MQMR=0;NS=1;NUMALT=1;ODDS=7.37776;PAIRED=0;PAIREDR=0;PA..O=0;PQA=0;PQR=0;PRO=0;QA=74;QR=0;RO=0;RPP=7.35324;RPPR=0;RUN=1;SAF=2;SAP=7.35324;..SAR=0;SRF=0;SRP=0;SRR=0;TYPE=snp GT:DP:RO:QR:AO:QA:GL 1/1:2:0:0:2:74:-7.03,-0.60206,0

I get this output

chr4 110541280 G A 0.000386598
chr17 17474586 C G 1
chr1 58910197 T G 1

After checking your usage, you can eliminate second awk call as well

awk '/^chr/{match($0,/;AO=([^;]*);.*;RO=([^;]*);/,v);
        r=v[1]/(v[1]+v[2]);
        print $1,$2,$4,$5,r; 
        print r > "$VCFDIR/allelefreqs.txt"}'

Upvotes: 2

Jose Ricardo Bustos M.
Jose Ricardo Bustos M.

Reputation: 8164

you can try

awk '
  /^chr/ {
     nrecords = split($8, records, ";")       # "AB=0" "ABP=0" "AC=0" ...
     for(i=1; i<=nrecords; i++){
       split(records[i], record, "=");        # dict["AB"]="0" ...
       dict[record[1]] = record[2];           # store each record in column 8
     }
     frequency_allele = dict["AO"]/(dict["AO"]+dict["RO"]);
     print $1,$2,$4,$5,frequency_allele;
}' file

with input

chr4    110541280   .   G   A   0   .   AB=0;ABP=0;AC=0;AF=0;AN=2;AO=3;CIGAR=1X;DP=7760;DPB=7760;DPRA=0;EPP=9.52472;EPPR=16838.4;GTI=0;LEN=1;MEANALT=1;MQM=60;MQMR=60;NS=1;NUMALT=1;ODDS=10694.5;PAIRED=0;PAIREDR=0;PAO=0;PQA=0;PQR=0;PRO=0;QA=111;QR=280123;RO=7757;RPP=9.52472;RPPR=16838.4;RUN=1;SAF=3;SAP=9.52472;SAR=0;SRF=7757;SRP=16847.1;SRR=0;TYPE=snp   GT:DP:RO:QR:AO:QA:GL    0/0:7760:7757:280123:3:111:0,-10,-10
chr17 17474586 . C G 197.853 . AB=0;ABP=0;AC=2;AF=1;AN=2;AO=8;CIGAR=1X;DP=8;DPB=8;DPRA=0;EPP=20.3821;EPPR=0;GTI=0;LEN=1;MEANALT=1;MQM=60;MQMR=0;NS=1;NUMALT=1;ODDS=15.6955;PAIRED=0;PAIREDR=0;PAO=0;PQA=0;PQR=0;PRO=0;QA=296;QR=0;RO=0;RPP=20.3821;RPPR=0;RUN=1;SAF=0;SAP=20.3821;SAR=8;SRF=0;SRP=0;SRR=0;TYPE=snp GT:DP:RO:QR:AO:QA:GL 1/1:8:0:0:8:296:-10,-2.40824,0

you get

chr4 110541280 G A 0.000386598
chr17 17474586 C G 1

Upvotes: 1

Related Questions