User Tools

Site Tools


genetica:preproc_models

This is an old revision of the document!


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 Data::Dump qw(dump);
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){
	#Ejemplo
	#(@{$genmap_previous{$marker}} {qw/chr position allele1 freq1 freq2 allele2 chi2 pvalue oddratio/}) = $data_line{$marker} =~ /^\s*(\d+)\s+rs\d+\s+(\d+)\s+([A,T,C,G]+)\s+(\d+\.*\d*e*-*\d*)\s+(\d+\.*\d*e*-*\d*)\s+([A,T,C,G]+)\s+(\d+\.*\d*e*-*\d*)\s+(\d+\.*\d*e*-*\d*)\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+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;
		#print "$affa1, $t1, $t2, $t3\n";
		$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";
 
#{($genmap{$a}->{chr} <=> $genmap{$b}->{chr}) or ($genmap{$a}->{fposition} <=> $genmap{$b}->{fposition})}
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.

genetica/preproc_models.1368618384.txt.gz · Last modified: 2020/08/04 10:48 (external edit)