genetica:zscore_code
Al respecto ver http://www.perlmonks.org/?node_id=1015649
use strict; use warnings; use File::Slurp qw(read_file); use Math::CDF qw(qnorm pnorm); use List::MoreUtils qw(uniq); my $ofile = "meta1.txt"; my @ifiles = @ARGV; my %ipairs; my @lpairs; foreach my $ifile (@ifiles){ (my $fk) = $ifile =~ /^(.*)\_sets.*/; my %ldata = reverse map {/^(.*(rs\d{1,20}\s+rs\d{1,20}).*)$/} grep {/.*rs\d{1,20}\s+rs\d{1,20}.*/} read_file $ifile; foreach my $dline (sort keys %ldata){ push @lpairs, $dline; ($ipairs{$fk}{$dline}{'head'}, $ipairs{$fk}{$dline}{'effect'}, $ipairs{$fk}{$dline}{'pvalue'}) = $ldata{$dline} =~ /^(.*)\s+(\d\.\d+)\s+\d\.\d+\s+(\d\.\d+)$/; } } @lpairs = uniq @lpairs; open OF, ">$ofile"; my $head = "CHR1 CHR2 SNP1 SNP2 P N"; print OF "$head\n"; foreach my $pair (@lpairs) { my $n = 0; my $z = 0; my $hl; my $pvalue = 0; my $fk; foreach $fk (%ipairs) { my $pvt = 0; if( exists $ipairs{$fk} && exists $ipairs{$fk}{$pair} && exists $ipairs{$fk}{$pair}{'pvalue'}){ #if($ipairs{$fk}{$pair}{'pvalue'}){ $pvt = $ipairs{$fk}{$pair}{'pvalue'}; if($pvt){ unless($hl){ $hl = $ipairs{$fk}{$pair}{'head'}; } $n++; $z+= qnorm($ipairs{$fk}{$pair}{'pvalue'}); } } } if($n>2){ $z = $z/sqrt($n); $pvalue = pnorm($z); } if ($pvalue) { printf "$pair -> %.4f\n", $pvalue; printf OF "$hl %.4f $n\n", $pvalue; } } close OF;
genetica/zscore_code.txt · Last modified: 2020/08/04 10:58 by 127.0.0.1