Differences

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

--- genetica:preproc_models [2013/05/15 15:16] – [Meta analisis con plink] osotolongo
+++ genetica:preproc_models [2013/05/21 10:45] – [Alternativa a Fisher: Barnard Test] osotolongo
@@ Line 39: / Line 39: @@
 use strict; use warnings;
 use File::Slurp qw(read_file);
-use Data::Dump qw(dump);
 use Text::NSP::Measures::2D::Fisher::twotailed;
@@ Line 63: / Line 62: @@
 		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 67: @@
 		$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 93: / Line 91: @@
 </code>
-extrae el modelo recesivo, calcula los valores deseados y los añade a los ya existentes en el archivo //admurcia_model.REC.recalc.model//.
+extrae el modelo recesivo, calcula los valores deseados y los añade a los ya existentes en el archivo //admurcia_model.REC.recalc.model//. Además usa el archivo **.bim** original para asignar a cada marcador su posicion fisica ya que //plink// la pide para el meta-analisis
 Para correr todas las DBs juntas hago,
@@ Line 108: / Line 106: @@
 Para terminar basta con correr plink sobre todas las bases procesadas,
 <code>
-$ a=$(ls *REC*); plink --meta-analysis $a --out recesivo
+$ a=$(ls *REC*); plink --meta-analysis $a + study --out recesivo
 </code>
 y
 <code>
-$ a=$(ls *DOM*); plink --meta-analysis $a --out dominante
+$ 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.