User Tools

Site Tools


genetica:preproc_models

Differences

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

Link to this comparison view

genetica:preproc_models [2013/05/16 12:18]
osotolongo [Calculando los odds ratio y demas]
genetica:preproc_models [2020/08/04 10:58]
Line 1: Line 1:
-====== Metanalisis de modelos de plink ====== 
  
-===== Corriendo los modelos ====== 
- 
-Para sacar el analisis de modelos que esta implementado en plink se ejecuta algo como 
- 
-<code> 
-$ plink plink --bfile archivo --model --allow-no-sex --out archivo 
-</code> 
- 
-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, 
- 
-<code bash> 
-$ for x in `ls ../Variomics/*.bed`; do plink --bfile ${x%.bed} --model --allow-no-sex --out $(basename ${x%.bed}); done 
-</code> 
- 
-Esto deja una serie de archivos de resultado, 
- 
-<code> 
-[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 
-</code> 
- 
-===== 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, 
- 
-<code perl models.pl> 
-#!/usr/bin/perl 
- 
-# 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 Text::NSP::Measures::2D::Fisher::twotailed; 
- 
-my $plist = shift; 
-my $model = shift; 
-my $bimfile = shift; 
-die "usage: models.pl model_file model_tla bim_file\n" unless ($plist && $model && $bimfile);  
- 
-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 %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 %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/($affa1+$affa2); 
- $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)); 
- } 
-} 
- 
-my $head = "CHR\tSNP\tBP\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$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"; 
-} 
- 
-close ODF; 
-</code> 
- 
-La orden, 
- 
-<code> 
-$ models.pl admurcia_model.model REC ../Variomics/ADMURimpQC2.bim 
-</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//. 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, 
-<code> 
-$ for x in `ls *.model | grep -v recalc`; do models.pl $x REC ../Variomics/${x%.model}.bim; done 
-</code> 
-y 
-<code> 
-$ for x in `ls *.model | grep -v recalc`; do models.pl $x DOM ../Variomics/${x%.model}.bim; done 
-</code> 
- 
-dado que me interesan solo los modelos dominate y recesivo. 
-===== Meta analisis con plink ===== 
-Para terminar basta con correr plink sobre todas las bases procesadas, 
-<code> 
-$ a=$(ls *REC*); plink --meta-analysis $a + study --out recesivo 
-</code> 
-y 
-<code> 
-$ a=$(ls *DOM*); plink --meta-analysis $a + study --out dominante 
-</code> 
genetica/preproc_models.txt · Last modified: 2020/08/04 10:58 (external edit)