welch16
2019-09-10 db8fcaa6e55028efe111653924b301ccf7892813
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
 
library(magrittr)
library(tidyverse)
 
### we are going to 
 
metadata <- read_tsv(here::here("rawdata","metadata","emca4_snps.tsv")) %>% 
    rename(rs_id = SNP)
 
indr <- here::here("data","filtered")
outdr <- here::here("data","clean")
 
cases <- file.path(indr,"genotype_CGEMS_CASE.rds") %>% readRDS()
controls <- file.path(indr,"genotype_CGEMS_CONTROL.rds") %>% readRDS()
 
rename_variables <- function(data)
{
    set_names(data, c("rs_id","specimen","individual","allele1","allele2"))    
}
 
## remove specimen column because it is not necessary for our 
## purposes
cases %<>% rename_variables() %>% select(-specimen)
controls %<>% rename_variables() %>% select(-specimen)
 
cases %>% group_by(rs_id) %>% tally()
controls %>% group_by(rs_id) %>% tally()
 
cases %>% group_by(individual) %>% tally() %>% nrow()
controls %>% group_by(individual) %>% tally() %>% nrow()
 
## calculate allele frequency to determine the risk allele
## defined as the minor allele
 
## parse case and control information as a binary variables
cases %<>% mutate( status = 1)
controls %<>% mutate( status = 0)
cgems <- bind_rows(cases,controls)
 
## as expected, in this df, the allele 1 NA's are the same individuals 
## than in allele2
 
cgems %<>% 
    nest(-rs_id, .key = "ind_data") %>% 
    mutate(
        alleles = map(ind_data, ~ unique(c(.$allele1,.$allele2))) %>% map(sort),
        allele1 = map_chr(alleles,1),
        allele2 = map_chr(alleles,2),
        n = map_int(ind_data,nrow),
        n_NA = map(ind_data, filter,is.na(allele1)) %>% map_int(nrow),
        n_AA = map2(ind_data, alleles, ~ filter(.x, allele1 == .y[1] & allele2 == .y[2])) %>% 
            map_int(nrow),
        n_aa = map2(ind_data, alleles, ~ filter(.x, allele1 == .y[2] & allele2 == .y[2])) %>% 
            map_int(nrow),
        n_Aa = n - n_NA - n_AA - n_aa,
        allele1_freq = ( 2 * n_AA + n_Aa) / (2 * n),
        risk_allele = if_else(allele1_freq > 0.5, allele2, allele1)) %>%   
    select(-alleles,-allele1,-allele2) 
 
cgems_summary <- cgems %>% 
    select(-ind_data) %>% 
    mutate( raf = pmin(allele1_freq, 1 - allele1_freq)) %>% 
    select(rs_id,risk_allele,raf,n,n_NA)
 
cgems_summary %<>% left_join(metadata,by = "rs_id") %>% arrange(-n_NA) %>% 
    mutate(
        alleles2 = str_split(alleles,"\\/") %>% map(unlist),
        is_neg_strand = ! map2_lgl(risk_allele,alleles2, ~ .x %in% .y),
        is_neg_strand = if_else(is_neg_strand,"yes","no")) %>% 
    select(-alleles2)
    
cgems_summary
 
cgems %<>% select(rs_id,ind_data,risk_allele) %>% unnest()
 
cgems
 
cgems %<>% 
    mutate(
        contains_risk_allele = allele1 == risk_allele | allele2 == risk_allele,
        contains_risk_allele = case_when(
            is.na(contains_risk_allele) ~ NA_real_,
            contains_risk_allele ~ 1,
            ! contains_risk_allele ~ 0
            ))
 
cgems_mat <- cgems %>%
    select(rs_id,individual,status,contains_risk_allele) %>% 
    spread(rs_id,contains_risk_allele)
    
## 
cgems_snps <- cgems_mat %>% as.data.frame() %>% 
    column_to_rownames("individual") %>% 
    select(-status) %>% as.matrix()
 
cgems_summary %>% mutate( p_NA = n_NA / n) %>% head(2)
 
cgems_summary %>% filter(is_neg_strand == "yes") %>% as.data.frame %>% 
    select(rs_id,risk_allele,alleles)
 
cgems_summary %>% write_tsv(file.path(indr,"2019_07_22_snp_summary.tsv"))
 
snps_to_keep <- cgems_summary %>% filter(!(n_NA  > 400 | raf  <= .05) ) %>% 
    pull(rs_id)
    
## remove snps with low risk allele frequency or higher number of missing values
cgems_mat %<>% select(individual,status,one_of(snps_to_keep))
cgems_snps <- cgems_snps[,snps_to_keep]
 
nsnps <- ncol(cgems_snps)
 
pmiss_sample <- apply( cgems_snps , 1 , function(x)sum(is.na(x))) / nsnps
 
### only keep samples with less than 10% of missing genotypes
## this means we will remove 35 samples
samples_to_keep <- names(pmiss_sample)[pmiss_sample < .1]
 
cgems_mat %<>% filter(individual %in% samples_to_keep)
cgems_snps <- cgems_snps[samples_to_keep,]
 
cgems_mat %>% count(status)
 
#### still similar
 
# # A tibble: 2 x 2
#   status     n
#    <dbl> <int>
# 1      0  1189
# 2      1  1210
 
cgems_mat_by_snps <- cgems_mat %>% gather(snp,allele,-status,-individual)
cgems_mat_by_snps %<>% na.omit() %>% nest(-snp)
 
cgems_mat_by_snps %<>% 
    mutate(
        model_snp = map(data, ~ glm(status ~ 0 + allele, data = . ,family = binomial(link = "logit"))))
        
cgems_mat_by_snps %<>% 
    mutate(
        snp_coeff = map_dbl(model_snp, ~ coef(.)[1]),
        odds_ratio = exp(snp_coeff),
        risk_allele = if_else(odds_ratio  > 1,"yes","no"))
        
cgems_mat_by_snps %>% count(risk_allele)
 
cgems_mat_by_snps %>% select(-data,-model_snp) %>% saveRDS(file.path(outdr,"2019_09_10_snp_odds_ratio.rds"))