Famdenovo.BRCA <- function(BRCA.inp, denovo_ratio1, denovo_ratio2) {
  
  csv <- BRCA.inp
  csv.fam <- data.frame(csv$FamilyID, csv$BRCA1carrier, csv$BRCA2carrier, csv$ID, csv$Gender, csv$FatherID, csv$MotherID, csv$AffectedBreast, csv$AffectedOvary, csv$AgeBreast, csv$AgeOvary, csv$AgeBreastContralateral, csv$Twins, csv$ethnic)
  colnames(csv.fam) <- c("FamilyID", "BRCA1carrier", "BRCA2carrier", "ID", "Gender", "FatherID", "MotherID", "AffectedBreast", "AffectedOvary", "AgeBreast", "AgeOvary", "AgeBreastContralateral", "Twins", "ethnic")
  csv.fam0 <- csv.fam
  
  #denovo_ratioAJ<-0.5*denovo_ratio
  #denovo_ratioAJ <- denovo_ratio
  if (csv.fam0$ethnic[1]=="AJ"){
    mrate1 <- 0.00609756097560976 * denovo_ratio1
    mrate2 <- 0.00679723502304147 * denovo_ratio2
  }else{
    mrate1 <- 0.0005829 * denovo_ratio1
    mrate2 <- 0.000676 * denovo_ratio2
  }
  families <- csv.fam$FamilyID[!duplicated(csv.fam$FamilyID)]
  agewarninglist=list()
  out2 <- t(c(9, -9.9, "BRCA2", 9.9))  # the id is set to -9.9 and the denovo_probability is set to 9.9, which is an example, will be deleted in the output
  colnames(out2) <- c("FamilyID", "ID", "gene", "denovo_probability")
  out1 <- t(c(9, -9.9, "BRCA1", 9.9))
  colnames(out1) <- c("FamilyID", "ID", "gene", "denovo_probability")
  
  for(n in 1:length(families)) {             # loop by family
    cat(n, '\n')
    csv.fam <- csv.fam0[csv.fam0$FamilyID == families[n], 4:ncol(csv.fam0)]             
    csv.fam <- csv.fam[!is.na(csv.fam$Gender), ]
    colnames(csv.fam) <- c("ID", "Gender", "FatherID", "MotherID", "AffectedBreast", "AffectedOvary", "AgeBreast", "AgeOvary", "AgeBreastContralateral", "Twins", "ethnic")
    csv.fam[,11] <- as.character(csv.fam[,11])
    ## recode individual id; NA = 0; otherwise an integer
    # pick individuals with NA father
    fid <- rep(0, times = nrow(csv.fam))
    index_f <- substr(csv.fam[,3], nchar(csv.fam[,3])-1, nchar(csv.fam[,3])) == "_0"
    fid[index_f] <- 0                   # father NA
    # pick individuals with NA mother
    mid <- rep(0, nrow(csv.fam))
    index_m <- substr(csv.fam[,4], nchar(csv.fam[,4])-1, nchar(csv.fam[,4])) == "_0"
    mid[index_m] <- 0                   # mother NA
    # recode the kid
    id <- rep(0, times = nrow(csv.fam))
    id <- 1:nrow(csv.fam)
    id <- as.double(id)
    # transform the remaining ids in father and mother
    for(i in 1:nrow(csv.fam)) {
      tp_index <- csv.fam[i,3] == csv.fam[,1]
      if(sum(tp_index) >= 1) {
        fid[i] <- id[tp_index][1]
      }
      tp1_index <- csv.fam[i,4] == csv.fam[,1]
      if(sum(tp1_index) >= 1) {
        mid[i] <- id[tp1_index][1]
      }
    }
    CGI.fam <- data.frame(id, csv.fam$Gender, fid, mid, csv.fam[,5:ncol(csv.fam)])
    colnames(CGI.fam) <- c("ID", "Gender", "FatherID", "MotherID", "AffectedBreast", "AffectedOvary", "AgeBreast", "AgeOvary", "AgeBreastContralateral", "Twins", "ethnic")
    data(brca.fam,BRCApenet.metaDSL.2008,death.othercauses, compriskSurv, CBRCApenet.2012, BrOvJointDsn.2014)
    ## loop for mutation carriers in the family
    temp <- csv.fam0[csv.fam0$FamilyID == families[n], ]
    temp <- temp[!is.na(temp$Gender), ]
    temp$BRCA1carrier[is.na(temp$BRCA1carrier)] <- -9
    temp$BRCA2carrier[is.na(temp$BRCA2carrier)] <- -9
    index1 <- temp$BRCA1carrier == 1 
    index2 <- temp$BRCA2carrier == 1
    # BRCA2 mutation carriers
    if(sum(index2) > 0) {
      for(k in 1:sum(index2)) {
        #is_familial2 <- 0  # mark familial state as 0
        # mutation carrier is the counselee
        cs0 <- temp$ID[index2][k] # counselee id
        index_counselee <- csv.fam$ID == cs0
        cs <- CGI.fam$ID[index_counselee]# recoded counselee id for brcapro
        ## if father or mother not available, next
        if(CGI.fam$FatherID[index_counselee] == 0 || CGI.fam$MotherID[index_counselee] == 0) {
          next
        }
        ## if the father or the mother is a carrier, mark it as familial
        ##if(temp$BRCA2carrier[temp$ID == temp$FatherID[temp$ID == cs0]] == 1 || temp$BRCA2carrier[temp$ID == temp$MotherID[temp$ID == cs0]] == 1) {
        ##  is_familial2 = 1
        ##}
        
        ## calculate p1 = P(Gc=1|Gm=0,Gf=0,D,P)
        tp_modified_brcapro <- modified_brcapro(CGI.fam, counselee.id = cs, 
                                                mrate1 = mrate1, mrate2 = mrate2, 
                                                f_index = which(CGI.fam$ID == CGI.fam$FatherID[index_counselee]), 
                                                m_index = which(CGI.fam$ID == CGI.fam$MotherID[index_counselee]), 
                                                gene=2, mod = 2, print = F)
        #agewarninglist[length(agewarninglist)+1]<-tp_modified_brcapro[[2]]
        #p1 <- tp_modified_brcapro[[1]]@probs[3] + tp_modified_brcapro[[1]]@probs[4]
        p1 <- 1-tp_modified_brcapro[1]
        
        ## calculate p2 = P(Gf=0|D,P) 
        tp_modified_brcapro <- modified_brcapro(CGI.fam, counselee.id = CGI.fam$FatherID[index_counselee], 
                                                mrate1 = mrate1, mrate2 = mrate2, 
                                                f_index = NULL, 
                                                m_index = NULL, 
                                                gene=2, mod = 0, print = F)
        #agewarninglist[length(agewarninglist)+1]<-tp_modified_brcapro[[2]]
        #p2 <- 1-tp_modified_brcapro[[1]]@probs[3]-tp_modified_brcapro[[1]]@probs[4]
        p2 <- tp_modified_brcapro[1]
        
        ## calculate p3 = P(Gm=0|D,P)
        tp_modified_brcapro <- modified_brcapro(CGI.fam, counselee.id = CGI.fam$MotherID[index_counselee], 
                                                mrate1 = mrate1, mrate2 = mrate2, 
                                                f_index = NULL, 
                                                m_index = NULL, 
                                                gene=2, mod = 0, print = F)
        #agewarninglist[length(agewarninglist)+1]<-tp_modified_brcapro[[2]]
        #p3 <- 1-tp_modified_brcapro[[1]]@probs[3]-tp_modified_brcapro[[1]]@probs[4]
        p3 <- tp_modified_brcapro[1]
        
        ## calculate p4 = P(Gc=1|D,P)
        tp_modified_brcapro <- modified_brcapro(CGI.fam, counselee.id = cs, 
                                                mrate1 = mrate1, mrate2 = mrate2, 
                                                f_index = NULL, 
                                                m_index = NULL, 
                                                gene=2, mod = 0, print = F)
        #agewarninglist[length(agewarninglist)+1]<-tp_modified_brcapro[[2]]
        #p4 = tp_modified_brcapro[[1]]@probs[3] + tp_modified_brcapro[[1]]@probs[4]
        p4 = 1-tp_modified_brcapro[1]
        
        ## the de novo probability p_BRCA2_denovo
        #print(c(p1=p1,p2=p2,p3=p3,p4=p4))
        p_BRCA2_denovo = p1*p2*p3/p4
        tp <- t(c(families[n], cs0, "BRCA2_denovo_prob", p_BRCA2_denovo))
        colnames(tp) <- colnames(out2)
        out2 <- rbind(out2, tp)
      }
    } #else {
      #cat("no BRCA2 mutation carriers\n")
    #}
    # BRCA1 mutation carriers
    if(sum(index1) > 0) {
      for(k in 1:sum(index1)) {
        #is_familial1 <- 0
        # mutation carrier is the counselee
        cs0 <- temp$ID[index1][k] # counselee id
        index_counselee <- csv.fam$ID == cs0
        cs <- CGI.fam$ID[index_counselee]# recoded counselee id for brcapro
        ## if father or mother not available, next
        if(CGI.fam$FatherID[index_counselee] == 0 || CGI.fam$MotherID[index_counselee] == 0) {
          next
        }
        ## if the father or mother is a carrier, mark is_familial1 = 1
        #if(temp$BRCA1carrier[temp$ID == temp$FatherID[temp$ID == cs0]] == 1 || temp$BRCA1carrier[temp$ID == temp$MotherID[temp$ID == cs0]] == 1) {
        #  is_familial1 = 1
        #}
        
        ## calculate p1 = P(Gc=1|Gm=0,Gf=0,D,P)
        tp_modified_brcapro <- modified_brcapro(CGI.fam, counselee.id = cs, mrate1 = mrate1, mrate2 = mrate2, 
                                                f_index = which(CGI.fam$ID == CGI.fam$FatherID[index_counselee]), 
                                                m_index = which(CGI.fam$ID == CGI.fam$MotherID[index_counselee]), 
                                                gene = 1, mod = 1, print = F)
        #agewarninglist[length(agewarninglist)+1]<-tp_modified_brcapro[[2]]
        #p1 <- tp_modified_brcapro[[1]]@probs[2] + tp_modified_brcapro[[1]]@probs[4]
        p1 <- 1-tp_modified_brcapro[1]
        #if(p1 < 0) {
        #  write.table(tp_modified_brcapro[[1]]@probs, file = "minus_0.txt", quote = F, col.names = F, row.names = F, append = F)
        #}
        
        ## calculate p2 = P(Gf=0|D,P) 
        tp_modified_brcapro <- modified_brcapro(CGI.fam, counselee.id = CGI.fam$FatherID[index_counselee], 
                                                mrate1 = mrate1, mrate2 = mrate2, 
                                                f_index = NULL, 
                                                m_index = NULL, 
                                                gene = 1, mod = 0, print = F)
        #agewarninglist[length(agewarninglist)+1]<-tp_modified_brcapro[[2]]
        p2 <- tp_modified_brcapro[1]
        #p2 <- 1-tp_modified_brcapro[[1]]@probs[2]-tp_modified_brcapro[[1]]@probs[4]
        #if(p2 < 0) {
        #  write.table(c("p2", p2), file = "minus_0.txt", append = T)
        #}
        
        ## calculate p3 = P(Gm=0|D,P)
        tp_modified_brcapro <- modified_brcapro(CGI.fam, counselee.id = CGI.fam$MotherID[index_counselee], 
                                                mrate1 = mrate1, mrate2 = mrate2, 
                                                f_index = NULL, 
                                                m_index = NULL, 
                                                gene = 1, mod = 0, print = F)
        #agewarninglist[length(agewarninglist)+1]<-tp_modified_brcapro[[2]]
        p3 <- tp_modified_brcapro[1]
        #p3 <- 1-tp_modified_brcapro[[1]]@probs[2]-tp_modified_brcapro[[1]]@probs[4]
        #print(tp_modified_brcapro[[1]]@probs)
        #if(p3 < 0) {
        #  write.table(c("p3", p3), file = "minus_0.txt", append = T)
        #}
        
        ## calculate p4 = P(Gc=1|D,P)
        tp_modified_brcapro <- modified_brcapro(CGI.fam, counselee.id = cs, 
                                                mrate1 = mrate1, mrate2 = mrate2,
                                                f_index = NULL, 
                                                m_index = NULL, 
                                                gene = 1, mod = 0, print = F)
        #agewarninglist[length(agewarninglist)+1]<-tp_modified_brcapro[[2]]
        p4 <- 1-tp_modified_brcapro[1]
        #print(tp_modified_brcapro[[1]]@probs)
        #p4 = tp_modified_brcapro[[1]]@probs[2] + tp_modified_brcapro[[1]]@probs[4]
        #if(p4 < 0) {
        #  write.table(c("p4", p4), file = "minus_0.txt", append = T)
        #}
        
        ## the de novo probability p_BRCA1_denovo
        ##print(c(p1=p1,p2=p2,p3=p3,p4=p4))
        p_BRCA1_denovo = p1*p2*p3/p4
        tp <- t(c(families[n], cs0, "BRCA1_denovo_prob", p_BRCA1_denovo))
        colnames(tp) <- colnames(out1)
        out1 <- rbind(out1, tp)
      }
    } #else {
      #cat("no BRCA1 mutation carriers\n")
    #}
  }
  out <- rbind(out1, out2)
  if(nrow(out) == 2) {
    cat("no BRCA1 or BRCA2 mutation carriers with parents available in input\n")
  } else {
    out <- out[out[,2] != -9.9 & out[,4] != 9.9, ]
  }
  x = as.matrix(out)
  if (ncol(x)==1) x = t(x)
  return(x)
}