Trabajaremos con el paquete SNPassoc

library(SNPassoc)
## Registered S3 method overwritten by 'SNPassoc':
##   method            from       
##   summary.haplo.glm haplo.stats

Carganos nuestro dataframe

library(readxl)
dbpol <- read_excel("C:/Users/fidel/OneDrive - UNIVERSIDAD AUTONOMA DE SINALOA/COLABS/DRa Guadron/dbpol.xlsx")

#lo visualizamos
#en la columnas casecontrol 1= caso y 0= control
#las columnas restantes son los polimorfismos
head(dbpol)
## # A tibble: 6 × 3
##   casecontrol rs429358 rs7412
##         <dbl> <chr>    <chr> 
## 1           1 TC       CC    
## 2           1 TT       CC    
## 3           1 TC       CC    
## 4           1 TC       TC    
## 5           1 TC       CC    
## 6           1 TT       CC

Se preparan los datos, para identificación de los polimorfismos, vamos a utilizar la función setupSNP del paquete SNPassoc

dbpol.s <- setupSNP(data=dbpol, colSNPs=2:ncol(dbpol), sep="")

Analizamos cada uno de los polimorfismos:

Primero el polimorfismo: rs429358

summary(dbpol.s$rs429358)
## Genotypes: 
##     frequency percentage
## T/T        74   50.68493
## C/T        57   39.04110
## C/C        15   10.27397
## 
## Alleles: 
##   frequency percentage
## T       205   70.20548
## C        87   29.79452
## 
## HWE (p value): 0.4304851

Observamos que se tiene frecuencia para los diferentes genotipos y para los respectivos alelos. Así como el equilibrio HW del polimorfismo rs429358.

Ahora analizamos para el polimorfismo: rs7412

summary(dbpol.s$rs7412)
## Genotypes: 
##     frequency percentage
## C/C       132  90.410959
## T/C        10   6.849315
## T/T         4   2.739726
## 
## Alleles: 
##   frequency percentage
## C       274  93.835616
## T        18   6.164384
## 
## HWE (p value): 0.000586435

Observamos frecuencia para los diferentes genotipos y para los respectivos alelos. Así como el equilibrio HW del polimorfismo rs7412.

Esto también lo podemos ver de manera visual:
tanto para el polimorfismo rs429358:

plot(dbpol.s$rs429358)

y para el polimorfismo rs7412.

plot(dbpol.s$rs7412)

también podemos hacerlo con un grafico de pastel:

plot(dbpol.s$rs429358, type=pie)

plot(dbpol.s$rs7412, type=pie)

Vamos a describir los alelos y hacer en una tabla el equilibrio HW

summary(dbpol.s)
##          alleles major.allele.freq HWE      missing (%)
## rs429358 T/C     70.2              0.430485 0          
## rs7412   C/T     93.8              0.000586 0
hwe <- tableHWE(dbpol.s)
head(hwe)
##          HWE (p value)
## rs429358   0.430485128
## rs7412     0.000586435

Ahora lo haremos considerando casos y controles

hwe2 <- tableHWE(dbpol.s, casecontrol)

hwe2
##          all.groups     X0     X1
## rs429358     0.4305 1.0000 1.0000
## rs7412       0.0006 0.0086 0.0291

recordemos que 0= control y 1= casos

Valos a calcular ahora HWE en la muestra complena pero sin controles

snpNHWE <- hwe2[,1]>0.05 & hwe2[,2]<0.05
rownames(hwe2)[snpNHWE]
## character(0)
hwe2[snpNHWE,]
##      all groups 0 1

Si el se rompe el equilibrio entonces lo correcto es eliminar el polimorfismo, aqui hay un codigo para filtrar esto y quedarnos solo con el polimorfismo que no rompe el equilibrio HW

## ----remove_SNPs--------------------------------------------------------------
snps.ok <- rownames(hwe2)[hwe2[,2]>=0.001]
pos <- which(colnames(dbpol.s)%in%snps.ok, useNames = FALSE)
dbpol.s <- setupSNP(dbpol.s, pos, sep="")

Ahora que ya tenemos listo nuestra base de datos y hemos determinado HWE en nuestros polimorfismos, vamos a hacer el analisis de asociación, generalmente se usa chi cuadrado, para determinar entre casos y controles. Lo vamos a hacer primero para el polimorfismo rs429358

## ----association--------------------------------------------------------------
association(casecontrol ~ rs429358, data = dbpol.s)
## 
## SNP: rs429358  adjusted by: 
##               0     %   1    %    OR lower upper   p-value   AIC
## Codominant                                                      
## T/T          32  86.5  42 38.5  1.00             1.171e-06 141.1
## C/T           5  13.5  52 47.7  7.92  2.84 22.12                
## C/C           0   0.0  15 13.8        0.00                      
## Dominant                                                        
## T/T          32  86.5  42 38.5  1.00             1.385e-07 141.5
## C/T-C/C       5  13.5  67 61.5 10.21  3.69 28.27                
## Recessive                                                       
## T/T-C/T      37 100.0  94 86.2  1.00             1.238e-02 160.0
## C/C           0   0.0  15 13.8        0.00                      
## Overdominant                                                    
## T/T-C/C      32  86.5  57 52.3  1.00             9.957e-05 154.1
## C/T           5  13.5  52 47.7  5.84  2.12 16.11                
## log-Additive                                                    
## 0,1,2        37  25.3 109 74.7  8.57  3.21 22.90 1.171e-06 139.4

Ahora para el polimorfismo: rs7412

association(casecontrol ~ rs7412, data = dbpol.s)
## 
## SNP: rs7412  adjusted by: 
##               0    %   1    %   OR lower upper p-value   AIC
## Codominant                                                  
## C/C          33 89.2  99 90.8 1.00              0.5269 170.0
## T/C           2  5.4   8  7.3 1.33  0.27  6.60              
## T/T           2  5.4   2  1.8 0.33  0.05  2.46              
## Dominant                                                    
## C/C          33 89.2  99 90.8 1.00              0.7728 169.2
## T/C-T/T       4 10.8  10  9.2 0.83  0.24  2.84              
## Recessive                                                   
## C/C-T/C      35 94.6 107 98.2 1.00              0.2833 168.1
## T/T           2  5.4   2  1.8 0.33  0.04  2.41              
## Overdominant                                                
## C/C-T/T      35 94.6 101 92.7 1.00              0.6802 169.1
## T/C           2  5.4   8  7.3 1.39  0.28  6.84              
## log-Additive                                                
## 0,1,2        37 25.3 109 74.7 0.75  0.32  1.75  0.5108 168.9

#AHORA SIGUE AJUSTAR EL MODELO PARA LAS DIFERENTES VARIABLES :)