require(BayesMendel) modified_brcapro <- function (family, counselee.id = 1, race = "Unknown", germline.testing = NULL, marker.testing = NULL, oophorectomy = NULL, mastectomy = NULL, params = brcaparams(), print = FALSE, imputeAges = FALSE, imputeRelatives = FALSE, mrate1, mrate2, f_index, m_index, gene, mod) # mrate1: de novo mutation rate of BRCA1, mrate2: de novo mutation rate of BRCA2; f_index: index of the counselee's father; m_index: index of the counselee's mother; gene: {1,2}, 1 = BRCA1, 2 = BRCA2; mod: {0,1,2}, 0 = no modification; 1: modify the likelihood of BRCA1, 2: BRCA2(2); fout: name of output file { warnmsg = "" if (imputeRelatives == TRUE & imputeAges == FALSE) { imputeAges = TRUE warnmsg = paste(warnmsg, "Warning: imputeAges was input as false, but has been set to TRUE. If imputeRelatives=TRUE then by default imputeAges must also be TRUE.") } if (nchar(warnmsg) > 0) { print(warnmsg) } family <- CheckFamStructure("brcapro", family, counselee.id, germline.testing, marker.testing, oophorectomy, mastectomy, imputeAges, imputeRelatives, params) if (is.character(family)) { return(family) } if (sum(family$Twins) > 0) { inputfamily <- family if (!is.null(germline.testing)) { inputfamily <- data.frame(inputfamily, germline.testing) } if (!is.null(marker.testing)) { inputfamily <- data.frame(inputfamily, marker.testing) } if (!is.null(oophorectomy)) { inputfamily <- data.frame(inputfamily, oophorectomy) } if (!is.null(mastectomy)) { inputfamily <- data.frame(inputfamily, mastectomy) } } originalcounselee.id <- counselee.id if (family$Twins[family$ID == counselee.id] > 0) { twins <- family[family$Twins == family$Twins[family$ID == counselee.id], ] counselee.id <- twins$ID[1] } proband.current.age <- max(family[which(family$ID == counselee.id), c("AgeBreast", "AgeOvary", "AgeBreastContralateral")]) nonAJ <- list(c(1 - 0.0005829, 0.0005829), c(1 - 0.000676, 0.000676)) AJ <- list(c(1 - 0.00609756097560976, 0.00609756097560976), c(1 - 0.00679723502304147, 0.00679723502304147)) Italian <- list(c(1 - 0.001673779, 0.001673779), c(1 - 0.00132622, 0.00132622)) Other <- params$allef allef <- list() if (any(family$ethnic == "AJ" & !is.na(family$ethnic))) allef["AJ"] <- list(AJ) if (any(family$ethnic == "nonAJ" & !is.na(family$ethnic))) allef["nonAJ"] <- list(nonAJ) if (any(family$ethnic == "Italian" & !is.na(family$ethnic))) allef["Italian"] <- list(Italian) if (any(family$ethnic == "Other" & !is.na(family$ethnic))) allef["Other"] <- list(Other) if (race == "Unknown") { params$penetrance <- params$penetrance } else { if (ethnic == "AJ") { params$penetrance$fFX[, "B00"] <-BRCApenet.AJ.2001.2008$fFX[, "B00"] params$penetrance$fFY[, "B00"] <- BRCApenet.AJ.2001.2008$fFY[, "B00"] params$penetrance$fMX[, "B00"] <- BRCApenet.AJ.2001.2008$fMX[, "B00"] params$penetrance$fMY[, "B00"] <- BRCApenet.AJ.2001.2008$fMY[, "B00"] } else { if (race == "nonAJ") { params$penetrance$fFX[, "B00"] <- BRCApenet.nonAJ.2001.2008$fFX[, "B00"] params$penetrance$fFY[, "B00"] <- BRCApenet.nonAJ.2001.2008$fFY[, "B00"] params$penetrance$fMX[, "B00"] <- BRCApenet.nonAJ.2001.2008$fMX[, "B00"] params$penetrance$fMY[, "B00"] <- BRCApenet.nonAJ.2001.2008$fMY[, "B00"] } else { if (race == "Hispanic") { params$penetrance$fFX[, "B00"] <- BRCAbaseline.race.2008$fFX[, "Hispanic"] params$penetrance$fFY[, "B00"] <- BRCAbaseline.race.2008$fFY[, "Hispanic"] params$penetrance$fMX[, "B00"] <- BRCAbaseline.race.2008$fMX[, "Hispanic"] params$penetrance$fMY[, "B00"] <- BRCAbaseline.race.2008$fMY[, "Hispanic"] } else { if (race == "NativeAmerican") { params$penetrance$fFX[, "B00"] <- BRCAbaseline.race.2008$fFX[, "NativeAmerican"] params$penetrance$fFY[, "B00"] <- BRCAbaseline.race.2008$fFY[, "NativeAmerican"] params$penetrance$fMX[, "B00"] <- BRCAbaseline.race.2008$fMX[, "NativeAmerican"] params$penetrance$fMY[, "B00"] <- BRCAbaseline.race.2008$fMY[, "NativeAmerican"] } else { if (race == "White") { params$penetrance$fFX[, "B00"] <- BRCAbaseline.race.2008$fFX[, "White"] params$penetrance$fFY[, "B00"] <- BRCAbaseline.race.2008$fFY[, "White"] params$penetrance$fMX[, "B00"] <- BRCAbaseline.race.2008$fMX[, "White"] params$penetrance$fMY[, "B00"] <- BRCAbaseline.race.2008$fMY[, "White"] } } } } } } psize <- dim(family)[1] nloci <- 2 ngen <- 3^nloci maxages <- apply(family[, c("AgeBreast", "AgeOvary", "AgeBreastContralateral")], 1, max) sex <- family$Gender currfamily <- family if (sum(family$Twins) > 0) { currfamily <- inputfamily } femaleproband <- ifelse(family$Gender[family$ID == counselee.id] == 0, 1, 0) if (family$AffectedBreast[family$ID == counselee.id] > 0) { firstdx.under40 <- ifelse(family$AgeBreast[family$ID == counselee.id] < 40, 1, 0) agedx <- family$AgeBreast[family$ID == counselee.id] if (femaleproband == 1 & firstdx.under40 == 1) { Breast <- params$CBCpenetrance$fFX.Under40 Ovarian <- params$penetrance$fFY } if (femaleproband == 1 & firstdx.under40 == 0) { Breast <- params$CBCpenetrance$fFX.Over40 Ovarian <- params$penetrance$fFY } if (femaleproband == 0 & firstdx.under40 == 1) { Breast <- params$CBCpenetrance$fMX.Under40 Ovarian <- params$penetrance$fMY } if (femaleproband == 0 & firstdx.under40 == 0) { Breast <- params$CBCpenetrance$fMX.Over40 Ovarian <- params$penetrance$fMY } CBC <- matrix(0, nrow = 110, ncol = 9) for (cccc in 1:ncol(Breast)) { CBC[agedx:110, cccc] <- Breast[1:length(agedx:110), cccc] } } else { if (femaleproband == 1) { CBC <- params$penetrance$fFX Ovarian <- params$penetrance$fFY } if (femaleproband == 0) { CBC <- params$penetrance$fMX Ovarian <- params$penetrance$fMY } } CBRCApenet <- list(CBC, Ovarian) names(CBRCApenet) <- c("fX", "fY") colnames(CBRCApenet$fX) <- colnames(CBRCApenet$fY) <- colnames(params$penetrance$fFX) if (ncol(params$penetrance[[1]]) == 9) { params$penetrance <- MakePenetPostIntervention(params$penetrance, params$PostOophorectomyHR, params$PostMastectomyHR, contralateral = FALSE) } else { if (ncol(params$penetrance[[1]]) != 27) { stop("Penetrance matrices must have 9 columns, or 27 if they include interventions") } } if (ncol(params$CBCpenetrance[[1]]) == 9) { params$CBCpenetrance <- MakePenetPostIntervention(params$CBCpenetrance, params$PostOophorectomyHR, params$PostMastectomyHR, contralateral = TRUE) } else { if (ncol(params$CBCpenetrance[[1]]) != 27) { stop("Penetrance matrices must have 9 columns, or 18 if they include interventions") } } missingage = any(family$AgeBreast == 1 | family$AgeOvary == 1 | (family$AgeBreastContralateral == 1 & family$AffectedBreast == 2)) if (missingage & imputeAges == FALSE) { agewarning = "Warning: Age imputation has been turned off, but there are unknown ages of some unaffected family members. You may want to get more information about family member ages and re-run the calculation, or set imputeAges=TRUE." post = runPeeling(family = family, params = params, ngen = ngen, oophorectomy = oophorectomy, mastectomy = mastectomy, germline.testing = germline.testing, marker.testing = marker.testing, psize = psize, counselee.id = counselee.id, allef = allef,mod=mod,gene=gene,nloci=nloci,mrate1=mrate1, mrate2=mrate2, f_index=f_index,m_index=m_index) } if (missingage & imputeAges == TRUE) { agewarning = "Warning: Unknown ages of some unaffected and affected family members have been imputed. You may want to get more information about family member ages and re-run the calculation." ifamily = ImputeAge(fff = family, params = params, model = "brcapro") postList = NULL nIter = params$nIter if (length(ifamily$mem_bc) > 0 | length(ifamily$mem_oc) > 0 | length(ifamily$mem_bcc) > 0 | length(ifamily$mem_mast) > 0 | length(ifamily$mem_ooph) > 0) { for (iIter in 1:nIter) { nf = ifamily$fff noophorectomy = oophorectomy nmastectomy = mastectomy if (length(ifamily$mem_bc) > 0) { nf$AgeBreast[ifamily$mem_bc] = ifamily$age_bc[iIter, ifamily$mem_bc] } if (length(ifamily$mem_oc) > 0) { nf$AgeOvary[ifamily$mem_oc] <- ifamily$age_oc[iIter, ifamily$mem_oc] } if (length(ifamily$mem_bcc) > 0) { nf$AgeBreastContralateral[ifamily$mem_bcc] <- ifamily$age_bcc[iIter, ifamily$mem_bcc] } if (length(ifamily$mem_mast) > 0) { nf$AgeMastectomy[ifamily$mem_mast] = ifamily$age_mast[iIter, ifamily$mem_mast] nmastectomy$AgeMastectomy[ifamily$mem_mast] = ifamily$age_mast[iIter, ifamily$mem_mast] } if (length(ifamily$mem_ooph) > 0) { nf$AgeOophorectomy[ifamily$mem_ooph] = ifamily$age_ooph[iIter, ifamily$mem_ooph] noophorectomy$AgeOophorectomy[ifamily$mem_ooph] = ifamily$age_ooph[iIter, ifamily$mem_ooph] } #postList = rbind(postList, runPeeling(family = nf, # params = params, ngen = ngen, oophorectomy = noophorectomy, # mastectomy = nmastectomy, germline.testing = germline.testing, # marker.testing = marker.testing, psize = psize, # counselee.id = counselee.id, allef = allef)) postList = rbind(postList, runPeeling(family = nf, params = params, ngen = ngen, oophorectomy = oophorectomy, mastectomy = mastectomy, germline.testing = germline.testing, marker.testing = marker.testing, psize = psize, counselee.id = counselee.id, allef = allef,mod=mod,gene=gene,nloci=nloci,mrate1=mrate1, mrate2=mrate2, f_index=f_index,m_index=m_index)) } post = apply(postList, 2, mean) } if (length(ifamily$mem_bc) == 0 & length(ifamily$mem_oc) == 0 & length(ifamily$mem_bcc) == 0 & length(ifamily$mem_mast) == 0 & length(ifamily$mem_ooph) == 0) { #post = runPeeling(family = ifamily$fff, params = params, # ngen = ngen, oophorectomy = oophorectomy, mastectomy = mastectomy, # germline.testing = germline.testing, marker.testing = marker.testing, # psize = psize, counselee.id = counselee.id, allef = allef) post = runPeeling(family = ifamily$fff, params = params, ngen = ngen, oophorectomy = oophorectomy, mastectomy = mastectomy, germline.testing = germline.testing, marker.testing = marker.testing, psize = psize, counselee.id = counselee.id, allef = allef,mod=mod,gene=gene,nloci=nloci,mrate1=mrate1, mrate2=mrate2, f_index=f_index,m_index=m_index) } } if (!missingage) { agewarning = "" post = runPeeling(family = family, params = params, ngen = ngen, oophorectomy = oophorectomy, mastectomy = mastectomy, germline.testing = germline.testing, marker.testing = marker.testing, psize = psize, counselee.id = counselee.id, allef = allef,mod=mod,gene=gene,nloci=nloci,mrate1=mrate1, mrate2=mrate2, f_index=f_index,m_index=m_index) } if (!missingage | imputeAges == FALSE) { family = family } if (missingage & imputeAges == TRUE) { family = ifamily$fff } ll <- NULL post <- matrix(post, nrow = 3, ncol = 3, byrow = FALSE) post <- post[-3, -3] rownames(post) <- c("BRCA10", "BRCA11") colnames(post) <- c("BRCA20", "BRCA21") ages.cancer <- c(family[which(family$ID == counselee.id), "AgeBreast"], family[which(family$ID == counselee.id), "AgeOvary"], family[which(family$ID == counselee.id), "AgeBreastContralateral"]) current.age <- max(ages.cancer) ages <- seq(proband.current.age + params$age.by, ifelse(params$age.to >= proband.current.age + params$age.by, params$age.to, 110), by = params$age.by) p0 <- post["BRCA10", "BRCA20"] p2 <- sum(post["BRCA10", ]) - p0 p1 <- sum(post[, "BRCA20"]) - p0 p12 <- 1 - p1 - p2 - p0 if (print == TRUE) { cat("The probability of being a carrier is", 1 - p0, "\n an BRCA1 carrier", p1, "\n an BRCA2 carrier", p2, "\n") } return(probs=c(p0=p0,p1=p1,p2=p2,p12=p12)) }