## ----eval = FALSE------------------------------------------------------------- # install.packages("polymapR") # library(polymapR) ## ----eval = FALSE------------------------------------------------------------- # geno <- read.csv("fitPoly_4x_output_2343_SNPs.csv") ## ----echo = FALSE------------------------------------------------------------- library(polymapR) data("gp_vignette_data") individuals <- gp_vignette_data$individuals parent1 <- gp_vignette_data$parent1 parent2 <- gp_vignette_data$parent2 geno.sub <- gp_vignette_data$geno.sub chk1.sub <- gp_vignette_data$chk1.sub pardose.sub <- gp_vignette_data$pardose.sub maxP.chk <- gp_vignette_data$maxP.chk SN_SN_P1.sub <- gp_vignette_data$SN_SN_P1.sub P1_homologues <- gp_vignette_data$P1_homologues P2_homologues <- gp_vignette_data$P2_homologues LGHomDf_P1 <- gp_vignette_data$LGHomDf_P1 LGHomDf_P2 <- gp_vignette_data$LGHomDf_P2 LGHomDf_P1a <- gp_vignette_data$LGHomDf_P1a P1_SxS_Assigned <- gp_vignette_data$P1_SxS_Assigned P2_SxS_Assigned <- gp_vignette_data$P2_SxS_Assigned ## ----eval = FALSE------------------------------------------------------------- # knitr::kable(head(geno)) ## ----echo = FALSE------------------------------------------------------------- knitr::kable(head(geno.sub)) ## ----eval = FALSE------------------------------------------------------------- # parent1.reps <- c("P1","P1a") # parent2.reps <- "P2" # individuals <- setdiff(unique(geno$SampleName),c(parent1.reps,parent2.reps)) ## ----eval = FALSE------------------------------------------------------------- # chk1 <- checkF1(input_type = "probabilistic", # probgeno_df = geno, # parent1 = parent1.reps, # parent2 = parent2.reps, # F1 = individuals, # polysomic = TRUE, # disomic = FALSE, # mixed = FALSE, # ploidy = 4) ## ----eval = FALSE------------------------------------------------------------- # chk1$checked_F1[1:5,] ## ----echo = FALSE------------------------------------------------------------- chk1.sub[1:5,] ## ----eval = FALSE------------------------------------------------------------- # pardose <- polymapR:::assign_parental_dosage(chk = chk1, # probgeno_df = geno) # knitr::kable(head(pardose)) ## ----echo = FALSE------------------------------------------------------------- knitr::kable(head(pardose.sub)) ## ----eval = FALSE------------------------------------------------------------- # length(which(chk1$checked_F1$qall_mult == 0)) #859 markers ## ----eval = FALSE------------------------------------------------------------- # remove.mark <- chk1$checked_F1[chk1$checked_F1$qall_mult==0,"MarkerName"] # geno1 <- geno[!geno$MarkerName %in% remove.mark,] ## ----eval = FALSE------------------------------------------------------------- # gpo <- gp_overview(probgeno_df = geno1) ## ----out.width = "500px", echo = FALSE, fig.align="center"-------------------- knitr::include_graphics("figures/gp_overview.png") ## ----eval = FALSE------------------------------------------------------------- # geno1 <- gpo$probgeno_df ## ----eval = FALSE------------------------------------------------------------- # maxP.chk <- check_maxP(probgeno_df = geno1) ## ----echo = FALSE------------------------------------------------------------- print(maxP.chk) ## ----out.width = "550px", echo = FALSE, fig.align="center"-------------------- knitr::include_graphics("figures/maxPdist.png") ## ----eval = FALSE------------------------------------------------------------- # geno2 <- screen_for_duplicate_individuals.gp(probgeno_df = geno1, # ploidy = 4, # parent1 = parent1.reps, # parent2 = parent2.reps, # F1 = individuals) ## ----out.width = "550px", echo = FALSE, fig.align="center"-------------------- knitr::include_graphics("figures/dup_indivs.png") ## ----echo = FALSE------------------------------------------------------------- write("\nNo duplicates found\n",stdout()) ## ----------------------------------------------------------------------------- nc <- parallel::detectCores() - 2 ## ----eval = FALSE------------------------------------------------------------- # chk1 <- checkF1(input_type = "probabilistic", # probgeno_df = geno2, # parent1 = parent1.reps, # parent2 = parent2.reps, # F1 = individuals, # polysomic = TRUE, # disomic = FALSE, # mixed = FALSE, # ploidy = 4) ## ----eval = FALSE------------------------------------------------------------- # SN_SN_P1 <- linkage.gp(probgeno_df = geno2, # chk = chk1, # markertype1 = c(1,0), # target_parent = "P1", # LOD_threshold = 3, # ncores = nc) ## ----eval = FALSE------------------------------------------------------------- # head(SN_SN_P1) ## ----echo = FALSE------------------------------------------------------------- head(SN_SN_P1.sub) ## ----eval = FALSE------------------------------------------------------------- # par(mfrow = c(2,1)) # P1_homologues <- cluster_SN_markers(linkage_df = SN_SN_P1, # LOD_sequence = seq(3,12,0.5), # LG_number = 5, # ploidy = 4, # parentname = "P1", # plot_clust_size = FALSE, # min_clust_size = 3) ## ----out.width = "600px", echo = FALSE, fig.align="center"-------------------- knitr::include_graphics("figures/p1_clustering_gp.png") ## ----------------------------------------------------------------------------- sort(table(P1_homologues[['7']]$cluster),decreasing = T) length(table(P1_homologues[['7']]$cluster)) ## ----eval = FALSE------------------------------------------------------------- # SN_SN_P2 <- linkage.gp(probgeno_df = geno2, # chk = chk1, # markertype1 = c(1,0), # target_parent = "P2", # LOD_threshold = 3, # ncores = nc) ## ----eval = FALSE------------------------------------------------------------- # P2_homologues <- cluster_SN_markers(linkage_df = SN_SN_P2, # LOD_sequence = seq(3,12,0.5), # LG_number = 5, # ploidy = 4, # parentname = "P2", # plot_clust_size = F, # min_clust_size = 3) ## ----out.width = "600px", echo = FALSE, fig.align="center"-------------------- knitr::include_graphics("figures/p2_clustering_gp.png") ## ----------------------------------------------------------------------------- sort(table(P2_homologues[['6']]$cluster),decreasing = T) length(table(P2_homologues[['6']]$cluster)) ## ----eval = FALSE------------------------------------------------------------- # SN_SS_P1 <- linkage.gp(probgeno_df = geno2, # chk = chk1, # markertype1 = c(1,0), # markertype2 = c(1,1), # target_parent = "P1", # ncores = nc) # # SN_SS_P2 <- linkage.gp(probgeno_df = geno2, # chk = chk1, # markertype1 = c(1,0), # markertype2 = c(1,1), # target_parent = "P2", # ncores = nc) ## ----eval = FALSE------------------------------------------------------------- # LGHomDf_P1 <- bridgeHomologues(cluster_stack = P1_homologues[["7"]], # cluster_stack2 = P2_homologues[["7"]], # linkage_df = SN_SS_P1, # linkage_df2 = SN_SS_P2, # LOD_threshold = 5, # LG_number = 5) ## ----out.width = "450px", echo = FALSE, fig.align="center"-------------------- knitr::include_graphics("figures/P1_bridges.png") ## ----------------------------------------------------------------------------- table(LGHomDf_P1$LG,LGHomDf_P1$homologue) ## ----eval = FALSE------------------------------------------------------------- # LGHomDf_P2 <- bridgeHomologues(cluster_stack = P2_homologues[["6"]], # cluster_stack2 = P1_homologues[["6"]], # linkage_df = SN_SS_P2, # linkage_df2 = SN_SS_P1, # LOD_threshold = 5, # LG_number = 5) ## ----out.width = "450px", echo = FALSE, fig.align="center"-------------------- knitr::include_graphics("figures/P2_bridges.png") ## ----------------------------------------------------------------------------- table(LGHomDf_P2$LG,LGHomDf_P2$homologue) ## ----eval = FALSE------------------------------------------------------------- # LGHomDf_P1a <-cluster_per_LG(LG = 3, # linkage_df = SN_SN_P1[SN_SN_P1$phase == "coupling",], # LG_hom_stack = LGHomDf_P1, # LOD_sequence = 3:10, # modify_LG_hom_stack = TRUE, # network.layout = "stacked", # nclust_out = 4, # label.offset=1.2) ## ----out.width = "450px", echo = FALSE, fig.align="center"-------------------- knitr::include_graphics("figures/cluster_per_LG_P1.3.png") ## ----------------------------------------------------------------------------- table(LGHomDf_P1a$LG,LGHomDf_P1a$homologue) ## ----------------------------------------------------------------------------- head(LGHomDf_P1a) ## ----eval = FALSE------------------------------------------------------------- # P1_SxS_Assigned <- assign_linkage_group(linkage_df = SN_SS_P1, # LG_hom_stack = LGHomDf_P1a, # SN_colname = "marker_a", # unassigned_marker_name = "marker_b", # phase_considered = "coupling", # LG_number = 5, # LOD_threshold = 3, # ploidy = 4) # # P2_SxS_Assigned <- assign_linkage_group(linkage_df = SN_SS_P2, # LG_hom_stack = LGHomDf_P2, # SN_colname = "marker_a", # unassigned_marker_name = "marker_b", # phase_considered = "coupling", # LG_number = 5, # LOD_threshold = 3, # ploidy = 4) ## ----------------------------------------------------------------------------- LGHomDf_P2c <- consensus_LG_names(modify_LG = LGHomDf_P2, template_SxS = P1_SxS_Assigned, modify_SxS = P2_SxS_Assigned) ## ----eval = FALSE------------------------------------------------------------- # save(LGHomDf_P1a,LGHomDf_P2c, file = "LGHomDf_stacks.Rdata") #for example.. ## ----eval = FALSE------------------------------------------------------------- # P2_SxS_Assigned <- assign_linkage_group(linkage_df = SN_SS_P2, # LG_hom_stack = LGHomDf_P2c, #this is changed # SN_colname = "marker_a", # unassigned_marker_name = "marker_b", # phase_considered = "coupling", # LG_number = 5, # LOD_threshold = 3, # ploidy = 4) ## ----eval = FALSE------------------------------------------------------------- # marker_assignments_P1 <- homologue_lg_assignment(input_type = "probabilistic", # probgeno_df = geno2, # chk = chk1, # assigned_list = list(P1_SxS_Assigned), # assigned_markertypes = list(c(1,1)), # LG_hom_stack = LGHomDf_P1a, # target_parent = "P1", # other_parent = "P2", # ploidy = 4, # pairing = "random", # convert_palindrome_markers = FALSE, # LG_number = 5, # LOD_threshold = 3, # write_intermediate_files = FALSE) # # marker_assignments_P2 <- homologue_lg_assignment(input_type = "probabilistic", # probgeno_df = geno2, # chk = chk1, # assigned_list = list(P2_SxS_Assigned), # assigned_markertypes = list(c(1,1)), # LG_hom_stack = LGHomDf_P2c, # target_parent = "P2", # other_parent = "P1", # ploidy = 4, # pairing = "random", # convert_palindrome_markers = FALSE, # LG_number = 5, # LOD_threshold = 3, # write_intermediate_files = FALSE) ## ----eval = FALSE------------------------------------------------------------- # marker_assignments <- check_marker_assignment(marker_assignments_P1,marker_assignments_P2) ## ----eval = FALSE------------------------------------------------------------- # saveRDS(marker_assignments, file = "marker_assignments.RDS") ## ----eval = FALSE------------------------------------------------------------- # all_linkages_list_P1 <- finish_linkage_analysis(input_type = "probabilistic", # marker_assignment = marker_assignments$P1, # probgeno_df = geno2, # chk = chk1, # target_parent = "P1", # other_parent = "P2", # convert_palindrome_markers = FALSE, # ploidy = 4, # pairing = "random", # LG_number = 5, # ncores = nc) # # all_linkages_list_P2 <- finish_linkage_analysis(input_type = "probabilistic", # marker_assignment = marker_assignments$P2, # probgeno_df = geno2, # chk = chk1, # target_parent = "P2", # other_parent = "P1", # convert_palindrome_markers = TRUE, # ploidy = 4, # pairing = "random", # LG_number = 5, # ncores = nc) ## ----eval = FALSE------------------------------------------------------------- # linkages <- list() # for(lg in names(all_linkages_list_P1)){ # linkages[[lg]] <- rbind(all_linkages_list_P1[[lg]], all_linkages_list_P2[[lg]]) # } ## ----eval = FALSE------------------------------------------------------------- # saveRDS(linkages, file = "linkages.RDS") ## ----eval = FALSE------------------------------------------------------------- # linkages <- readRDS("linkages.RDS") ## ----eval = FALSE------------------------------------------------------------- # integrated.maplist <- MDSMap_from_list(linkages) ## ----eval = FALSE------------------------------------------------------------- # phased.maplist <- create_phased_maplist(input_type = "probabilistic", # maplist = integrated.maplist, # chk = chk1, # ploidy = 4, # marker_assignment.1 = marker_assignments$P1, # marker_assignment.2 = marker_assignments$P2)