From 73ddd72f911a5af55b9482612d8f439b74c38e55 Mon Sep 17 00:00:00 2001 From: Sam Gurr Date: Mon, 23 Dec 2024 15:06:30 -0500 Subject: [PATCH] updated parentage script --- RAnalysis/Scripts/Popgen/.Rhistory | 1004 ++++++++++----------- RAnalysis/Scripts/Popgen/05_parentage.Rmd | 165 +++- 2 files changed, 628 insertions(+), 541 deletions(-) diff --git a/RAnalysis/Scripts/Popgen/.Rhistory b/RAnalysis/Scripts/Popgen/.Rhistory index cd51720..eac8bbf 100644 --- a/RAnalysis/Scripts/Popgen/.Rhistory +++ b/RAnalysis/Scripts/Popgen/.Rhistory @@ -1,512 +1,512 @@ -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 20, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -2, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.5, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE +F0BF1O.GenoM_final, # GenoM +F0BF1O.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1 # n.bootstrap ) -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 25, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -2, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.5, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE +# F0 AND F1 OFFPSRING ALL :::::::::::::::::::::::::::::::::::::::::: +F0BF1O.sequoia <- BootPed( +F0BF1O.GenoM_final, # GenoM +F0BF1O.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +100, # n.SNP +1 # n.bootstrap ) -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 25, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -2, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.5, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE +# F0 AND F1 OFFPSRING ALL :::::::::::::::::::::::::::::::::::::::::: +F0BF1O.sequoia <- BootPed( +F0BF1O.GenoM_final, # GenoM +F0BF1O.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.01,# min.LLR +100, # n.SNP +1 # n.bootstrap ) -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 25, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -2, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.5, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE +# F0 AND F1 OFFPSRING ALL :::::::::::::::::::::::::::::::::::::::::: +F0BF1O.sequoia <- BootPed( +F0BF1O.GenoM_final, # GenoM +F0BF1O.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.001,# min.LLR +100, # n.SNP +1 # n.bootstrap ) -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -2, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.5, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE +# F0 AND F1 OFFPSRING ALL :::::::::::::::::::::::::::::::::::::::::: +F0BF1O.sequoia <- BootPed( +F0BF1O.GenoM_final, # GenoM +F0BF1O.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1 # n.bootstrap ) -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -1, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.5, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE +# F0 AND F1 OFFPSRING ALL :::::::::::::::::::::::::::::::::::::::::: +F0BF1O.sequoia <- BootPed( +F0BF1O.GenoM_final, # GenoM +F0BF1O.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +10 # n.bootstrap ) -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -0.2, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.5, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE +# F0 AND F1 OFFPSRING ALL :::::::::::::::::::::::::::::::::::::::::: +F0BF1O.sequoia <- BootPed( +F0BF1O.GenoM_final, # GenoM +F0BF1O.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +20 # n.bootstrap ) -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -0.2, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.5, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE +# low +F1BF2O_LOW.sequoia <- BootPed( +F1BF2O_LOW.GenoM_final, +F1BF2O_LOW.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +10 # n.bootstrap ) -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -5, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.5, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE +# low +F1BF2O_LOW.sequoia <- BootPed( +F1BF2O_LOW.GenoM_final, # GenoM +F1BF2O_LOW.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +20 # n.bootstrap ) -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -10, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.5, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE +# low +F1BF2O_LOW.sequoia <- BootPed( +F1BF2O_LOW.GenoM_final, # GenoM +F1BF2O_LOW.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +50 # n.bootstrap ) -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -5, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.1, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE -) -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -5, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.2, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE -) -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -5, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.2, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE -) -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -2, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.2, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE -) -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -2, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.2, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE -) -test2 <- F0BF1J.bootstrap.sequoia$PedigreePar %>% select(id, dam, LLRdam, sire, LLRsire) %>% na.omit() -test2 -# Call the cumulative dataframe that we will write to in the for loop below -df_total <- data.frame() # start dataframe -pedigree.table <- data.frame(matrix(nrow = 1, ncol = 7)) # create dataframe to save cumunalitively during for loop -colnames(pedigree.table) <- c('id', 'dam', 'LLRdam', 'sire' , 'LLRsire') # names for comuns in the for loop -# Call the cumulative dataframe that we will write to in the for loop below -df_total <- data.frame() # start dataframe -pedigree.table <- data.frame(matrix(nrow = 1, ncol = 6)) # create dataframe to save cumunalitively during for loop -colnames(pedigree.table) <- c('id', 'dam', 'LLRdam', 'sire' , 'LLRsire') # names for comuns in the for loop -for (i in 1:100) { -randomSNPs[i] <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs[i]] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -2, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.2, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE -) -pedigree_data <- data.frame(F0BF1J.bootstrap.sequoia$PedigreePar %>% select(id, dam, LLRdam, sire, LLRsire) %>% na.omit()) -df_total <- rbind(df_total,pedigree_data) #bind to a cumulative list dataframe -print(df_total) # print to monitor progress -} -randomSNPs[[1]] <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -randomSNPs_[`1`] <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -randomSNPs_[1] <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -randomSNPs_`1` <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -randomSNPs_1 <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), 50, replace = TRUE)))) -randomSNPs -for (i in 1:100) { -x <- (50+i) - i # assign i to 50, # randomSNPs == 50 - repeat 100 times (for loop) -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), x, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs[i]] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -2, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.2, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE -) -pedigree_data <- data.frame(F0BF1J.bootstrap.sequoia$PedigreePar %>% select(id, dam, LLRdam, sire, LLRsire) %>% na.omit()) -df_total <- rbind(df_total,pedigree_data) #bind to a cumulative list dataframe -print(df_total) # print to monitor progress -} -# Call the cumulative dataframe that we will write to in the for loop below -df_total <- data.frame() # start dataframe -for (i in 1:100) { -x <- (50+i) - i # assign i to 50, # randomSNPs == 50 - repeat 100 times (for loop) -randomSNPs <- array(sample(c(sample(1:ncol(F0BF1J.matrix), x, replace = TRUE)))) -bootstrap_GenoM <- F0BF1J.matrix[,randomSNPs] -F0BF1J.bootstrap.sequoia <- sequoia( -GenoM = bootstrap_GenoM, -LifeHistData = F0BF1J.LifeHistData.editted, #F0BF1J.LifeHistData, -SeqList = NULL, -# most comp intensive, run everything -Module = "par", -Err = 1e-04, -# threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated -Tfilter = -2, -# minimum LLR required for acceptance of proposed relationship relative to the next most likely -Tassign = 0.2, -# MaxSibshipSize = 100, -# breeding system complexity -Complex = "full", -Herm = "A", -UseAge = "extra", -CalcLLR = TRUE -) -pedigree_data <- data.frame(F0BF1J.bootstrap.sequoia$PedigreePar %>% select(id, dam, LLRdam, sire, LLRsire) %>% na.omit()) -df_total <- rbind(df_total,pedigree_data) #bind to a cumulative list dataframe -print(df_total) # print to monitor progress -} -df_total -# were we able to capture at least one pedigree for each of the F1 offspring? -length(colnames(All.vcf@gt[,63:139])) # number of F1 juvniles in the vcf /genind/ genoM matrix -length(unique(df_total$id)) -df_total -# get the est fit model -# The larger the absolute value of the LLR, the stronger the evidence for the preferred model -df_total$absLLRtotal <- abs(df_total$LLRdam) + abs(df_total$LLRsire) -df_total -df_total.byLLR <- df_total[,.SD[which.max(absLLRtotal)],by=id] -df_total -df_total.byLLR <- df_total %>% group_by(id) %>% dplyr::top_n(1, absLLRtotal) -df_total.byLLR -df_total.byLLR -pedigree_data <- df_total -# get the est fit model -# The larger the absolute value of the LLR, the stronger the evidence for the preferred model -pedigree_data$absLLRtotal <- abs(pedigree_data$LLRdam) + abs(pedigree_data$LLRsire) -pedigree_data.byLRRtotal <- pedigree_data %>% -group_by(id) %>% -dplyr::top_n(1, absLLRtotal) -pedigree_data.byLRRtotal -F0BF1J.pedigree <- pedigree_data %>% -group_by(id) %>% # group by id - there are mutliple occurances of pedigree assign for ids! -dplyr::top_n(1, absLLRtotal) # take the best hit essentiall, highest LLR total -F0BF1J.pedigree <- pedigree_data %>% -group_by(id) %>% # group by id - there are mutliple occurances of pedigree assign for ids! -dplyr::top_n(1, absLLRtotal) %>% # take the best hit essentiall, highest LLR total -dplyr::select(id, dam, sire) # select the target columns for next sequoia run -F0BF1J.pedigree -# Pairs - dataframe with columns ID1 and ID2 containing all possible combinations -F0BF1J.Pairs_raw <- data.frame(F0BF1J.ListIDs) %>% -dplyr::rename(ID1 = F0BF1J.ListIDs) %>% # rename -dplyr::mutate(ID2 = ID1) %>% # duplicate -tidyr::expand(ID1, ID2) %>% # get all possible ocmbinations -dplyr::mutate(AgeDif = case_when( -# ID1 and ID2 are F0 broodstock -(grepl("F0", ID1, ignore.case = TRUE) & -grepl("F0", ID2, ignore.case = TRUE)) ~ 0, -# ID1 and ID2 are F1 Juveniles -(!grepl("F0", ID1, ignore.case = TRUE) & -!grepl("F0", ID2, ignore.case = TRUE)) ~ 0, -# ID1 is F1 juvenile and ID2 is F0 broodstock -(!grepl("F0", ID1, ignore.case = TRUE) & -grepl("F0", ID2, ignore.case = TRUE)) ~ 1, -# IF1 is F0 broodstock and ID2 is F1 juvenile -(grepl("F0", ID1, ignore.case = TRUE) & -!grepl("F0", ID2, ignore.case = TRUE)) ~ -1, -)) -F0BF1J.Pairs <- F0BF1J.Pairs_raw[F0BF1J.Pairs_raw$ID1 != F0BF1J.Pairs_raw$ID2,] -# RUN IT! -CalcPairLL( -Pairs = F0BF1J.Pairs, -GenoM = F0BF1J.matrix, -Pedigree = F0BF1J.pedigree, -LifeHistData = F0BF1J.LifeHistData.editted, # working on it -AgePrior = TRUE, # working on it -SeqList = NULL, -Complex = "full", -Herm = "A", # "A" (distinguish between dam and sire role, default if at least 1 individual with sex=4), or "B" (no distinction between dam and sire role). -Err = 1e-04, -ErrFlavour = "version2.9", -Tassign = 0.5, -Tfilter = -2, -quiet = FALSE, -Plot = TRUE -) -F0BF1J.pedigree -F0BF1J.matrix -F0BF1J.Pairs -F0BF1J.LifeHistData.editted -# RUN IT! -CalcPairLL( -Pairs = F0BF1J.Pairs, -GenoM = F0BF1J.matrix, -Pedigree = F0BF1J.pedigree, -LifeHistData = F0BF1J.LifeHistData.editted, # working on it -AgePrior = TRUE, # working on it -SeqList = NULL, -Complex = "full", -Herm = "A", # "A" (distinguish between dam and sire role, default if at least 1 individual with sex=4), or "B" (no distinction between dam and sire role). -Err = 1e-04, -ErrFlavour = "version2.9", -Tassign = 0.5, -Tfilter = -2, -quiet = FALSE, -Plot = TRUE -) -class(F0BF1J.pedigree) -class(as.data.frame(F0BF1J.pedigree)) -CalcPairLL( -Pairs = F0BF1J.Pairs, -GenoM = F0BF1J.matrix, -Pedigree = as.data.frame(F0BF1J.pedigree), -LifeHistData = F0BF1J.LifeHistData.editted, # working on it -AgePrior = TRUE, # working on it -SeqList = NULL, -Complex = "full", -Herm = "A", # "A" (distinguish between dam and sire role, default if at least 1 individual with sex=4), or "B" (no distinction between dam and sire role). -Err = 1e-04, -ErrFlavour = "version2.9", -Tassign = 0.5, -Tfilter = -2, -quiet = FALSE, -Plot = TRUE -) -# view cl names to know what to subset by -colnames(All.vcf@gt) -# F3 juveniles are 178:250 -colnames(All.vcf@gt[,251:393]) # third generation of juveniles (offpsring of pH-specific parents, grandchildren of F1 pH parents, great graandchildren of F0) -# F3 juveniles are 178:250 -colnames(All.vcf@gt[,251:393]) # third generation of juveniles (offpsring of pH-specific parents, grandchildren of F1 pH parents, great graandchildren of F0) -# view cl names to know what to subset by -colnames(All.vcf@gt) -# F3 juveniles are 178:250 -colnames(All.vcf@gt[,251:392]) # third generation of juveniles (offpsring of pH-specific parents, grandchildren of F1 pH parents, great graandchildren of F0) -# view cl names to know what to subset by -colnames(All.vcf@gt) -# F1 broodstock 27:62 -colnames(All.vcf@gt[,27:62]) -# view cl names to know what to subset by -colnames(All.vcf@gt) -# F2 broodstock -colnames(All.vcf@gt[,140:177]) -F0BF1JF1BF2J.vcf <- All.vcf[,c(1:26,63:139,178:250, 27:62])] # now 175 samples! -F0BF1JF1BF2J.vcf <- All.vcf[,c(1:26,63:139,178:250,27:62])] # now 175 samples! +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +10)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +20)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +40)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-0.1,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-0.01,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.01,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.001,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.01,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.01,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.01,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.01,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.01,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.01,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.01,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.01,# min.LLR +30, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +25, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# moderate +F1BF2O_MOD.sequoia <- BootPed( +F1BF2O_MOD.GenoM_final, +F1BF2O_MOD.LifeHistData, # LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# low +F2BF3O_LOW.sequoia <- BootPed( +F2BF3O_LOW.GenoM_final, +F2BF3O_LOW.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +1)# n.bootstrap +# low +F2BF3O_LOW.sequoia <- BootPed( +F2BF3O_LOW.GenoM_final, +F2BF3O_LOW.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +10)# n.bootstrap +# moderate +F2BF3O_MOD.sequoia <- BootPed( +F2BF3O_MOD.GenoM_final, +F2BF3O_MOD.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +10)# n.bootstrap +# high +F2BF3O_HIGH.sequoia <- BootPed( +F2BF3O_HIGH.GenoM_final, +F2BF3O_HIGH.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +10)# n.bootstrap +# high +F2BF3O_HIGH.sequoia <- BootPed( +F2BF3O_HIGH.GenoM_final, +F2BF3O_HIGH.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +25, # n.SNP +10)# n.bootstrap +# high +F2BF3O_HIGH.sequoia <- BootPed( +F2BF3O_HIGH.GenoM_final, +F2BF3O_HIGH.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +25, # n.SNP +20)# n.bootstrap +# high +F2BF3O_HIGH.sequoia <- BootPed( +F2BF3O_HIGH.GenoM_final, +F2BF3O_HIGH.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +20)# n.bootstrap +# high +F2BF3O_HIGH.sequoia <- BootPed( +F2BF3O_HIGH.GenoM_final, +F2BF3O_HIGH.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +50)# n.bootstrap +# moderate +F2BF3O_MOD.sequoia <- BootPed( +F2BF3O_MOD.GenoM_final, +F2BF3O_MOD.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +5)# n.bootstrap +# moderate +F2BF3O_MOD.sequoia <- BootPed( +F2BF3O_MOD.GenoM_final, +F2BF3O_MOD.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +25, # n.SNP +5)# n.bootstrap +# moderate +F2BF3O_MOD.sequoia <- BootPed( +F2BF3O_MOD.GenoM_final, +F2BF3O_MOD.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +25, # n.SNP +5)# n.bootstrap +# moderate +F2BF3O_MOD.sequoia <- BootPed( +F2BF3O_MOD.GenoM_final, +F2BF3O_MOD.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +5)# n.bootstrap +# moderate +F2BF3O_MOD.sequoia <- BootPed( +F2BF3O_MOD.GenoM_final, +F2BF3O_MOD.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +5)# n.bootstrap +# moderate +F2BF3O_MOD.sequoia <- BootPed( +F2BF3O_MOD.GenoM_final, +F2BF3O_MOD.LifeHistData,# LifeHistData +-2,# threshold.LLR +0.1,# min.LLR +50, # n.SNP +5)# n.bootstrap diff --git a/RAnalysis/Scripts/Popgen/05_parentage.Rmd b/RAnalysis/Scripts/Popgen/05_parentage.Rmd index a2d25b6..9804652 100644 --- a/RAnalysis/Scripts/Popgen/05_parentage.Rmd +++ b/RAnalysis/Scripts/Popgen/05_parentage.Rmd @@ -536,7 +536,7 @@ based on the number of SNPs and number of iterations. ```{r pedigree function 'BootPed'} -BootPed <- function(GenoM, LifeHistData, n.SNPs, n.bootstrap) { +BootPed <- function(GenoM, LifeHistData, threshold.LLR, min.LLR, n.SNPs, n.bootstrap) { pedigree.raw <- data.frame(matrix(ncol = 7)) # create dataframe to save cumunalitively during for loop colnames(pedigree.raw) <- c('id', 'loci', 'n.loci', 'dam', 'LLRdam', 'sire' , 'LLRsire') # names for comuns in the for loop @@ -548,6 +548,7 @@ BootPed <- function(GenoM, LifeHistData, n.SNPs, n.bootstrap) { bootstrap_GenoM <- GenoM[,randomSNPs] SNP_IDs <- paste0(colnames(bootstrap_GenoM), collapse=", ") # run seqquoia for these 50 SNPs + bootstrap.sequoia <- sequoia( GenoM = bootstrap_GenoM, LifeHistData = LifeHistData, @@ -556,9 +557,9 @@ BootPed <- function(GenoM, LifeHistData, n.SNPs, n.bootstrap) { Module = "par", Err = 1e-04, # threshold for log 10 likelihood (LLR) between a proposed relationship versus unrelated - Tfilter = -2, + Tfilter = threshold.LLR, # minimum LLR required for acceptance of proposed relationship relative to the next most likely - Tassign = 0.2, + Tassign = min.LLR, # MaxSibshipSize = 100, # breeding system complexity Complex = "full", @@ -604,7 +605,7 @@ pedigree.summary <- pedigree.merge.dams.sires %>% # truncates the 'raw' file above for just a single id per row, as that with the GREATEST absLLRtotal pedigree.final <- pedigree.summary %>% group_by(id) %>% # group by id - there are mutliple occurances of pedigree assign for ids! - dplyr::top_n(1, absLLRtotal) %>% # take the best hit essentiall, highest LLR total + dplyr::top_n(1, absLLRtotal) %>% # take the best hit essential, highest LLR total dplyr::select(id, dam, sire, absLLRtotal, Fertilization_pairs, loci) # select the target columns for next sequoia run # Now we want to ensure we have pedigree results for all target offspring @@ -626,7 +627,7 @@ if ( length(unique(pedigree.raw$id)) < offspring.n) {# TRUE = not enough bootstr nrow(pedigree.final %>% filter(Fertilization_pairs %in% 'TRUE')), "/", nrow(pedigree.final), - "true fertilized pairs")) + " true fertilized pairs")) print(paste0("Edit the n.SNP or n.bootstrap call for better pedigree")) @@ -654,6 +655,13 @@ return(list(raw = pedigree.raw, ## BootPed Run +**Important notes about inputs** + +* **n.SNP:** It's important to note that while fewer loci may result in more pedigree assignments, the accuracy of these assignments may be compromised. For optimal results, it is recommended to use a sufficient number of high-quality, independent SNPs to balance between assignment rate and accuracy + - (a) Reduced stringency: With fewer loci, the likelihood thresholds for assigning relationships become less stringent, allowing more potential relationships to be considered + - (b) Increased false positives: A smaller number of loci may lead to an increased rate of false-positive assignments, as there is less genetic information to discriminate between true and spurious relationships + - (c) Overestimation of relatedness: Fewer loci can result in overestimation of relatedness, particularly for more distant relationships like cousins or grandparent-grandchild pairs + * first sanity check that the matrix == life history data, edit if necessary ```{r BootPed prep} @@ -710,6 +718,8 @@ sort(rownames(F1BF2O_MOD.GenoM_final)) == sort(F1BF2O_MOD.LifeHistData$ID) inputs: - GenoM = F0BF1O.GenoM_final - LifeHistData = F0BF1O.LifeHistData + - threshold.LLR = -2 + - min.LLR = 0.1 - n.SNP = 50 - n.bootstrap = 100 @@ -721,10 +731,12 @@ output: # F0 AND F1 OFFPSRING ALL :::::::::::::::::::::::::::::::::::::::::: F0BF1O.sequoia <- BootPed( - F0BF1O.GenoM_final, - F0BF1O.LifeHistData, - 50, - 100 + F0BF1O.GenoM_final, # GenoM + F0BF1O.LifeHistData, # LifeHistData + -2,# threshold.LLR + 0.1,# min.LLR + 50, # n.SNP + 100 # n.bootstrap ) # data files @@ -768,10 +780,12 @@ output: # low F1BF2O_LOW.sequoia <- BootPed( - F1BF2O_LOW.GenoM_final, - F1BF2O_LOW.LifeHistData, - 30, - 100 + F1BF2O_LOW.GenoM_final, # GenoM + F1BF2O_LOW.LifeHistData, # LifeHistData + -2,# threshold.LLR + 0.1,# min.LLR + 50, # n.SNP + 50 # n.bootstrap ) # run on 12/21/2024 # [1] "Predigree summary: 27/46 offspring with pedigree assignment; 12/27true fertilized pairs" @@ -800,10 +814,11 @@ write.csv(F1BF2O_LOW.pedigree.final, "RAnalysis/Output/Popgen/parentage_sequoia/ # moderate F1BF2O_MOD.sequoia <- BootPed( F1BF2O_MOD.GenoM_final, - F1BF2O_MOD.LifeHistData, - 30, - 100 - ) + F1BF2O_MOD.LifeHistData, # LifeHistData + -2,# threshold.LLR + 0.1,# min.LLR + 50, # n.SNP + 1)# n.bootstrap # run on 12/21/2024 # [1] "Predigree summary: 46/46 offspring with pedigree assignment; 21/47 TRUE fertilized pairs" # [1] "Proceed, all offspring have at least 1 pedigree assignment" @@ -832,10 +847,11 @@ write.csv(F1BF2O_MOD.pedigree.final, "RAnalysis/Output/Popgen/parentage_sequoia/ # high F1BF2O_HIGH.sequoia <- BootPed( F1BF2O_HIGH.GenoM_final, - F1BF2O_HIGH.LifeHistData, - 30, - 200 - ) + F1BF2O_HIGH.LifeHistData,# LifeHistData + -2,# threshold.LLR + 0.1,# min.LLR + 50, # n.SNP + 1)# n.bootstrap # run on 12/21/2024 # [1] "Predigree summary: 33/33 offspring with pedigree assignment; 18/33 TRUE fertilized pairs" # [1] "Proceed, all offspring have at least 1 pedigree assignment" @@ -879,26 +895,97 @@ output: # F2 AND F3 OFFPSRING (treatment specific) ::::::::::::::::::::::::: # low -BootPed( - F2BF3O_LOW.GenoM_final, - F2BF3O_LOW.LifeHistData, - 50, - 100 - ) +F2BF3O_LOW.sequoia <- BootPed( + F2BF3O_LOW.GenoM_final, + F2BF3O_LOW.LifeHistData,# LifeHistData + -2,# threshold.LLR + 0.1,# min.LLR + 50, # n.SNP + 10)# n.bootstrap +# run on 12/21/2024 +# [1] "Predigree summary: 27/46 offspring with pedigree assignment; 12/27true fertilized pairs" +# [1] "Edit the n.SNP or n.bootstrap call for better pedigree" + +# data files +F2BF3O_LOW.pedigree.raw <- as.data.frame(F2BF3O_LOW.sequoia$raw) +rownames(F2BF3O_LOW.pedigree.raw) <- NULL # renomve rownames + +F2BF3O_LOW.pedigree.summary <- as.data.frame(F2BF3O_LOW.sequoia$summary) +rownames(F2BF3O_LOW.pedigree.summary) <- NULL # renomve rownames +F2BF3O_LOW.pedigree.summary <- apply(F2BF3O_LOW.pedigree.summary,2,as.character) # need to do this to output the list of batches! + +F2BF3O_LOW.pedigree.final <- as.data.frame(F2BF3O_LOW.sequoia$final) +rownames(F2BF3O_LOW.pedigree.final) <- NULL # renomve rownames + +# output +write.csv(F2BF3O_LOW.pedigree.raw, "RAnalysis/Output/Popgen/parentage_sequoia/F2BF3O/Low/F2BF3O_Low_raw.csv") +write.csv(F2BF3O_LOW.pedigree.summary, "RAnalysis/Output/Popgen/parentage_sequoia/F2BF3O/Low/F2BF3O_Low_summary.csv") +write.csv(F2BF3O_LOW.pedigree.final, "RAnalysis/Output/Popgen/parentage_sequoia/F2BF3O/Low/F2BF3O_Low_final.csv") + + + + + # moderate -BootPed( - F2BF3O_MOD.GenoM_final, - F2BF3O_MOD.LifeHistData, - 50, - 100 - ) +F2BF3O_MOD.sequoia <- BootPed( + F2BF3O_MOD.GenoM_final, + F2BF3O_MOD.LifeHistData,# LifeHistData + -2,# threshold.LLR + 0.1,# min.LLR + 50, # n.SNP + 5)# n.bootstrap +# run on 12/21/2024 +# [1] "Predigree summary: 46/46 offspring with pedigree assignment; 21/47 TRUE fertilized pairs" +# [1] "Proceed, all offspring have at least 1 pedigree assignment" + +# data files +F2BF3O_MOD.pedigree.raw <- as.data.frame(F2BF3O_MOD.sequoia$raw) +rownames(F2BF3O_MOD.pedigree.raw) <- NULL # renomve rownames + +F2BF3O_MOD.pedigree.summary <- as.data.frame(F2BF3O_MOD.sequoia$summary) +rownames(F2BF3O_MOD.pedigree.summary) <- NULL # renomve rownames +F2BF3O_MOD.pedigree.summary <- apply(F2BF3O_MOD.pedigree.summary,2,as.character) # need to do this to output the list of batches! + +F2BF3O_MOD.pedigree.final <- as.data.frame(F2BF3O_MOD.sequoia$final) +rownames(F2BF3O_MOD.pedigree.final) <- NULL # renomve rownames + +# output +write.csv(F2BF3O_MOD.pedigree.raw, "RAnalysis/Output/Popgen/parentage_sequoia/F2BF3O/Moderate/F2BF3O_Mod_raw.csv") +write.csv(F2BF3O_MOD.pedigree.summary, "RAnalysis/Output/Popgen/parentage_sequoia/F2BF3O/Moderate/F2BF3O_Mod_summary.csv") +write.csv(F2BF3O_MOD.pedigree.final, "RAnalysis/Output/Popgen/parentage_sequoia/F2BF3O/Moderate/F2BF3O_Mod_final.csv") + + + + + + # high -BootPed( - F2BF3O_HIGH.GenoM_final, - F2BF3O_HIGH.LifeHistData, - 50, - 100 - ) +F2BF3O_HIGH.sequoia <- BootPed( + F2BF3O_HIGH.GenoM_final, + F2BF3O_HIGH.LifeHistData,# LifeHistData + -2,# threshold.LLR + 0.1,# min.LLR + 50, # n.SNP + 50)# n.bootstrap +# run on 12/21/2024 +# [1] "Predigree summary: 33/33 offspring with pedigree assignment; 18/33 TRUE fertilized pairs" +# [1] "Proceed, all offspring have at least 1 pedigree assignment" + +# data files +F2BF3O_HIGH.pedigree.raw <- as.data.frame(F2BF3O_HIGH.sequoia$raw) +rownames(F2BF3O_HIGH.pedigree.raw) <- NULL # renomve rownames + +F2BF3O_HIGH.pedigree.summary <- as.data.frame(F2BF3O_HIGH.sequoia$summary) +rownames(F2BF3O_HIGH.pedigree.summary) <- NULL # renomve rownames +F2BF3O_HIGH.pedigree.summary <- apply(F2BF3O_HIGH.pedigree.summary,2,as.character) # need to do this to output the list of batches! + +F2BF3O_HIGH.pedigree.final <- as.data.frame(F2BF3O_HIGH.sequoia$final) +rownames(F2BF3O_HIGH.pedigree.final) <- NULL # renomve rownames + +# output +write.csv(F2BF3O_HIGH.pedigree.raw, "RAnalysis/Output/Popgen/parentage_sequoia/F2BF3O/High/F2BF3O_High_raw.csv") +write.csv(F2BF3O_HIGH.pedigree.summary, "RAnalysis/Output/Popgen/parentage_sequoia/F2BF3O/High/F2BF3O_High_summary.csv") +write.csv(F2BF3O_HIGH.pedigree.final, "RAnalysis/Output/Popgen/parentage_sequoia/F2BF3O/High/F2BF3O_High_final.csv") ``` ## `CalcPairLL`