4_GOCS.R

# Long-term selection in layers (GOCS)
# Ivan Pocrnic

# Clean
rm(list = ls())
# Load packages
library("AlphaSimR")
library("tidyverse")
# Load fje
source("../functions.R")
# Set WD
args = commandArgs(trailingOnly=TRUE)
if (length(args) < 1) {
  stop("Must provide replicate number [1], number of sires [2] and OCS degrees [3] !!!")
}
rep = args[1]
noSires = as.numeric(args[2])
degrees = as.numeric(args[3])
scenarioMain = paste("Additive", noSires, sep="_")
scenarioName = paste("Rep_", rep, "_Additive_OCSNOMP_", degrees, sep="")
dire = paste("../../data",scenarioMain,scenarioName, sep = "/")
unlink(dire, recursive = TRUE)
dir.create(path = dire, recursive = TRUE, showWarnings = FALSE)
setwd(dir = dire)

noDams=1080

# Create blupf90 parameter files and bash scripts
prepare_par()

# For the sake of comparison, use the same Burn-in as for the Standard scenario

bdirname = paste("Rep_", rep, "_Additive", sep="")
bdir = paste("..",bdirname, "burnin.RData", sep = "/")
load(bdir)


#### Selection: 10 generations of genomic selection (10 x 2 years) ####

# Parents in the 1st gen of selection are from the last generation of burnin:
GSStartMales   = StartMalesBurnin
GSStartFemales = StartFemalesBurnin

# Year is the last year of the burnin:
year = year_burnin

# Separately keep selection candidates dams to use for the 1st gen of OCS
# For this scenario (OCS for both males and females) these dams are not selected, but are from the last gen. burnin selection candidates (p3f)
# Kind of bizarre situation: Those selected in burn-in, are also starting females for this scenario  
# This has to be done when switching from BLUP burn-in to GS-OCS scenarios 

damg = vector("list",30)
damg[[6]] = p3f



# Only for scenarios using burin-in population created by older version of AlphaSimR

# if(simParam$version < ){
#}

genMap_new = vector("list",founderPop@nChr)

for(i in 1:founderPop@nChr){
  genMap_new[[i]] = founderPop@genMap[[1]]
  names(genMap_new[[i]]) = paste(i,1:length(genMap_new[[i]]),sep="_")
}
founderPop@genMap = genMap_new
str(founderPop@genMap)

SP$switchFemaleMap(genMap_new)
SP$switchMaleMap(genMap_new)

GSStartMales@iid = as.integer(GSStartMales@id)
GSStartFemales@iid = as.integer(GSStartFemales@id)

str(SP$femaleMap)
str(SP$maleMap)


for(gen in 6:15){
    Program = paste("OCSNOMP_", degrees, sep="")
    cat("Working on the round:",Program, ":", gen,"\n")

    # Q1:
    year = year + 1

    # Mate initial sires and dams:
    p_start = selectCross(pop=c(GSStartMales,GSStartFemales), nFemale = noDams, nMale = noSires, nCrosses = noDams, nProgeny = 20)

    # Selected first 4 M & 9 F hatched per mating (per full-sib family):
    p2m = selectWithinFam(p_start,nInd = 4, famType = "B", sex = "M", use = "rand")
    p2f = selectWithinFam(p_start,nInd = 9, famType = "B", sex = "F", use = "rand")

    # Prepare genotypes I:
    # Each time select all relevant genotypes (ugly code now!)
    # Also need if/else due to transition phase from burnin to genomics
    if(gen > 6){
      # Selected dams from the G-4 keep genotypes:
      run_prepareGENO(p1fs)
      # Selected sires from the G-4 keep genotypes:
      run_prepareGENO(p2ms)
      # Selected dams from the G-3 keep genotypes:
      run_prepareGENO(p2fs)
      # Selected sires from the G-3 keep genotypes:
      run_prepareGENO(p3ms)
      # Selected dams from the G-2 keep genotypes:
      run_prepareGENO(p3fs)
      # Selected sires from the G-2 keep genotypes:
      run_prepareGENO(p4ms)
      # Selected dams from the G-1 keep genotypes:
      run_prepareGENO(GSStartFemales)
      # Selected sires from the G-1 keep genotypes:
      run_prepareGENO(GSStartMales)
      # All male candidates from this generation get genotypes:
      run_prepareGENO(p2m)
    } else {
      # Selected dams from the G-1 keep genotypes:
      run_prepareGENO(GSStartFemales)
      # Selected sires from the G-1 keep genotypes:
      run_prepareGENO(GSStartMales)
      # All male candidates from this generation get genotypes:
      run_prepareGENO(p2m)
    }


    #### First batch of selection candidates ####
    candidatesgroup = 1
    # Generate records:
    RecSys = RecSysMale(RecSys, p2m)
    RecSys = RecSysFemale(RecSys, p2f)

    # Event @ 25 weeks (6m) - p2f get T1 phenotypes
    # Event @ 52 weeks (12m) - damg[[gen]] / p5f get T2 phenotypes (From Previous Generation)
    # Event @ 100 weeks (23m) - p3f get T3 phenotypes (From Previous Generation)

    # Q2:

    OldDir = getwd()
    Dir = paste("GBlup", Program, gen, candidatesgroup, sep = "_")
    unlink(paste(OldDir,Dir, sep = "/"), recursive = TRUE)
    dir.create(path = Dir, showWarnings = FALSE)
    setwd(dir = Dir)
    # getwd()

    # Prepare pedigree and datafile for blupf90
    # Removes phenotypes that are not really available at given timepoint
    run_prepare(RecSys)

    # Run blupf90, and update GEBV in recording system
    RecSys = run_gblup(RecSys)

    setwd(dir = OldDir)

    # Clean marker files from the last round of selection:
    unlink("mrk*")
    unlink("markeri*")

    # Set EBV for selection candidates:
    p2m@ebv = as.matrix(RecSys[RecSys$IId %in% p2m@id, c("EbvT1", "EbvT2","EbvT3")])
    p2f@ebv = as.matrix(RecSys[RecSys$IId %in% p2f@id, c("EbvT1", "EbvT2","EbvT3")])

    # Summarize Selection Candidates from Batch 1
    CSumm = PullSumm(CSumm,p2m,"M")
    CSumm = PullSumm(CSumm,p2f,"F")
    CSummTogether = PullSummTogether(CSummTogether, c(p2m,p2f))

    # WARNING: For the females (p2f) from batch 1 we get EBV on available T1 and T2
    # In reality they don't have T2 at this time point
    # Why doesn't matter: We don't use them till later when they actually get T2.

    # Select next generations F & M based on OCS
    
    # Select populations to be included in OCS
    candidates_plan=c(damg[[gen]],p2m)
    
    # Run OCS:
    run_OCS_NOMP(candidates_plan, RecSys, degrees, noSires, noDams)
    system(command = paste("cp AlphaMate.log AlphaMate.log",gen,year,candidatesgroup,sep = "_"))
    
    # Read crossPlan
    crossPlan = read.table(file="ContributorsModeOptTarget1.txt", header=TRUE)

    # Population of selected males based on ocs:
    p2ms = p2m[p2m@id %in% filter(crossPlan, Gender == 1)[[1]],]
    #a1=(filter(crossPlan, Gender == 1) %>% dplyr::select(Id))
    
    # Population of selected females based on ocs:
    p1fs_tmp=damg[[gen]]
    p1fs = p1fs_tmp[p1fs_tmp@id %in% filter(crossPlan, Gender == 2)[[1]],]
  
    p33 = selectCross(pop=c(p2ms,p1fs), nFemale = noDams, nMale = noSires, nCrosses = noDams, nProgeny = 20)
  
    # Selected first 4 M & 9 F hatched per mating (per full-sib family):
    p3m = selectWithinFam(p33,nInd = 4, famType = "B", sex = "M", use = "rand")
    p3f = selectWithinFam(p33,nInd = 9, famType = "B", sex = "F", use = "rand")

    # Prepare genotypes II:
    # Each time select all relevant genotypes (ugly code now!)
    # Also need if/else due to transition phase from burnin to genomics
    if(gen > 6){
      # Selected dams from the G-4 keep genotypes:
      run_prepareGENO(p2fs)
      # Selected sires from the G-4 keep genotypes:
      run_prepareGENO(p3ms)
      # Selected dams from the G-3 keep genotypes:
      run_prepareGENO(p3fs)
      # Selected sires from the G-3 keep genotypes:
      run_prepareGENO(p4ms)
      # Selected dams from the G-2 keep genotypes:
      run_prepareGENO(GSStartFemales)
      # Selected sires from the G-2 keep genotypes:
      run_prepareGENO(GSStartMales)
      # Selected dams from the G-1 keep genotypes:
      run_prepareGENO(p1fs)
      # Selected sires from the G-1 keep genotypes:
      run_prepareGENO(p2ms)
      # All male candidates from this generation get genotypes:
      run_prepareGENO(p3m)
    } else {
      # Selected dams from the G-2 keep genotypes:
      run_prepareGENO(GSStartFemales)
      # Selected sires from the G-2 keep genotypes:
      run_prepareGENO(GSStartMales)
      # Selected dams from the G-1 keep genotypes:
      # Some of those can be equal to initial dams so remove duplicates:
      tmp_gd = p1fs[!(p1fs@id) %in% GSStartFemales@id,]
      run_prepareGENO(tmp_gd)
      # Selected sires from the G-1 keep genotypes:
      run_prepareGENO(p2ms)
      # All male candidates from this generation get genotypes:
      run_prepareGENO(p3m)
    }


    ##### Second batch of selection candidates ####
    candidatesgroup = 2
    # Generate records:
    RecSys = RecSysMale(RecSys, p3m)
    RecSys = RecSysFemale(RecSys, p3f)

    # Event @ 25 weeks (6m) - p3f get T1 phenotypes
    # Event @ 52 weeks (12m) - p2f get T2 phenotypes
    # Event @ 100 weeks (23m) - p4f get T3 phenotypes (From Previous Generation)

    OldDir = getwd()
    Dir = paste("GBlup", Program, gen, candidatesgroup, sep = "_")
    unlink(paste(OldDir,Dir, sep = "/"), recursive = TRUE)
    dir.create(path = Dir, showWarnings = FALSE)
    setwd(dir = Dir)
    #getwd()

    # Prepare pedigree and datafile for blupf90
    # Removes phenotypes that are not really available at given timepoint
    run_prepare(RecSys)

    # Run blupf90, and update GEBV in recording system
    RecSys = run_gblup(RecSys)

    setwd(dir = OldDir)

    # Clean marker files from the last round of selection:
    unlink("mrk*")
    unlink("markeri*")

    # Set EBV for selection candidates:
    p3m@ebv = as.matrix(RecSys[RecSys$IId %in% p3m@id, c("EbvT1", "EbvT2","EbvT3")])
    p3f@ebv = as.matrix(RecSys[RecSys$IId %in% p3f@id, c("EbvT1", "EbvT2","EbvT3")])

    # Summarize Selection Candidates from Batch 2
    CSumm = PullSumm(CSumm,p3m,"M")
    CSumm = PullSumm(CSumm,p3f,"F")
    CSummTogether = PullSummTogether(CSummTogether, c(p3m,p3f))

    # WARNING: For the females (p3f) from batch 2 we get EBV on available T1 and T2
    # In reality they don't have T2 at this time point
    # Why doesn't matter: We don't use them till later when they actually get T2.

    # Select next generations F & M based on OCS
    
    # Select populations to be included in OCS
    candidates_plan=c(p2f,p3m)
    
    # Run OCS:
    run_OCS_NOMP(candidates_plan, RecSys, degrees, noSires, noDams)
    system(command = paste("cp AlphaMate.log AlphaMate.log",gen,year,candidatesgroup,sep = "_"))
    
    # Read crossPlan
    crossPlan = read.table(file="ContributorsModeOptTarget1.txt", header=TRUE)
    
    # Population of selected males based on ocs:
    p3ms = p3m[p3m@id %in% filter(crossPlan, Gender == 1)[[1]],]

    # Population of selected females based on ocs:
    p2fs = p2f[p2f@id %in% filter(crossPlan, Gender == 2)[[1]],]
    
    
    # Q3:
    year = year + 1
    
    # Mate young OCS selected males with selected females (From Previous Generation)
    # Those females were already selected in previous generation and saved as damg[[gen]]
    p44=selectCross(pop=c(p3ms,p2fs), nFemale = noDams, nMale = noSires, nCrosses = noDams, nProgeny = 20)
    
    # Selected first 4 M & 9 F hatched per mating (per full-sib family):
    p4m = selectWithinFam(p44,nInd = 4, famType = "B", sex = "M", use = "rand")
    p4f = selectWithinFam(p44,nInd = 9, famType = "B", sex = "F", use = "rand")


    # Prepare genotypes III:
    # Each time select all relevant genotypes (ugly code now!)
    # Also need if/else due to transition phase from burnin to genomics
    if(gen > 6){
      # Selected dams from the G-4 keep genotypes:
      run_prepareGENO(p3fs)
      # Selected sires from the G-4 keep genotypes:
      run_prepareGENO(p4ms)
      # Selected dams from the G-3 keep genotypes:
      run_prepareGENO(GSStartFemales)
      # Selected sires from the G-3 keep genotypes:
      run_prepareGENO(GSStartMales)
      # Selected dams from the G-2 keep genotypes:
      run_prepareGENO(p1fs)
      # Selected sires from the G-2 keep genotypes:
      run_prepareGENO(p2ms)
      # Selected dams from the G-1 keep genotypes:
      run_prepareGENO(p2fs)
      # Selected sires from the G-1 keep genotypes:
      run_prepareGENO(p3ms)
      # All male candidates from this generation get genotypes:
      run_prepareGENO(p4m)
    } else {
      # Selected dams from the G-3 keep genotypes:
      run_prepareGENO(GSStartFemales)
      # Selected sires from the G-3 keep genotypes:
      run_prepareGENO(GSStartMales)
      # Selected dams from the G-2 keep genotypes:
      # Some of those can be equal to initial dams so remove duplicates:
      # tmp_gd = p1fs[!(p1fs@id) %in% GSStartFemales,]
      run_prepareGENO(tmp_gd)
      # Selected sires from the G-2 keep genotypes:
      run_prepareGENO(p2ms)
      # Selected dams from the G-1 keep genotypes:
      run_prepareGENO(p2fs)
      # Selected sires from the G-1 keep genotypes:
      run_prepareGENO(p3ms)
      # All male candidates from this generation get genotypes:
      run_prepareGENO(p4m)
    }


    #### Third batch of selection candidates ####
    candidatesgroup = 3
    # Generate records:
    RecSys = RecSysMale(RecSys, p4m)
    RecSys = RecSysFemale(RecSys, p4f)

    # Event @ 25 weeks (6m) - p4f get T1 phenotypes
    # Event @ 52 weeks (12m) - p3f get T2 phenotypes
    # Event @ 100 weeks (23m) - damg[[gen]] / p5f get T3 phenotypes (From Previous Generation)

    OldDir = getwd()
    Dir = paste("GBlup", Program, gen, candidatesgroup, sep = "_")
    unlink(paste(OldDir,Dir, sep = "/"), recursive = TRUE)
    dir.create(path = Dir, showWarnings = FALSE)
    setwd(dir = Dir)
    #getwd()

    # Prepare pedigree and datafile for blupf90
    # Removes phenotypes that are not really available at given timepoint
    run_prepare(RecSys)

    # Run blupf90, and update GEBV in recording system
    RecSys = run_gblup(RecSys)

    setwd(dir = OldDir)

    # Clean marker files from the last round of selection:
    unlink("mrk*")
    unlink("markeri*")

    # Set EBV for selection candidates:
    p4m@ebv = as.matrix(RecSys[RecSys$IId %in% p4m@id, c("EbvT1", "EbvT2","EbvT3")])
    p4f@ebv = as.matrix(RecSys[RecSys$IId %in% p4f@id, c("EbvT1", "EbvT2","EbvT3")])

    # Summarize Selection Candidates from Batch 3
    CSumm = PullSumm(CSumm,p4m,"M")
    CSumm = PullSumm(CSumm,p4f,"F")
    CSummTogether = PullSummTogether(CSummTogether, c(p4m,p4f))

    # WARNING: For the females (p4f) from batch 3 we get EBV on available T1 and T2
    # In reality they dont have T2 at this time point
    # Why doesn't matter: We dont use them till later when they actualy get T2.


    # Select next generations F & M based on OCS
    
    # Select populations to be included in OCS
    candidates_plan=c(p3f,p4m)
    
    # Run OCS:
    run_OCS_NOMP(candidates_plan, RecSys, degrees, noSires, noDams)
    system(command = paste("cp AlphaMate.log AlphaMate.log",gen,year,candidatesgroup,sep = "_"))
    
    # Read crossPlan
    crossPlan = read.table(file="ContributorsModeOptTarget1.txt", header=TRUE)
    
    # Population of selected males based on ocs:
    p4ms = p4m[p4m@id %in% filter(crossPlan, Gender == 1)[[1]],]
    
    # Population of selected females based on ocs:
    p3fs = p3f[p3f@id %in% filter(crossPlan, Gender == 2)[[1]],]
    
    # Q4:
    
    
    # Mate young OCS selected males with selected females (From Previous Generation)

    p55=selectCross(pop=c(p4ms,p3fs), nFemale = noDams, nMale = noSires, nCrosses = noDams, nProgeny = 20)
    

    # Selected first 4 M & 9 F hatched per mating (per full-sib family):
    p5m = selectWithinFam(p55,nInd = 4, famType = "B", sex = "M", use = "rand")
    p5f = selectWithinFam(p55,nInd = 9, famType = "B", sex = "F", use = "rand")


    # Prepare genotypes IV:
    # Each time select all relevant genotypes (ugly code now!)
    # Also need if/else due to transition phase from burnin to genomics
    if(gen > 6){
      # Selected dams from the G-4 keep genotypes:
      run_prepareGENO(GSStartFemales)
      # Selected sires from the G-4 keep genotypes:
      run_prepareGENO(GSStartMales)
      # Selected dams from the G-3 keep genotypes:
      run_prepareGENO(p1fs)
      # Selected sires from the G-3 keep genotypes:
      run_prepareGENO(p2ms)
      # Selected dams from the G-2 keep genotypes:
      run_prepareGENO(p2fs)
      # Selected sires from the G-2 keep genotypes:
      run_prepareGENO(p3ms)
      # Selected dams from the G-1 keep genotypes:
      run_prepareGENO(p3fs)
      # Selected sires from the G-1 keep genotypes:
      run_prepareGENO(p4ms)
      # All male candidates from this generation get genotypes:
      run_prepareGENO(p5m)
    } else {
      # Selected dams from the G-4 keep genotypes:
      run_prepareGENO(GSStartFemales)
      # Selected sires from the G-4 keep genotypes:
      run_prepareGENO(GSStartMales)
      # Selected dams from the G-3 keep genotypes:
      # Some of those can be equal to initial dams so remove duplicates:
      # tmp_gd = p1fs[!(p1fs@id) %in% GSStartFemales,]
      run_prepareGENO(tmp_gd)
      # Selected sires from the G-3 keep genotypes:
      run_prepareGENO(p2ms)
      # Selected dams from the G-2 keep genotypes:
      run_prepareGENO(p2fs)
      # Selected sires from the G-2 keep genotypes:
      run_prepareGENO(p3ms)
      # Selected dams from the G-1 keep genotypes:
      run_prepareGENO(p3fs)
      # Selected sires from the G-1 keep genotypes:
      run_prepareGENO(p4ms)
      # All male candidates from this generation get genotypes:
      run_prepareGENO(p5m)
    }

    #### Forth batch of selection candidates ####
    candidatesgroup = 4
    # Generate records:
    RecSys = RecSysMale(RecSys, p5m)
    RecSys = RecSysFemale(RecSys, p5f)

    # Event @ 25 weeks (6m) - p5f get T1 phenotypes
    # Event @ 52 weeks (12m) - p4f get T2 phenotypes
    # Event @ 100 weeks (23m) - p2f get T3 phenotypes

    OldDir = getwd()
    Dir = paste("GBlup", Program, gen, candidatesgroup, sep = "_")
    unlink(paste(OldDir,Dir, sep = "/"), recursive = TRUE)
    dir.create(path = Dir, showWarnings = FALSE)
    setwd(dir = Dir)
    #getwd()

    # Prepare pedigree and datafile for blupf90
    # Removes phenotypes that are not really available at given timepoint
    run_prepare(RecSys)

    # Run blupf90, and update GEBV in recording system
    RecSys = run_gblup(RecSys)

    setwd(dir = OldDir)

    # Clean marker files from the last round of selection:
    unlink("mrk*")
    unlink("markeri*")

    # Set EBV for selection candidates:
    p5m@ebv = as.matrix(RecSys[RecSys$IId %in% p5m@id, c("EbvT1", "EbvT2","EbvT3")])
    p5f@ebv = as.matrix(RecSys[RecSys$IId %in% p5f@id, c("EbvT1", "EbvT2","EbvT3")])

    # Summarize Selection Candidates from Batch 4
    CSumm = PullSumm(CSumm,p5m,"M")
    CSumm = PullSumm(CSumm,p5f,"F")
    CSummTogether = PullSummTogether(CSummTogether, c(p5m,p5f))

    # WARNING: For the females (p5f) from batch 4 we get EBV on available T1 and T2
    # In reality they don't have T2 at this time point
    # Why it doesn't matter: We don't use them till later when they actually get T2

    # This are initial parents for the next round
    # Select next generations F & M based on OCS
    
    # Select populations to be included in OCS
    candidates_plan=c(p4f,p5m)
    
    # Run OCS:
    run_OCS_NOMP(candidates_plan, RecSys, degrees, noSires, noDams)
    system(command = paste("cp AlphaMate.log AlphaMate.log",gen,year,candidatesgroup,sep = "_"))
    
    # Read crossPlan
    crossPlan = read.table(file="ContributorsModeOptTarget1.txt", header=TRUE)

    # Population of selected males based on ocs:
    GSStartMales = p5m[p5m@id %in% filter(crossPlan, Gender == 1)[[1]],]

    # Population of selected females based on ocs:
    GSStartFemales = p4f[p4f@id %in% filter(crossPlan, Gender == 2)[[1]],]

    # These are actually selected in Q1 next year after they got T2 Phenotype
    # But same as before, we pretend as we already have their EBV
    
    damg[[gen+1]] = p5f

}

save.image("results.RData")