R/genotype_combinations_f.R In heritEWAS: Identify Heritable Methylation Marks

Defines functions genotype_combinations_f

```genotype_combinations_f <- function(fam, ncores = 1, cl = NULL, verbose) {

# Make a copy of the original ID
fam\$ID.orig <- fam\$indiv

# Trim the pedigrees w.r.t typed people
fam <- trim_pedigree(fam, reorder = TRUE)

# Check if only one typed person in the family
if (nrow(fam) == 1) {
out <- data.frame(p = c(1, 1), x = c(0, 1))
names(out)[2] <- fam\$ID.orig[1]
return(out)
}

# Find the founders, typed people and their intersection
founders <- fam\$indiv[fam\$mother == 0 & fam\$father == 0]
typed <- fam\$indiv[fam\$typed == 1]
common <- intersect(founders, typed)

# Possible genotypes for typed people
typed.geno <- generate_typed_geno(fam)

if (verbose) {
cat("Number of genotype combinations to check:", nrow(typed.geno), "\n")
}

# Possible genotypes for founders to order 1
founder.geno.deg1 <- diag(length(founders))

# Possible genotype for remaining people who are neither founders nor typed
remaining <- setdiff(1:nrow(fam), c(typed, founders))
if (length(remaining) == 0) {
remaining.geno <- matrix(0, 0, 0)
} else {
remaining.geno <- geno.comb(0:1, length(remaining))
}

# Create matrix of transmission probablities
trans.matrix <- create_trans_matrix()

# Calculate the transmission probabilities
trans.prob <- function(g, trans.matrix, fam) {
# Given a vector g = (g1, ..., gn) of unphased genotypes,
# one for each member of fam with each gi = 0, 1 or 2,
# calculate the transmission probabilities, i.e. P(g)/Prior(g).
# This function assumes the people in fam are numbered consecutively
# and either the mum and dad IDs are both 0 or neither are.

# Reorder the people in fam consecutively
# fam <- convert.IDs(cbind(family = 1, fam),
#                    convert.IDs.numeric = TRUE)[, -1]
gm <- g[fam\$mother]
gd <- g[fam\$father]
go <- g[fam\$mother != 0 & fam\$father != 0]
y <- 1 + 9 * go + 3 * gm + gd
p <- trans.matrix\$p[y]
prod(p)
}

# Calculate the probability for each genotype combinaton
p.trans.deg1 <- function(gtyped) {
gtyped <- as.numeric(gtyped)
g <- rep(-1, nrow(fam))
g[typed] <- gtyped
trans <- 0

# If the founders and typed people overlap then restrict founder.geno.deg1
# so it is consistent with gtyped
fgd1 <- founder.geno.deg1

if (length(common) > 0) {
keep <- rep(TRUE, nrow(fgd1))
for (i in 1:length(common)) {
j <- which(founders == common[i])[1]
k <- which(typed == common[i])[1]
keep <- keep & (fgd1[, j] == gtyped[k])
}
if (sum(keep) == 0) {
return(0)
}
fgd1 <- fgd1[keep, ]

if (sum(keep) == 1) fgd1 <- t(as.matrix(as.vector(fgd1)))
}

if (length(remaining) == 0) {
for (i in 1:nrow(fgd1)) {
g[founders] <- fgd1[i, ]
trans <- trans + trans.prob(g, trans.matrix, fam)
}
} else if (nrow(remaining.geno) > 0) {
for (i in 1:nrow(fgd1)) {
for (j in 1:nrow(remaining.geno)) {
g[founders] <- fgd1[i, ]
g[remaining] <- as.numeric(remaining.geno[j, ])
trans <- trans + trans.prob(g, trans.matrix, fam)
}
}
}
# Prior(g) = 2 * p^k * (1-p)^(2*f-k) when k==1
trans <- 2 * trans
trans
}

if (ncores == 1) {
prob <- apply(typed.geno, 1, p.trans.deg1)
} else {
#cl <- parallel::makeCluster(ncores)
prob <- parallel::parApply(cl, typed.geno, 1, p.trans.deg1)
}

# Output file
out <- data.frame(prob[prob != 0], typed.geno[prob != 0, ])
names(out) <- c("p", fam\$ID.orig[typed])
row.names(out) <- NULL

out

}
```

Try the heritEWAS package in your browser

Any scripts or data that you put into this service are public.

heritEWAS documentation built on July 1, 2020, 6:02 p.m.