-
Notifications
You must be signed in to change notification settings - Fork 0
/
gwa.gait2.solarius.RNA.R
71 lines (53 loc) · 2.44 KB
/
gwa.gait2.solarius.RNA.R
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
##' Calculate the association of SNPs that come out the first screening with gwa.gait2.MatrixEQTL.RNA in GAIT2
##'
##' This function applies to the output of the function "gwa.gait2.solariusRNA.getData"
##' @title Refine the results with solarius of RNAseq GWAs
##' @param out the output of function "gwa.gait2.solariusRNA.getData"
##' @param chr the chromosome
##' @return data.frame with the results of the gwas, with all gens
##' @author Angel Martinez-Perez, \email{[email protected]}
##' @examples
##' out <- gwa.gait2.solariusRNA.getData(chr=22, inDir='gwaRNAseq')
##' gwa.gait2.solarius.RNA(out, chr=22)
##'
gwa.gait2.solarius.RNA <- function(out , chr)
{
trait <- out[1,2]
stopifnot(chr %in% 1:22)
coresMax <- 60
stopifnot(class(out)=='data.frame')
df.RNA <- load.RNAseqGAIT2(gens=c('id',trait))
df <- plyr::join(df.Core,df.RNA, by = 'id', type = "left")
row.names(df) <- df$id
SNPs <- out$SNP
## load the genotype data from this SNPs
df.snps <- load.gwaGAIT2(chr=chr, snps = c('id',SNPs))
row.names(df.snps) <- df.snps$id
df.snps <- df.snps %>% select(-id) %>% as.matrix
df.map.all <- df.map.all2[SNPs, ]
df.map <- df.map.all %>% select(c('rsIDref','chromosome','position'))
df.map <- dplyr::rename(df.map, SNP=rsIDref)
## el solar da problemas con variables con nombre largo
cambiar <- nchar(SNPs) > 19
if(any(cambiar)){
originales <- SNPs[cambiar]
sinteticos <- paste0(paste0(sample(c(letters,LETTERS),6,replace=TRUE),collapse=''),1:sum(cambiar))
colnames(df.snps) <- plyr::mapvalues(colnames(df.snps), from=originales, to=sinteticos)
df.map$SNP <- plyr::mapvalues(df.map$SNP, from=originales, to=sinteticos)
}
## build the formula
formula.sol <- as.formula(paste(trait, '~ 1'))
cores <- min(coresMax,ncol(df.snps))
## call solarius
assoc.solarius <- solarAssoc(formula.sol, data= df ,snpcovdata=df.snps,snpmap=df.map, cores=cores,kinship=kin2)
## y finalmente deshacemos el cambio de nombre:
resultados <- assoc.solarius$snpf
if(any(cambiar))
{
resultados$SNP <- plyr::mapvalues(resultados$SNP, from=sinteticos, to=originales)
}
resultados <- merge(resultados,df.map.all[,-which(colnames(df.map.all) %in% c('chromosome','position'))],by.x='SNP',by.y='rsIDref')
## rename a column.name
resultados <- dplyr::rename(resultados, rsIDref= SNP)
return(resultados)
}