CODEWITHSUNDEEP
unixawkvcf-vcard
R_G
R_G

Reputation: 1

unix: counting samples with two variants or a homozygous variant per gene

I have a coding problem beyond my limited skills with unix power tools. I'm looking to count the number of sample with either: i) a homozygous variant in a gene (BB below); or ii) two variants in a gene (2x AB). For example, from:

Variant Gene    Sample1 Sample2 Sample3
   1    TP53    AA  BB  AB
   2    TP53    AB  AA  AB
   3    TP53    AB  AA  AA
   4    KRAS    AA  AB  AA
   5    KRAS    AB  AB  BB

I'm looking for:

Gene Two_variants Homozygous Either
TP53     2            1        3
KRAS     1            1        2

Any help would be much appreciated. Thanks.

R_G

Upvotes: 0

Views: 108

Answers (2)

perreal
perreal

Reputation: 97958

use warnings;
use strict;
my (@header, %data);
open(my $file, "<", "input") or die("$?");
while (<$file>) {
    @header = split, next if not @header;
    my @v = split;
    $data{$v[1]}->{$_}++ for (@v[2..$#v]);
}
close $file;
print "Gene Two_variants Homozygous Either\n";
for my $k (keys %data) {
    my ($var2, $homoz) = (int($data{$k}{AB}/2), $data{$k}{BB});
    my $sum = $var2 + $homoz;
    printf("%4s %8d %9d %8d\n", $k, $var2, $homoz, $sum) if $sum;
}

Upvotes: 0

gbrener
gbrener

Reputation: 5835

In GNU awk:

awk '/\<AB\>.+\<AB\>/ { arr[$2,"AB"] += 1 }
             /\<BB\>/ { arr[$2,"BB"] += 1 }
                  END { for ( elt in arr ) {
                          split ( elt, index_parts, SUBSEP )
                          genes[index_parts[1]] = 0
                        }
                        printf "%4s%13s%11s%7s\n", "Gene", "Two_variants", "Homozygous", "Either"
                        for ( gene in genes ) {
                          printf "%4s%6d%13d%9d\n", gene, arr[gene,"AB"], arr[gene,"BB"], arr[gene,"AB"] + arr[gene,"BB"]
                        }
                      }' input.txt

Upvotes: 1

Related Questions