-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathmakeSd.R
126 lines (106 loc) · 4.23 KB
/
makeSd.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
113
114
115
116
117
118
119
120
121
122
123
124
#' @aliases makeSd makeD
#' @rdname makeD
#' @export
makeSd <- function(pedigree, heterogametic,
DosageComp = c(NULL, "ngdc", "hori", "hedo", "hoha", "hopi"),
parallel = FALSE, ncores = getOption("mc.cores", 2L),
invertSd = TRUE, returnS = FALSE, det = TRUE, verbose = TRUE){
if(length(unique(pedigree[,4])) > 2) stop("Error: more than 2 sexes specified")
dc.model <- match.arg(DosageComp)
if(is.null(dc.model)){
warning("Assuming 'ngdc' dosage compensation model")
dc.model <- "ngdc"
}
if(dc.model == "hopi" | dc.model == "hori"){
warning("Assume sex chromosomal dominance allelic interactions do not occur under 'hopi' or 'hori'\n")
return(NULL)
}
Sout <- makeS(pedigree, heterogametic = heterogametic,
DosageComp = dc.model, returnS = TRUE)
# makeA() returns `dsCMatrix`, but S is `dgCMatrix` from above
## makeD()-like code below expects symmetric matrix ('dsCMatrix')
S <- forceSymmetric(Sout$S)
nPed <- numPed(pedigree[, 1:3])
damsex <- pedigree[unique(nPed[, 2])[-1], 4]
if(any(damsex == heterogametic)){
pedname <- names(pedigree)
pedigree <- pedigree[, c(1,3,2,4)]
names(pedigree) <- pedname
nPed <- numPed(pedigree[, 1:3])
}
sex <- rep(-998, dim(pedigree)[1L])
sex[homs <- which(pedigree[,4] != heterogametic)] <- 1
sex[hets <- which(pedigree[,4] == heterogametic)] <- 0
nhom <- sum(sex) # Number of individuals with homogametic sex chromosomes
N <- dim(nPed)[1L]
#FIXME turned off next check so can test parallel=TRUE on small pedigrees
# if(parallel){
# if(length(S@x)/ncores < 10){
# warning("pedigree too small - 'parallel' set to FALSE instead")
# parallel <- FALSE
# }
# }
if(!parallel){
if(verbose) cat("starting to make Sd...")
Cout <- .C("sdij", PACKAGE = "nadiv",
as.integer(nPed[, 2] - 1), # [[1]] dam ID/No.
as.integer(nPed[, 3] - 1), # [[2]] sire ID/No.
as.integer(S@i), # [[3]] S@i
as.integer(S@p), # [[4]] S@p
as.double(S@x), # [[5]] S@x
as.integer(N), # [[6]] No. in pedigree
as.double(rep(0, length(S@x))), # [[7]] Sd@x
as.integer(rep(0, length(S@i))), # [[8]] Sd@i
as.integer(rep(0, N)), # [[9]] Sd@p
as.integer(0), # [[10]] cnt/count
as.integer(sex)) # [[11]] sex
Sd <- sparseMatrix(i = Cout[[8]][1:Cout[[10]]],
p = Cout[[9]][1:(nhom+1)],
x = Cout[[7]][1:Cout[[10]]],
dims = c(nhom, nhom),
dimnames = list(as.character(pedigree[homs, 1]), NULL),
symmetric = TRUE, index1 = FALSE)
diag(Sd) <- 1 - Sout$inbreeding[homs]
if(!returnS) S <- NULL
rm("Cout")
} else{
#TODO
stop("code not yet written to parallelize function") #FIXME
# listA <- data.frame(Row = as.integer(rep(1:length(A@p[-1]), diff(A@p))), Column = as.integer(A@i + 1))
# wrap_dij <- function(x){
# sub_lA <- listA[min(x):max(x), 1:2]
# lA_r <- dim(sub_lA)[1]
# Cout <- .C("dijp", PACKAGE = "nadiv",
# as.integer(numeric.pedigree[, 2] - 1),
# as.integer(numeric.pedigree[, 3] - 1),
# as.integer(lA_r),
# as.integer(sub_lA[, 1] - 1),
# as.integer(sub_lA[, 2] - 1),
# as.integer(A@i),
# as.integer(A@p),
# as.double(A@x/2),
# as.double(rep(0, lA_r)))
# Cout[[9]]
# }
# if(verbose) cat("starting to make D...")
# Dijs <- parallel::pvec(seq(1, dim(listA)[1], 1), FUN = wrap_dij, mc.set.seed = FALSE, mc.silent = FALSE, mc.cores = ncores, mc.cleanup = TRUE)
# D <- sparseMatrix(i = A@i, p = A@p, x = Dijs, ...)
# if(!returnA) A <- NULL
# D <- drop0(D)
# diag(D) <- 2 - dA
}
if(verbose) cat(".done", "\n")
if(det) logDet <- determinant(Sd, logarithm = TRUE)$modulus[1] else logDet <- NULL
if(invertSd){
if(verbose) cat("starting to invert Sd...")
Sdinv <- as(solve(Sd), "dgCMatrix")
Sdinv@Dimnames <- Sd@Dimnames
if(verbose) cat(".done", "\n")
listSdinv <- sm2list(Sdinv, rownames = Sd@Dimnames[[1L]],
colnames = c("row", "column", "Sdinverse"))
return(list(S = S, Sd = Sd, logDet = logDet,
Sdinv = Sdinv, listSdinv = listSdinv))
} else{
return(list(S = S, Sd = Sd, logDet = logDet))
}
}