Differences

This shows you the differences between two versions of the page.

--- genetica:preproc_models [2013/05/15 15:40] – [Meta analisis con plink] osotolongo
+++ genetica:preproc_models [2013/05/23 07:17] – [Calculando los odds ratio y demas] osotolongo
@@ Line 37: / Line 37: @@
 # Copyright 2013 O. Sotolongo <osotolongo@fundacioace.com>
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
 use strict; use warnings;
 use File::Slurp qw(read_file);
-use Data::Dump qw(dump);
 use Text::NSP::Measures::2D::Fisher::twotailed;
@@ Line 48: / Line 57: @@
 my $ofile = $plist;
-$ofile =~ s/(.*)\.(.*)/$1\.$model\.recalc\.$2/;
+$ofile =~ s/(.*)\.(.*)/$1\.$model\.recalc2\.$2/;
+# cambiar a clustering con (?:PATTERN)
-my %input_data = reverse map {/^(.*\s+(rs\d+)\s+.*)$/} grep {/^(.*\s+$model\s+.*)$/} read_file $plist;
+my %input_data = reverse map {/^(\s*\d+\s+(\S+)\s+[A,T,C,G]+\s+[A,T,C,G]+\s+$model\s+.*)$/} grep {/^(.*\s+$model\s+.*)$/} read_file $plist;
-my %bim_data = map {/^\d+\s+(rs\d+)\s+\d+\.*\d*\s+(\d+)\s+[A,T,C,G]+\s+[A,T,C,G]+\s*$/} grep {/.*rs\d+.*/} read_file $bimfile;
+my %bim_data = map {/^\d+\s+(.*)\s+\d+\.*\d*\s+(\d+)\s+[A,T,C,G]+\s+[A,T,C,G]+\s*$/} read_file $bimfile;
 my %cells;
 foreach my $marker (sort keys %input_data){
-	(@{$cells{$marker}} {qw/chr allele1 allele2 affa1 affa2 unaffa1 unaffa2 chi2 df pval0/}) = $input_data{$marker} =~ /^\s*(\d+)\s+rs\d+\s+([A,T,C,G])+\s+([A,T,C,G])+\s+$model\s+(\d+)\/(\d+)\s+(\d+)\/(\d+)\s+(\d+\.*\d*e*-*\d*|NA)\s+(\d+\.*\d*e*-*\d*|NA)\s+(\d+\.*\d*e*-*\d*|NA)\s*$/;
+	(@{$cells{$marker}} {qw/chr allele1 allele2 affa1 affa2 unaffa1 unaffa2 chi2 df pval0/}) = $input_data{$marker} =~ /^\s*(\d+)\s+$marker\s+([A,T,C,G])+\s+([A,T,C,G])+\s+$model\s+(\d+)\/(\d+)\s+(\d+)\/(\d+)\s+(\d+\.*\d*e*-*\d*|NA)\s+(\d+\.*\d*e*-*\d*|NA)\s+(\d+\.*\d*e*-*\d*|NA)\s*$/;
 	if(exists($cells{$marker}{affa1}) && exists($cells{$marker}{affa2}) && exists($cells{$marker}{unaffa1}) && exists($cells{$marker}{unaffa2})){
 		my $affa1 = $cells{$marker}{affa1};
@@ Line 63: / Line 72: @@
 		my $t2 = $affa1 + $unaffa1;
 		my $t3 = $t1 + $unaffa1 + $unaffa2;
-		#print "$affa1, $t1, $t2, $t3\n";
 		$cells{$marker}{pvalue} = calculateStatistic(n11=>$affa1, n1p=>$t1, np1=>$t2, npp=>$t3);
 		$affa1 = 0.1 unless ($affa1);
@@ Line 69: / Line 77: @@
 		$unaffa1 = 0.1 unless ($unaffa1);
 		$unaffa2 = 0.1 unless ($unaffa2);
-		$cells{$marker}{afreq} = $affa1/$affa2;
+		$cells{$marker}{afreq} = $affa1/($affa1+$affa2);
-		$cells{$marker}{ufreq} = $unaffa1/$unaffa2;
+		$cells{$marker}{ufreq} = $unaffa1/($unaffa1+$unaffa2);
 		$cells{$marker}{oddsratio} = ($affa1*$unaffa2)/($affa2*$unaffa1);
 		$cells{$marker}{stderr} = sqrt((1/$affa1)+(1/$affa2)+(1/$unaffa1)+(1/$unaffa2));
@@ Line 81: / Line 89: @@
 foreach my $marker (sort {($cells{$a}->{chr} <=> $cells{$b}->{chr}) or ($a cmp $b)} keys %input_data){
-	print ODF "$cells{$marker}{chr}\t$marker\t$bim_data{$marker}\t$cells{$marker}{allele1}\t$cells{$marker}{allele2}\t$cells{$marker}{afreq}\t$cells{$marker}{ufreq}\t$cells{$marker}{chi2}\t$cells{$marker}{df}\t$cells{$marker}{pval0}\t$cells{$marker}{oddsratio}\t$cells{$marker}{stderr}\t$cells{$marker}{pvalue}\n";
+	if(exists($bim_data{$marker})){
+		print ODF "$cells{$marker}{chr}\t$marker\t$bim_data{$marker}\t$cells{$marker}{allele1}\t$cells{$marker}{allele2}\t$cells{$marker}{afreq}\t$cells{$marker}{ufreq}\t$cells{$marker}{chi2}\t$cells{$marker}{df}\t$cells{$marker}{pval0}\t$cells{$marker}{oddsratio}\t$cells{$marker}{stderr}\t$cells{$marker}{pvalue}\n";
+	}
 }
@@ Line 114: / Line 124: @@
 $ a=$(ls *DOM*); plink --meta-analysis $a + study --out dominante
 </code>
+Los resultados salen ordenados por cromosoma y SNP. Para ordenar por //p-value// puede hacerse,
+<code>
+$ sort -g -k7,7 dominante.meta > dominante.meta.sorted
+</code>
+para cortar por un valor especifico de //p-value//,
+<code>
+$ awk {'if($7<5e-6 || $7=="P") print'} dominante.meta.sorted > dominante.meta.sorted.cut
+</code>
+y si estamos interesados en algun SNP especifico,
+<code>
+$ echo "DB `head -n 1 ADNIimpQC2.REC.recalc.model`" > 20p.rec.cool.txt
+$ for x in *REC*; do echo "${x%.REC.recalc.model} `grep -w "rs714235" $x`" >> 20p.rec.cool.txt; done
+</code>
+====== Alternativa a Fisher: Barnard Test ======
+Como al blandengue no le gusta el test de Fisher, he buscado como hacer que el script me calcule los test de Barnard. No he encontrado ningun modulo en Perl pero como todo ya esta inventado en esta vida hay un paquete en R que lo calcula. Asi que lo que he hecho es usar el modulo **Statistics::R** para calcular el test en R desde Perl. ;-)
+Esto es un poco enredado pero se resume en sustituir la llamada al test de Fisher con,
+<code perl>
+my $shit = Statistics::R->new();
+$shit->run(q`library(Barnard)`);
+$shit->run(qq`barnardw.test($affa1, $affa2, $unaffa1, $unaffa2, verbose=FALSE)->bt`);
+$shit->run(q`bt["p.value"][[1]][[2]] -> x`);
+$cells{$marker}{pvalue} = shit->get('x');
+</code>
+**OJO:** Esto demora muchisimo. Despues de 5 días de calculo sobre una base de datos y modelo dominate tuve que para la ejecucion. Lo lanzare cuando tenga tiempo y la base de datos a usar este definida.