User Tools

Site Tools


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