forked from elijahedmondson/NASA-GRSD-HS
-
Notifications
You must be signed in to change notification settings - Fork 0
/
HS.cox.RMA.chrom.R
112 lines (74 loc) · 4.93 KB
/
HS.cox.RMA.chrom.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
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
HS.cox.RMA.chrom = function(perms, chr, pheno, pheno.col, days.col, probs, K, addcovar,
markers, snp.file, outdir = "~/Desktop/",
tx = "", sanger.dir = "~/Desktop/R/QTL/WD/HS.sanger.files/")
{
begin <- Sys.time()
file.prefix = paste0("Bootstrap.", tx, ".", pheno.col, "(Chr", chr, ")")
plot.title = paste("Bootstrap ", tx, " ", pheno.col, " ", "(Chr ", chr, "):", sep = "")
print(paste(plot.title, Sys.time()))
load(file = paste0(sanger.dir, "probs", chr, ".Rdata"))
samples = intersect(rownames(pheno), rownames(probs))
samples = intersect(samples, rownames(addcovar))
samples = intersect(samples, rownames(K[[1]]))
stopifnot(length(samples) > 0)
pheno = pheno[samples,,drop = FALSE]
addcovar = addcovar[samples,,drop = FALSE]
probs = probs[samples,,,drop = FALSE]
#samples2 = markers$SNP_ID[which(markers$Chr == chr & markers$Mb_NCBI38 > (peakMB - 4000000) & markers$Mb_NCBI38 < (peakMB + 4000000))]
samples2 = markers$SNP_ID[which(markers$Chr == chr)]
probs = probs[,,samples2, drop = FALSE]
K = K[[chr]][samples,samples,drop = FALSE]
setwd(outdir)
##
permutations = matrix(1, nrow = perms, ncol = 5, dimnames = list(1:perms, c("LOD", "min", "max", "average", "#Markers")))
rm(samples, samples2)
for(p in 1:perms) {
LODtime = Sys.time()
print(p)
#repeat {
tryCatch({
phenoperm = data.frame(pheno[sample(nrow(pheno), replace = TRUE), ],
check.names = FALSE, check.rows = FALSE)
sanger.dir = sanger.dir
### FROM DG ###
samples = sub("\\.[0-9]$", "", rownames(phenoperm))
probsperm = probs[samples,,]
rownames(probsperm) = make.unique(rownames(probsperm))
Kperm = K
Kperm = K[samples,samples,drop = FALSE]
rownames(Kperm) = make.unique(rownames(Kperm))
### Move the model into the loop REGRESSION MODEL ###
Kperm = Kperm[samples, samples]
chrs = chr
data = vector("list", length(chrs))
names(data) = chrs
rng = which(markers[,2] == chrs)
data = list(probsperm = probsperm, Kperm = Kperm,
markers = markers[rng,])
result = vector("list", length(data))
names(result) = names(data)
rm(probsperm, Kperm)
### WORK HORSE ###
### RETURN ANY POS FOR ENRIRE CHROMOSOME ###
load("/Users/elijah/Desktop/R/QTL/WD/hs.colors.Rdata")
result = GRSD.coxph4perms(data, chr = chr, pheno = phenoperm, pheno.col, days.col, addcovar, tx, sanger.dir)
top = max(-log10(result$pv))
MAX.LOD = result$POS[which(-log10(result$pv) == top)]
MegaBase = (min(MAX.LOD) + max(MAX.LOD))/2000000
print(paste(MegaBase, "Mb: LOD", round(top, digits = 2)))
#if (MAX.LOD > (peakMB - (window/2)) & MAX.LOD < (peakMB + (window/2))) break
}, error=function(e){})
#}
#print(paste0("Accepted locus: ", MegaBase, " Mb"))
# Save the locus.
permutations[p,] = c(top, min(MAX.LOD), max(MAX.LOD), ((min(MAX.LOD) + max(MAX.LOD))/2), length(MAX.LOD))
print(paste(round(difftime(Sys.time(), LODtime, units = 'mins'), digits = 2),
"minutes..."))
}
print(paste(round(difftime(Sys.time(), begin, units = 'hours'), digits = 2),
"hours elapsed during analysis"))
quant = quantile(permutations[,4], c(0.025,0.975))
print(paste("95% Confidence Interval for QTL:", min(quant), "-", max(quant)))
print(paste("Interval =", (round((((max(quant) - min(quant))/1000000)), digits = 2)), "Mb"))
return(permutations)
} #HS.cox.bootstrap