R/OLD/runFSC_step_agg.R

In stranda/holoSimCell: Linked demographic/genetic simulation

#' run fastsimcoal
#'
#' takes a population history and creates sequences based on the coalscent
#'
#' @export
runFSC_step_agg = function(
  ph = ph,				#A new pophist object - (pophist, Nvecs, tmat, struct, hab_suit, coalhist)
  l = landscape, 			#A new landscape object - (details, occupied, empty, sampled, hab_suit, sumrast, samplocsrast, samplocs)
  sample_n = NULL,		#Number of sampled individuals per population
  preLGMparms = NULL,		#This has parms for the refuge, preLGM size and timing
  label = "test",			#Label for FSC simulation files
  delete_files = TRUE,	#Logical - clear out .par, .arp, and other FSC outputs?
  num_cores = 1,			#Number of processors to use for FSC
  exec = "fsc25",			#Executable for FSC (needs to be in $PATH variable)
  loc_parms = NULL,		#Vector of locus parameters
  found_Ne = NULL,			#Founding population size, required for STEP change model		
  gmap = NULL         #gmap relating spatially aggregated cells for coalescent sim to fine-grained cells from forward sim
) {
  
  ############################################################	
  #Error checks
  if(is.null(found_Ne)) {
    stop("Must specify founding Ne if using the STEP change growth model")
  }
  ############################################################
  
  
  ############################################################
  #Create objects used downstream as FSC inputs
  #Initialize objects for FSC
  pops <- ph$pophist
  simhist <- ph$Nvecs
  tmat <- ph$tmat
  struct <- ph$struct
  coalhist <- ph$coalhist
  
  migmat <- vector("list", 1)
  migmat[[1]] <- round(tmat,5)
  diag(migmat[[1]]) <- 0
  
  allevents <- data.frame(time=NA, source = NA, sink = NA, prop_m = NA, new_size = NA, new_growth=NA, new_mig = NA)
  
  ############################################################
  
  
  ############################################################
  #Deal with empty populations (never colonized,
  # outside of suitable area at all times in the sim, 
  # or below threshold suitability)
  
  #Remove them from the migration matrix
  #empty_pops <- sort(unique(c(l$empty, which(rowSums(simhist) == 0))))
  #empty_pops <- l$empty
  #empty_pops <- sort(unique(c(l$empty, pops$pop[!pops$pop %in% coalhist$src & !pops$pop %in% coalhist$snk & !pops$pop %in% l$empty])))
  #empty_pops2 <- sort(unique(c(l$empty, which(rowSums(simhist) == 0))))    #Don't allow migration with populations that went extinct!
  #New Definition of empty populations... JDR 4/1/2020
  empty <- rep(FALSE, nrow(gmap))
  empty[l$empty] <- TRUE
  empty_gpops <- sort(unique(c(
    gmap$gpop[!gmap$gpop %in% gmap$gpop[empty == FALSE]],
    which(rowSums(ph2$Nvecs) == 0),
    gmap$gpop[!gmap$gpop %in% coalhist$src & !gmap$gpop %in% coalhist$snk])))

  if(length(empty_gpops) > 0) {
    simhist <- simhist[-empty_gpops,]
    #migmat[[1]][empty_gpops,] <- 0		#Turns migration off for all empty pops, and any pops that went extinct (e.g., GoM populations)
    #migmat[[1]][,empty_gpops] <- 0
    migmat[[1]] <- migmat[[1]][-empty_gpops,]
    migmat[[1]] <- migmat[[1]][,-empty_gpops]
  } 
  ############################################################
  
  
  ############################################################
  #Build objects with IDs fixed to account for empty populations and the different starting index (0 for FSC, 1 for forward sim)
  #Correct arrival times so you're thinking backwards in time
  #Using the new coalhist and landscape objects
  ngens <- struct["maxtime"]
  oldID <- sort(unique(c(coalhist$src, which(ph$Nvecs[,1] > 0))))
  newID <- c(1:length(oldID))-1
  coalhist$snk <- plyr::mapvalues(coalhist$snk, oldID, newID, warn_missing=FALSE)
  coalhist$src <- plyr::mapvalues(coalhist$src, oldID, newID, warn_missing=FALSE)
  coalhist$time <- ngens - coalhist$time
  ############################################################
  
  
  ##################################################
  # POPULATION INFORMATION
  ##################################################
  #Definitions for pop.info
  #h <- struct["xdim"]*struct["ydim"]   
  h <- max(gmap$gpop)
  
  pop_size <- simhist[,length(simhist[1,])]
  pop_size[pop_size == 0] <- 500			#!# What should we do here?  set to 1?  Set to K?  Set to max attained in this pop?  Can't be 0!!
  
  #Sample size
  sample_size <- rep(0,h)
  #Growth rate = 0
  growth_rate <- rep(0,h)
  
  if(length(empty_gpops > 0)) {
    sample_size <- sample_size[-empty_gpops]
    growth_rate <- growth_rate[-empty_gpops]
  }
  
  ##!## Slight change here... map first to the gpop column of the gmap object, then to the FSC id (old, new) - JDR 4/1/2020
  sample_ids <- plyr::mapvalues(l$sampled, gmap$pop, gmap$gpop, warn_missing=FALSE)
  sample_pops <- plyr::mapvalues(sample_ids, oldID, newID, warn_missing=FALSE)
  sample_size[sample_pops] <- sample_n
  
  pop_info <- fscPopInfo(pop.size = pop_size, sample.size = sample_size, growth.rate = growth_rate)	
  
  
  ##################################################
  # PARAMETERS FOR SIMULATED LOCI
  ##################################################
  if(loc_parms$marker == "snp") {
    locus_params <- fscLocusParams(locus.type = as.character(loc_parms$marker), num.loci = loc_parms$nloci, mut.rate = loc_parms$seq_length*loc_parms$mu)
    attr(locus_params, "opts") <- "-s 0"
  } else if(loc_parms$marker == "dna") {
    locus_params <- fscLocusParams(locus.type = as.character(loc_parms$marker), sequence.length = loc_parms$seq_length, num.chrom = loc_parms$nloci, mut.rate = loc_parms$mu)
  } else if(loc_parms$marker == "msat") {
    locus_params <- fscLocusParams(locus.type = as.character(loc_parms$marker), num.loci = 1, num.chrom = loc_parms$nloci, mut.rate = loc_parms$mu)
  }
  
  
  
  ##################################################	
  # HISTORY OF DEMOGRAPHIC EVENTS
  ##################################################
  usemigmat <- 0
  ev <- 1
  
  for(G in 0:ngens) {
    tmp_pops <- coalhist[which(coalhist$time == G),]
    #tmp_pops <- FSCpops[which(FSCpops$arrive == G),]
    tmp_leaving <- c()
    
    if(length(tmp_pops[,1]) > 0) {
      for(fp in 1:length(tmp_pops[,1])) {
        if(!is.na(tmp_pops$snk[fp])) {
          ##!## Adding this if statement to check whether this is the ORIGINAL founding event - JDR 4/1/2020
          if(tmp_pops$prop[fp] == 1) {
            allevents[ev,] <- c(G-1, tmp_pops$src[fp], tmp_pops$src[fp], 0, as.numeric(found_Ne)/as.numeric(pop_size[tmp_pops$src[fp]+1]), 0, usemigmat)
            ev <- ev+1
            tmp_leaving <- c(tmp_leaving, tmp_pops$src[fp]+1)
          }
          ##!## Adding in tmp_pops$prop for migrant proportion - JDR 4/1/2020
          allevents[ev,] <- c(G, tmp_pops$src[fp], tmp_pops$snk[fp], tmp_pops$prop[fp], 1, 0, usemigmat+1)
          ev <- ev+1
        }
      }
      
      migmat[[usemigmat+2]] <- migmat[[usemigmat+1]]
      migmat[[usemigmat+2]][tmp_leaving,] <- 0
      migmat[[usemigmat+2]][,tmp_leaving] <- 0
      usemigmat <- usemigmat + 1
      
      rm(tmp_pops)
      rm(tmp_leaving)
      
    }
  }
  
  ##################################################
  # COMBINE EVERYTHING, SORT, AND RUN FSC 
  ##################################################
  allevents <- allevents[order(allevents$time, allevents$new_mig, allevents$prop_m, allevents$source),]	
  
  hist_ev <- fscHistEv(num.gen = allevents$time, source.deme = allevents$source, sink.deme = allevents$sink, prop.migrants = allevents$prop_m, new.sink.size = allevents$new_size, new.sink.growth = allevents$new_growth, new.mig.mat = allevents$new_mig)
  hist_ev1 <- fscHistEv(num.gen = ngens, source.deme = newID[which(simhist[,1] > 0)], sink.deme = newID[which(simhist[,1] > 0)], prop.migrants = 0, new.sink.size = 1, new.sink.growth = 0, new.mig.mat = "nomig")
  
  if(length(newID[which(simhist[,1] > 0)]) == 1) {
    hist_ev2 <- fscHistEv(num.gen = preLGMparms$preLGM_t[1], newID[which(simhist[,1] > 0)], sink.deme = newID[which(simhist[,1] > 0)], prop.migrants = 0, new.sink.size = round(preLGMparms$preLGM_Ne/preLGMparms$ref_Ne,5), new.sink.growth = 0, new.mig.mat = "nomig")
    hist_ev <- rbind(hist_ev, hist_ev1, hist_ev2)
  } else {
    arb.preLGM.pop <- newID[which(simhist[,1] > 0)][1]
    hist_ev2 <- fscHistEv(num.gen = preLGMparms$preLGM_t[1]-1, source.deme = newID[which(simhist[,1] > 0)][-1], sink.deme = arb.preLGM.pop, prop.migrants = 1, new.sink.size = 1, new.sink.growth = 0, new.mig.mat = "nomig")
    hist_ev3 <- fscHistEv(num.gen = preLGMparms$preLGM_t[1], source.deme = arb.preLGM.pop, sink.deme = arb.preLGM.pop, prop.migrants = 0, new.sink.size = round(preLGMparms$preLGM_Ne/preLGMparms$ref_Ne,5), new.sink.growth = 0, new.mig.mat = "nomig")
    hist_ev <- rbind(hist_ev, hist_ev1, hist_ev2, hist_ev3)
  }
  
  fscout <- fastsimcoal(label = label, pop.info = pop_info, locus.params = locus_params, mig.rates = migmat, hist.ev = hist_ev, num.cores = num_cores, delete.files = delete_files, exec = exec)
  #Naming strata prior to output
  FSCid <- sapply(strsplit(fscout@data$ids,"_"), function(x){x[1]})
  FSCgridid <- plyr::mapvalues(FSCid,newID,oldID,warn_missing=FALSE)
  ##!## Adding this in to facilitate mapping back to population abbreviations
  focalpops <- gmap[gmap$pop %in% l$sampdf$cell,]
  origgridid <- plyr::mapvalues(FSCgridid, focalpops$gpop, focalpops$pop, warn_missing = FALSE)
  #Then map these old ID's back to the abbreviation of the sampled population
  FSCabbrev <- plyr::mapvalues(origgridid,l$sampdf$cell,l$sampdf$abbrev)
  fscout@data$strata <- FSCabbrev
  fscout@data$ids <- paste0(fscout@data$strata, "_", c(1:length(fscout@data[,1])))

  fscout
  
}