This is an old revision of the document!
Para sacar el analisis de modelos que esta implementado en plink se ejecuta algo como
$ plink plink --bfile archivo --model --allow-no-sex --out archivo
Esto deja el resultado en un archivo llamado archivo.model. Como aqui el objetivo final es un meta analisi vamos a ahcerlo para todas las DB del Variomics,
$ for x in `ls ../Variomics/*.bed`; do plink --bfile ${x%.bed} --model --allow-no-sex --out $(basename ${x%.bed}); done
Esto deja una serie de archivos de resultado,
[osotolongo@detritus models]$ ls -lah *.model -rw-rw-r-- 1 osotolongo osotolongo 479M May 15 13:18 ADMURimpQC2.model -rw-rw-r-- 1 osotolongo osotolongo 831M May 15 13:20 ADNIimpQC2.model -rw-rw-r-- 1 osotolongo osotolongo 692M May 15 13:22 GenADA_impQC2.model -rw-rw-r-- 1 osotolongo osotolongo 805M May 15 13:25 NIA_AD_impQC2.model -rw-rw-r-- 1 osotolongo osotolongo 597M May 15 13:26 TGEN_impQC2.model
Quiero extraer un determinado modelo (e.g.: recesivo → REC) de estos archivos y calcular los odds ratio, standard error y p-value. Para calcular las p-value voy a usar el Fisher's Exact Test dado que hay alelos de muy baja frecuencia a los que no puedo hacer un $\chi^2$.
Hice un scriptcillo para ayudarme,
#!/usr/bin/perl # Copyright 2013 O. Sotolongo <osotolongo@fundacioace.com> use strict; use warnings; use File::Slurp qw(read_file); use Text::NSP::Measures::2D::Fisher::twotailed; my $plist = shift; my $model = shift; my $ofile = $plist; $ofile =~ s/(.*)\.(.*)/$1\.$model\.recalc\.$2/; my %input_data = reverse map {/^(.*\s+(rs\d+)\s+.*)$/} grep {/^(.*\s+$model\s+.*)$/} read_file $plist; 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*$/; if(exists($cells{$marker}{affa1}) && exists($cells{$marker}{affa2}) && exists($cells{$marker}{unaffa1}) && exists($cells{$marker}{unaffa2})){ my $affa1 = $cells{$marker}{affa1}; my $affa2 = $cells{$marker}{affa2}; my $unaffa1 = $cells{$marker}{unaffa1}; my $unaffa2 = $cells{$marker}{unaffa2}; my $t1 = $affa1 + $affa2; my $t2 = $affa1 + $unaffa1; my $t3 = $t1 + $unaffa1 + $unaffa2; $cells{$marker}{pvalue} = calculateStatistic(n11=>$affa1, n1p=>$t1, np1=>$t2, npp=>$t3); $affa1 = 0.1 unless ($affa1); $affa2 = 0.1 unless ($affa2); $unaffa1 = 0.1 unless ($unaffa1); $unaffa2 = 0.1 unless ($unaffa2); $cells{$marker}{afreq} = $affa1/$affa2; $cells{$marker}{ufreq} = $unaffa1/$unaffa2; $cells{$marker}{oddsratio} = ($affa1*$unaffa2)/($affa2*$unaffa1); $cells{$marker}{stderr} = sqrt((1/$affa1)+(1/$affa2)+(1/$unaffa1)+(1/$unaffa2)); } } my $head = "CHR\tSNP\tA1\tA2\tF_A\tF_U\tCHISQ\tDF\tP0\tOR\tSE\tP"; open ODF, ">$ofile" or die "Could not open output file\n"; print ODF "$head\n"; 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$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"; } close ODF;
La orden,
$ models.pl admurcia_model.model REC
extrae el modelo recesivo, calcula los valores deseados y los añade a los ya existentes en el archivo admurcia_model.REC.recalc.model.
Para correr todas las DBs juntas hago,
$ for x in *.model; do models.pl $x REC; done
y
$ for x in *.model; do models.pl $x DOM; done
dado que me interesan solo los modelos dominate y recesivo.