This is an old revision of the document!
Table of Contents
Como preparar un metanalisi de los modelos de plink
Corriendo los modelos
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
Calculando los odds ratio y demas
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,
- models.pl
#!/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.