simulate_data.R

library(cellAlign)
require(ggplot2)
library(geiger)
library(sp)
library(rgeos)

generate_data = function(slope, sd, population_size=10000){
"
Generates a population of data that has specified slope and deviation. The Y intercept is set to 1. 
It also generates a population of data to be used for showing no increase expression (a flat line). 
This flat linecomes from the average expression value of the original population
Parameters:
  slope: what the slope of the population should be. Mimics increaseing (or decreasing) protein expression
  sd: the standard devation of the data. Measures the biological variability in data
  population_size: how many data points to generate. Defaults to 1,000
  
Return: 
  A list, where time_traj is the pseudotime values, 
  Y is the expression values for the population with slope and sd, 
  and flat_Y is the expression values for the flat line
"
  X = runif(population_size)
  X = as.data.frame(sort(X))
  X = t(X)
  
  Ynorm <-  slope * X +1+ rnorm(population_size,sd = sd)
  
  Yflat = rep(mean(Ynorm), population_size)
  Yflat = t(Yflat)
  colnames(Ynorm) = paste0("t", 1:population_size)
  rownames(Ynorm) = "gene"
  colnames(Yflat) = paste0("t", 1:population_size)
  rownames(Yflat) = "gene"
  
  X = t(X)
  rownames(X) = paste0("t", 1:population_size)
  colnames(X) = NULL
  
  data = list("time_traj" = X, "Y" = Ynorm, "flat_Y"=Yflat)
  
  return(data)

}
subsample = function(population_data, number_of_cells){
"
Subsamples from the population generated by the generate_data() function. 
It splits that data into quarters and samples rnadomly form each one, to get an even distribution across time. 
Parameters:
  population_data is a dataframe that is the result of generate_data() that we're going to subsample from.
  number_of_cells the number of cells to subsample
Returns:
  List where traj is the pseudotime trajectory for the subsample, and exp is the protein expression values for the subsample
"
  traj = population_data$time_traj
  exp = population_data$Y
  trajdf = as.data.frame(traj)
  
  #divide time into quarters
  q = split(as.vector(trajdf), factor(sort(rank(as.vector(trajdf))%%4)))
  first = q$`0`
  second = q$`1`
  third = q$`2`
  fourth = q$`3`
  
  #get an even number of sample from each
  remainder = number_of_cells %% 4
  s_size = number_of_cells %/% 4
  s_size_plus1 = s_size 
  s_size_plus2 = s_size 
  s_size_plus3 = s_size 
  
  if (remainder == 1){
    s_size_plus1 = s_size_plus1 + 1 
  }
  if (remainder == 2){
    s_size_plus1 = s_size_plus1 + 1 
    s_size_plus2 =  s_size_plus2+1 
  }
  if (remainder == 3){
    s_size_plus1 = s_size_plus1 + 1 
    s_size_plus2 =  s_size_plus2+1 
    s_size_plus3 = s_size_plus3 +1 
  }
  
  first_index = rownames(first)
  first_sample = sample(first_index, s_size_plus1)
  
  second_index = rownames(second)
  second_sample = sample(second_index, s_size_plus2)
  
  third_index = rownames(third)
  third_sample = sample(third_index, s_size_plus3)
  
  fourth_index = rownames(fourth)
  fourth_sample = sample(fourth_index, s_size)
  
  all_samples = c(first_sample, second_sample, third_sample, fourth_sample)
  
  #subset traj
  trajdf <- cbind(index = rownames(trajdf), trajdf)
  rownames(trajdf) <- 1:nrow(trajdf)
  traj_subset = trajdf[which((trajdf$index %in% all_samples)==TRUE),]
  rownames(traj_subset) <- traj_subset$index
  traj_subset$index=NULL
  traj_subset = as.matrix(traj_subset)
  colnames(traj_subset) = NULL
  
  
  #subset exp
  exp_subset = exp[ ,which((colnames(exp) %in% all_samples)==TRUE)]
  exp_subset = t(exp_subset)
  rownames(exp_subset)="gene"
  exp_subset = as.data.frame(exp_subset)
  
  subsample = list("traj"=traj_subset, "exp"=exp_subset)
  return(subsample)
}
interpolate = function(time_traj, expression, subsample=TRUE, number_of_points=200){
"
Interpolates the data along the trajectory.
Data is interpolated by equally spaced points along the pseudotime trajectory. 

Paraeters:
  time_traj: the pseudo time trajectory
  expression: the protein expression values
  subsample: boolean that says weather we're interpolating a line for the subsample or the entire population. Defaults to TRUE
  number_of_points: the number of points to interpolate (the points will all be equally spaced). Defualts to 200
Return:
  A list where interpolated is a dataframe containg traj and exp for the interpolated line, 
  and points is a dataframe containing traj and exp data for the points (not an interpolated line)
"
  Y = expression
  X = time_traj
  interGlobalLPS = cellAlign::interWeights(expDataBatch = Y, trajCond =X,
                                           winSz = 0.1, numPts = number_of_points)
  require(reshape2)
  whichgene = "gene"
  selectedLPS<-interGlobalLPS$interpolatedVals[whichgene,]
  
  dfLPSi = data.frame(traj = interGlobalLPS$traj, value=(selectedLPS), error=interGlobalLPS$error[whichgene,])
  
  if(subsample==TRUE){
    dfLPS = data.frame(traj = X, gene=t(Y[whichgene,]))
  }else{
    dfLPS = data.frame(traj = X, gene=Y[whichgene,])
  }
  dfLPSM = melt(dfLPS, id.vars = 'traj')
  
  result = list("interpolated" = dfLPSi, "points"=dfLPSM)
  return(result)

}
find_area = function(full_line_curve, sample_curve){
"
Finds the area between two curves
Parameters: 
  full_line_curve is the curve that you want to compare the sampled line to (weather it be the oracel or the flat line)
  sample_curve is the curve of the sample. 
  If you wanted to compare a full flate line to a full oracle line, you could use this curve for that as well.
  The important thing is that sample_curve is used for what has the shortest line.
Returns:
  total_area: an number that is the total area between two curves
"
  y1 = full_line_curve$interpolated$value
  y2 = sample_curve$interpolated$value
  x1 = full_line_curve$interpolated$traj
  x2 = sample_curve$interpolated$traj
  # plot(x1, y1)
  # plot(x2, y2)
  
  # convert to a sp object (spatial lines)
  l1 <- Line(matrix(c(x1, y1), nc = 2, byrow = F))
  l2 <- Line(matrix(c(x2, y2), nc = 2, byrow = F))
  ll1 <- Lines(list(l1), ID = "1")
  ll2 <- Lines(list(l2), ID = "1")
  sl1 <- SpatialLines(list(ll1), proj4string = CRS("+init=epsg:4269"))
  sl2 <- SpatialLines(list(ll2), proj4string = CRS("+init=epsg:4269"))
  # Calculate locations where spatial lines intersect
  # if there is no area to be caculated (ex when we are calculating slope of stright line with sd of 0)
  if(all.equal(y1,y2)==TRUE){
    total_area = 0
    return(total_area)
  }else if(is.null(gIntersection(sl1, sl2, byid = TRUE))){ #if there are no interecting points
    smallest = min(sample_curve$interpolated$traj)
    largest = max(sample_curve$interpolated$traj)
    x = c()
    x = append(x, largest, after = length(x))
    x = append(x, smallest, after = length(x))
    x = sort(x, decreasing = FALSE)
    
  }else{ #if there are intersecting points
    int.pts <- gIntersection(sl1, sl2, byid = TRUE)
    int.coords <- int.pts@coords
    
    # Add end points
    smallest = min(sample_curve$interpolated$traj)
    largest = max(sample_curve$interpolated$traj)
    x = int.coords[,1]
    x = append(x, largest, after = length(x))
    x = append(x, smallest, after = length(x))
    x = sort(x, decreasing = FALSE)
  }
  
  total_area = 0
  for (i in c(1:length(x))){ #for everypoint that is an intersection
    first_x_coor = x[i]
    second_x_coor = x[i+1]
    if(is.na(second_x_coor)){
      break
    }else{
      area = geiger:::.area.between.curves(x1, y1, y2, xrange = c(first_x_coor, second_x_coor))
      area = abs(area)
      total_area = total_area + area
    }
  }
  total_area

  return(total_area)
}
accuracy_perc = function(oracle_area_list, flat_area_list){
  difference = flat_area_list - oracle_area_list
  neg = sum(difference < 0) #flat wins
  pos = sum(difference > 0) #truth wins
  accuracy = pos/length(difference)
  return(accuracy)
}
simulate = function(slope, sd, number_of_cells, number_of_subsamples=1000, population_size=10000){
'
Generates data population with given slope, sd. 
It draws a given number of subsamples from that population and (for each) finds the area between the sample and oracle, and the sample and flat line.
The areas are then used to find the percent accuarcy. 
Return:
  percent accuracy
'
  data = generate_data(slope=slope, sd=sd, population_size=population_size) #generate data
  oracle_line = interpolate(time_traj = data$time_traj, expression = data$Y, subsample = FALSE, number_of_points = 1000) #interpolate oracle
  flat_line = interpolate(time_traj = data$time_traj, expression = data$flat_Y, subsample = FALSE, number_of_points = 1000) #interpolate flat line
  
  sample_vs_oracle_list = c()
  sample_vs_flat_list = c()
  for(i in 1:number_of_subsamples){
    population_subsample = subsample(population_data = data, number_of_cells = number_of_cells)
    sample_line = interpolate(time_traj=population_subsample$traj, expression=population_subsample$exp, subsample=TRUE, number_of_points = 1000)
   
    sample_vs_oracle = find_area(full_line_curve = oracle_line, sample_curve = sample_line)
    sample_vs_oracle_list[i] = sample_vs_oracle
      
    sample_vs_flat = find_area(full_line_curve = flat_line, sample_curve = sample_line)
    sample_vs_flat_list[i] = sample_vs_flat
  }
  
  perc_accuracy = accuracy_perc(oracle_area_list=sample_vs_oracle_list, flat_area_list=sample_vs_flat_list)
  return(perc_accuracy)
}
simulate_false_positives = function(slope, sd, number_of_cells, number_of_subsamples=1000, population_size=10000){
  '
Return the percent of sample incorrectly classified as having a trend (false positive)

'
  ##Generate data (flat line is now oracle)
  X = runif(population_size)
  X = as.data.frame(sort(X))
  X = t(X)
  
  #generate the sloped line
  Yslope =slope * X 
  #generate flat line that will intersect at the midpoint
  y1 = min(Yslope)
  y2 = max(Yslope)
  midpoint = (y1+y2)/2
  Yflat <-  0 *X + midpoint + rnorm(population_size,sd = sd) #give the flat line variance (since it's representing the true population)
  # plot(X, Yflat)
  # plot(X, Yslope)
  
  #now we need to format it correctly
  colnames(Yflat) = paste0("t", 1:population_size)
  rownames(Yflat) = "gene"
  colnames(Yslope) = paste0("t", 1:population_size)
  rownames(Yslope) = "gene"
  X = t(X)
  rownames(X) = paste0("t", 1:population_size)
  colnames(X) = NULL
  
  flat_data = list("time_traj" = X, "Y" = Yflat) #oracle in this case is flat line
  slope_data = list("time_traj" = X, "Y" = Yslope)
  
  # time_traj, expression, subsample=TRUE, number_of_points=200
  #interpolate
  flat_line = interpolate(time_traj=X, expression=Yflat, subsample=FALSE, number_of_points = 1000) 
  slope_line = interpolate(time_traj=X, expression=Yslope, subsample=FALSE, number_of_points = 1000) 
  
  
  sample_vs_flat_list = c()
  sample_vs_slope_list = c()
  for(i in 1:number_of_subsamples){
    population_subsample = subsample(population_data = flat_data, number_of_cells = number_of_cells)
    sample_line = interpolate(time_traj=population_subsample$traj, expression=population_subsample$exp, subsample=TRUE, number_of_points = 1000)
    
    sample_vs_flat = find_area(full_line_curve = flat_line, sample_curve = sample_line)
    sample_vs_flat_list[i] = sample_vs_flat
    
    sample_vs_slope = find_area(full_line_curve = slope_line, sample_curve = sample_line)
    sample_vs_slope_list[i] = sample_vs_slope
  }
  
  perc_classified_as_pos = accuracy_perc(oracle_area_list=sample_vs_slope_list, flat_area_list=sample_vs_flat_list) #how many are (incorrectly) classified as having a trend
  return(perc_classified_as_pos)
}

#Functions for figure 1
simulate_fig_1_sd = function(slope, number_of_cells, population_size){
  sd0 = simulate(slope=slope, sd=0, number_of_cells = number_of_cells, population_size=10000)
  sd0p25 = simulate(slope=slope, sd=.25, number_of_cells = number_of_cells, population_size=10000)
  sd0p5 = simulate(slope=slope, sd=.5, number_of_cells = number_of_cells, population_size=10000)
  sd0p7 = simulate(slope=slope, sd=.7, number_of_cells = number_of_cells, population_size=10000)
  sd1 = simulate(slope=slope, sd=1, number_of_cells = number_of_cells, population_size=10000)
  df = data.frame("accuracy"=c(sd0, sd0p25, sd0p5, sd0p7, sd1), "deviation"=c(0, .25, .5,.7,1), "slope"=c(slope), "cells"=c(number_of_cells))
  return(df)
}
simulate_fig_1_cells = function(slope){
  cells_7 = simulate_fig_1_sd(slope= slope, number_of_cells = 7, population_size=10000)
  cells_16 = simulate_fig_1_sd(slope= slope, number_of_cells = 16, population_size=10000)
  cells_20 = simulate_fig_1_sd(slope= slope, number_of_cells = 20, population_size=10000)
  cells_30 = simulate_fig_1_sd(slope= slope, number_of_cells = 30, population_size=10000)
  cells_100 = simulate_fig_1_sd(slope= slope, number_of_cells = 100, population_size=10000)
  all_df = do.call("rbind", list(cells_7, cells_16, cells_20, cells_30, cells_100))
  return(all_df)
  
}
#Functions for figure 2
simulate_fig2_FP_sd = function(slope, number_of_cells, population_size){
  sd0 = simulate_false_positives(slope=slope, sd=0, number_of_cells = number_of_cells, population_size=10000)
  sd0p25 = simulate_false_positives(slope=slope, sd=.25, number_of_cells = number_of_cells, population_size=10000)
  sd0p5 = simulate_false_positives(slope=slope, sd=.5, number_of_cells = number_of_cells, population_size=10000)
  sd0p7 = simulate_false_positives(slope=slope, sd=.7, number_of_cells = number_of_cells, population_size=10000)
  sd1 = simulate_false_positives(slope=slope, sd=1, number_of_cells = number_of_cells, population_size=10000)
  df = data.frame("accuracy"=c(sd0, sd0p25, sd0p5, sd0p7, sd1), "deviation"=c(0, .25, .5,.7,1), "slope"=c(slope), "cells"=c(number_of_cells))
  return(df)
}
simulate_fig2_FP_cells = function(slope){
  cells_7 = simulate_fig2_FP_sd(slope= slope, number_of_cells = 7, population_size=10000)
  cells_16 = simulate_fig2_FP_sd(slope= slope, number_of_cells = 16, population_size=10000)
  cells_20 = simulate_fig2_FP_sd(slope= slope, number_of_cells = 20, population_size=10000)
  cells_30 = simulate_fig2_FP_sd(slope= slope, number_of_cells = 30, population_size=10000)
  cells_100 = simulate_fig2_FP_sd(slope= slope, number_of_cells = 100, population_size=10000)
  all_df = do.call("rbind", list(cells_7, cells_16, cells_20, cells_30, cells_100))
  return(all_df)
  
}
#Functions for figure 3A
simulate_fig3_cells = function(slope, sd){
  cells_7 = simulate(sd=sd, slope=slope, number_of_cells=7, number_of_subsamples=1000, population_size=10000) #changed
  cells_16 = simulate(sd=sd, slope=slope, number_of_cells=16, number_of_subsamples=1000, population_size=10000)#changed
  cells_20 = simulate(sd=sd, slope=slope, number_of_cells=20, number_of_subsamples=1000, population_size=10000)#changed
  cells_30 = simulate(sd=sd, slope=slope, number_of_cells=30, number_of_subsamples=1000, population_size=10000)#changed
  cells_100 = simulate(sd=sd, slope=slope, number_of_cells=100, number_of_subsamples=1000, population_size = 10000)
  all_cells = data.frame("accuracy"=c(cells_7, cells_16, cells_20, cells_30, cells_100), "deviation"=c(sd), "slope"=c(slope), "cells"=c(7,16,20,30,100))
  return(all_cells)
}
#Functions for figure 3B
simulate_fig3_FP_cells = function(slope, sd){
  cells_7 = simulate_false_positives(sd=sd, slope=slope, number_of_cells=7, number_of_subsamples=1000, population_size = 10000)
  cells_16 = simulate_false_positives(sd=sd, slope=slope, number_of_cells=16, number_of_subsamples=1000, population_size = 10000)
  cells_20 = simulate_false_positives(sd=sd, slope=slope, number_of_cells=20, number_of_subsamples=1000, population_size = 10000)
  cells_30 = simulate_false_positives(sd=sd, slope=slope, number_of_cells=30, number_of_subsamples=1000, population_size = 10000)
  cells_100 = simulate_false_positives(sd=sd, slope=slope, number_of_cells=100, number_of_subsamples=1000, population_size = 10000)
  all_cells = data.frame("accuracy"=c(cells_7, cells_16, cells_20, cells_30, cells_100), "deviation"=c(sd), "slope"=c(slope), "cells"=c(7,16,20,30,100))
  return(all_cells)
}


####Generate Figure 1####
#For 10 reps, you'd do this 10 times
dir.create("data")
dir.create("data/Fig1")
dir.create("data/Fig1/rep1")

slope0p5 = simulate_fig_1_cells(slope=.5)
write.csv(slope0p5, "data/Fig1/rep1/slope0p5")
slope1 = simulate_fig_1_cells(slope=1)
write.csv(slope1, "data/Fig1/rep1/slope1")
slope2 = simulate_fig_1_cells(slope=2)
write.csv(slope2, "data/Fig1/rep1/slope2")
slope4 = simulate_fig_1_cells(slope=4)
write.csv(slope4, "data/Fig1/rep1/slope4")

dir.create("data/Fig1/rep2")
slope0p5 = simulate_fig_1_cells(slope=.5)
write.csv(slope0p5, "data/Fig1/rep2/slope0p5")
slope1 = simulate_fig_1_cells(slope=1)
write.csv(slope1, "data/Fig1/rep2/slope1")
slope2 = simulate_fig_1_cells(slope=2)
write.csv(slope2, "data/Fig1/rep2/slope2")
slope4 = simulate_fig_1_cells(slope=4)
write.csv(slope4, "data/Fig1/rep2/slope4")

dir.create("data/Fig1/rep3")
slope0p5 = simulate_fig_1_cells(slope=.5)
write.csv(slope0p5, "data/Fig1/rep3/slope0p5")
slope1 = simulate_fig_1_cells(slope=1)
write.csv(slope1, "data/Fig1/rep3/slope1")
slope2 = simulate_fig_1_cells(slope=2)
write.csv(slope2, "data/Fig1/rep3/slope2")
slope4 = simulate_fig_1_cells(slope=4)
write.csv(slope4, "data/Fig1/rep3/slope4")

dir.create("data/Fig1/rep4")
slope0p5 = simulate_fig_1_cells(slope=.5)
write.csv(slope0p5, "data/Fig1/rep4/slope0p5")
slope1 = simulate_fig_1_cells(slope=1)
write.csv(slope1, "data/Fig1/rep4/slope1")
slope2 = simulate_fig_1_cells(slope=2)
write.csv(slope2, "data/Fig1/rep4/slope2")
slope4 = simulate_fig_1_cells(slope=4)
write.csv(slope4, "data/Fig1/rep4/slope4")

dir.create("data/Fig1/rep5")
slope0p5 = simulate_fig_1_cells(slope=.5)
write.csv(slope0p5, "data/Fig1/rep5/slope0p5")
slope1 = simulate_fig_1_cells(slope=1)
write.csv(slope1, "data/Fig1/rep5/slope1")
slope2 = simulate_fig_1_cells(slope=2)
write.csv(slope2, "data/Fig1/rep5/slope2")
slope4 = simulate_fig_1_cells(slope=4)
write.csv(slope4, "data/Fig1/rep5/slope4")

dir.create("data/Fig1/rep6")
slope0p5 = simulate_fig_1_cells(slope=.5)
write.csv(slope0p5, "data/Fig1/rep6/slope0p5")
slope1 = simulate_fig_1_cells(slope=1)
write.csv(slope1, "data/Fig1/rep6/slope1")
slope2 = simulate_fig_1_cells(slope=2)
write.csv(slope2, "data/Fig1/rep6/slope2")
slope4 = simulate_fig_1_cells(slope=4)
write.csv(slope4, "data/Fig1/rep6/slope4")

dir.create("data/Fig1/rep7")
slope0p5 = simulate_fig_1_cells(slope=.5)
write.csv(slope0p5, "data/Fig1/rep7/slope0p5")
slope1 = simulate_fig_1_cells(slope=1)
write.csv(slope1, "data/Fig1/rep7/slope1")
slope2 = simulate_fig_1_cells(slope=2)
write.csv(slope2, "data/Fig1/rep7/slope2")
slope4 = simulate_fig_1_cells(slope=4)
write.csv(slope4, "data/Fig1/rep7/slope4")

dir.create("data/Fig1/rep8")
slope0p5 = simulate_fig_1_cells(slope=.5)
write.csv(slope0p5, "data/Fig1/rep8/slope0p5")
slope1 = simulate_fig_1_cells(slope=1)
write.csv(slope1, "data/Fig1/rep8/slope1")
slope2 = simulate_fig_1_cells(slope=2)
write.csv(slope2, "data/Fig1/rep8/slope2")
slope4 = simulate_fig_1_cells(slope=4)
write.csv(slope4, "data/Fig1/rep8/slope4")

dir.create("data/Fig1/rep9")
slope0p5 = simulate_fig_1_cells(slope=.5)
write.csv(slope0p5, "data/Fig1/rep9/slope0p5")
slope1 = simulate_fig_1_cells(slope=1)
write.csv(slope1, "data/Fig1/rep9/slope1")
slope2 = simulate_fig_1_cells(slope=2)
write.csv(slope2, "data/Fig1/rep9/slope2")
slope4 = simulate_fig_1_cells(slope=4)
write.csv(slope4, "data/Fig1/rep9/slope4")

dir.create("data/Fig1/rep10")
slope0p5 = simulate_fig_1_cells(slope=.5)
write.csv(slope0p5, "data/Fig1/rep10/slope0p5")
slope1 = simulate_fig_1_cells(slope=1)
write.csv(slope1, "data/Fig1/rep10/slope1")
slope2 = simulate_fig_1_cells(slope=2)
write.csv(slope2, "data/Fig1/rep10/slope2")
slope4 = simulate_fig_1_cells(slope=4)
write.csv(slope4, "data/Fig1/rep10/slope4")


####Generate Figure 2 (False Positives)####
#For ten reps youd do this ten times
dir.create("data/Fig2")

dir.create("data/Fig2/rep1")
slope0p5 = simulate_fig2_FP_cells(slope=.5)
write.csv(slope0p5, "data/Fig2/rep1/slope0p5")
slope1 = simulate_fig2_FP_cells(slope=1)
write.csv(slope1, "data/Fig2/rep1/slope1")
slope2 = simulate_fig2_FP_cells(slope=2)
write.csv(slope2, "data/Fig2/rep1/slope2")
slope4 = simulate_fig2_FP_cells(slope=4)
write.csv(slope4, "data/Fig2/rep1/slope4")

dir.create("data/Fig2/rep2")
slope0p5 = simulate_fig2_FP_cells(slope=.5)
write.csv(slope0p5, "data/Fig2/rep2/slope0p5")
slope1 = simulate_fig2_FP_cells(slope=1)
write.csv(slope1, "data/Fig2/rep2/slope1")
slope2 = simulate_fig2_FP_cells(slope=2)
write.csv(slope2, "data/Fig2/rep2/slope2")
slope4 = simulate_fig2_FP_cells(slope=4)
write.csv(slope4, "data/Fig2/rep2/slope4")

dir.create("data/Fig2/rep3")
slope0p5 = simulate_fig2_FP_cells(slope=.5)
write.csv(slope0p5, "data/generated_data2/fig1_FP/rep3/slope0p5")
slope1 = simulate_fig2_FP_cells(slope=1)
write.csv(slope1, "data/Fig2/rep3/slope1")
slope2 = simulate_fig2_FP_cells(slope=2)
write.csv(slope2, "data/Fig2/rep3/slope2")
slope4 = simulate_fig2_FP_cells(slope=4)
write.csv(slope4, "data/Fig2/rep3/slope4")

dir.create("data/Fig2/rep4")
slope0p5 = simulate_fig2_FP_cells(slope=.5)
write.csv(slope0p5, "data/Fig2/rep4/slope0p5")
slope1 = simulate_fig2_FP_cells(slope=1)
write.csv(slope1, "data/Fig2/rep4/slope1")
slope2 = simulate_fig2_FP_cells(slope=2)
write.csv(slope2, "data/Fig2/rep4/slope2")
slope4 = simulate_fig2_FP_cells(slope=4)
write.csv(slope4, "data/Fig2/rep4/slope4")

dir.create("data/Fig2/rep5")
slope0p5 = simulate_fig2_FP_cells(slope=.5)
write.csv(slope0p5, "data/Fig2/rep5/slope0p5")
slope1 = simulate_fig2_FP_cells(slope=1)
write.csv(slope1, "data/Fig2/rep5/slope1")
slope2 = simulate_fig2_FP_cells(slope=2)
write.csv(slope2, "data/Fig2/rep5/slope2")
slope4 = simulate_fig2_FP_cells(slope=4)
write.csv(slope4, "data/Fig2/rep5/slope4")

dir.create("data/Fig2/rep6")
slope0p5 = simulate_fig2_FP_cells(slope=.5)
write.csv(slope0p5, "data/Fig2/rep6/slope0p5")
slope1 = simulate_fig2_FP_cells(slope=1)
write.csv(slope1, "data/Fig2/rep6/slope1")
slope2 = simulate_fig2_FP_cells(slope=2)
write.csv(slope2, "data/Fig2/rep6/slope2")
slope4 = simulate_fig2_FP_cells(slope=4)
write.csv(slope4, "data/Fig2/rep6/slope4")

dir.create("data/Fig2/rep7")
slope0p5 = simulate_fig2_FP_cells(slope=.5)
write.csv(slope0p5, "data/Fig2/rep7/slope0p5")
slope1 = simulate_fig2_FP_cells(slope=1)
write.csv(slope1, "data/Fig2/rep7/slope1")
slope2 = simulate_fig2_FP_cells(slope=2)
write.csv(slope2, "data/Fig2/rep7/slope2")
slope4 = simulate_fig2_FP_cells(slope=4)
write.csv(slope4, "data/Fig2/rep7/slope4")

dir.create("data/Fig2/rep8")
slope0p5 = simulate_fig2_FP_cells(slope=.5)
write.csv(slope0p5, "data/Fig2/rep8/slope0p5")
slope1 = simulate_fig2_FP_cells(slope=1)
write.csv(slope1, "data/Fig2/rep8/slope1")
slope2 = simulate_fig2_FP_cells(slope=2)
write.csv(slope2, "data/Fig2/rep8/slope2")
slope4 = simulate_fig2_FP_cells(slope=4)
write.csv(slope4, "data/Fig2/rep8/slope4")

dir.create("data/Fig2/rep9")
slope0p5 = simulate_fig2_FP_cells(slope=.5)
write.csv(slope0p5, "data/Fig2/rep9/slope0p5")
slope1 = simulate_fig2_FP_cells(slope=1)
write.csv(slope1, "data/Fig2/rep9/slope1")
slope2 = simulate_fig2_FP_cells(slope=2)
write.csv(slope2, "data/Fig2/rep9/slope2")
slope4 = simulate_fig2_FP_cells(slope=4)
write.csv(slope4, "data/Fig2/rep9/slope4")

dir.create("data/Fig2/rep10")
slope0p5 = simulate_fig2_FP_cells(slope=.5)
write.csv(slope0p5, "data/Fig2/rep10/slope0p5")
slope1 = simulate_fig2_FP_cells(slope=1)
write.csv(slope1, "data/Fig2/rep10/slope1")
slope2 = simulate_fig2_FP_cells(slope=2)
write.csv(slope2, "data/Fig2/rep10/slope2")
slope4 = simulate_fig2_FP_cells(slope=4)
write.csv(slope4, "data/Fig2/rep10/slope4")

####Generate Figure 3A####
dir.create("data/Fig3A")
#slope_over_var=.5 ::  .5/1   .75/1.5  1/2  1.5/3   2/4   3/6
dir.create("data/Fig3A/slope_var_0p5")
  for(i in 1:10){
  combo1 = simulate_fig3_cells(slope=.5, sd=1)
  combo2 = simulate_fig3_cells(slope=.75, sd=1.5)
  combo3 = simulate_fig3_cells(slope=1, sd=2)
  combo4 = simulate_fig3_cells(slope=2, sd=4)
  combo5 = simulate_fig3_cells(slope=3, sd=6)
  slope_var_0p5 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3A/slope_var_0p5", "/rep", i, sep=""))
  write.csv(slope_var_0p5, "data/Fig3A/slope_var_0p5")
}

#slope_var = 1 ::   .5/.5  1/1  1.75/1.75  2/2   3/3
dir.create("data/Fig3A/slope_var_1")
  for(i in 1:10){
  combo1 = simulate_fig3_cells(slope=.5, sd=.5)
  combo2 = simulate_fig3_cells(slope=1, sd=1)
  combo3 = simulate_fig3_cells(slope=1.75, sd=1.75)
  combo4 = simulate_fig3_cells(slope=2, sd=2)
  combo5 = simulate_fig3_cells(slope=3, sd=3)
  slope_var_1 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3A/slope_var_1", "/rep", i, sep=""))
  write.csv(slope_var_1, file_path)
}

#Slope over Var 1.5 ::  # .75/.5  2.25/1.5  1.5/1   3/2   4.5/3  
dir.create("data/Fig3A/slope_var_1p5")
for(i in 1:10){
  combo1 = simulate_fig3_cells(slope=.75, sd=.5)
  combo2 = simulate_fig3_cells(slope=2.25, sd=1.5)
  combo3 = simulate_fig3_cells(slope=1.5, sd=1)
  combo4 = simulate_fig3_cells(slope=3, sd=2)
  combo5 = simulate_fig3_cells(slope=4.5, sd=3)
  slope_var_1p5 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3A/slope_var_1p5", "/rep", i, sep=""))=
  write.csv(slope_var_1p5, file_path)
}

#Slope over var = 2 ::     1/.5    1.5/.75   2/1   4/2   6/3
dir.create("data/Fig3A/slope_var_2")
for(i in 1:10){
  combo1 = simulate_fig3_cells(slope=1, sd=.5)
  combo2 = simulate_fig3_cells(slope=1.5, sd=.75)
  combo3 = simulate_fig3_cells(slope=2, sd=1)
  combo4 = simulate_fig3_cells(slope=4, sd=2)
  combo5 = simulate_fig3_cells(slope=6, sd=3)
  slope_var_2 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3A/slope_var_2", "/rep", i, sep=""))=
  write.csv(slope_var_2, file_path)
}
#Slope over var = 3 ::     1.5/.5   2.25/.75  3/1   6/2   9/3
dir.create("data/Fig3A/slope_var_3")
  for(i in 1:10){
  combo1 = simulate_fig3_cells(slope=1.5, sd=.5)
  combo2 = simulate_fig3_cells(slope=2.25, sd=.75)
  combo3 = simulate_fig3_cells(slope=3, sd=1)
  combo4 = simulate_fig3_cells(slope=6, sd=2)
  combo5 = simulate_fig3_cells(slope=9, sd=3)
  slope_var_3 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3A/slope_var_3", "/rep", i, sep=""))=
  write.csv(slope_var_3, file_path)
}
# Slope over var = 4  ::  # 2/.5    3/.75   4/1   8/2   12/3
dir.create("data/Fig3A/slope_var_4")
for(i in 1:10){
  combo1 = simulate_fig3_cells(slope=2, sd=.5)
  combo2 = simulate_fig3_cells(slope=3, sd=.75)
  combo3 = simulate_fig3_cells(slope=4, sd=1)
  combo4 = simulate_fig3_cells(slope=8, sd=2)
  combo5 = simulate_fig3_cells(slope=12, sd=3)
  slope_var_4 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3A/slope_var_4", "/rep", i, sep=""))=
  write.csv(slope_var_4, file_path)
}

#Slope over var = 6 ::    3/.5    4.5/.75   6/1   12/2  18/3
dir.create("data/Fig3A/slope_var_6")
  for(i in 1:10){
  combo1 = simulate_fig3_cells(slope=3, sd=.5)
  combo2 = simulate_fig3_cells(slope=4.5, sd=.75)
  combo3 = simulate_fig3_cells(slope=6, sd=1)
  combo4 = simulate_fig3_cells(slope=12, sd=2)
  combo5 = simulate_fig3_cells(slope=18, sd=3)
  slope_var_6 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3A/slope_var_6", "/rep", i, sep=""))=
  write.csv(slope_var_6, file_path)
}

#Slope over var = 10 ::   5/.5    7.5/.75   10/1  20/2   30/3
dir.create("data/Fig3A/slope_var_10")
  for(i in 1:10){
  combo1 = simulate_fig3_cells(slope=5, sd=.5)
  combo2 = simulate_fig3_cells(slope=7.5, sd=.75)
  combo3 = simulate_fig3_cells(slope=10, sd=1)
  combo4 = simulate_fig3_cells(slope=20, sd=2)
  combo5 = simulate_fig3_cells(slope=30, sd=3)
  slope_var_10 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3A/slope_var_10", "/rep", i, sep=""))
  write.csv(slope_var_10, file_path)
}

####Generate Figure 3B (False Positives)####
dir.create("data/Fig3B")
##slope/var = .5:  .25/.5    .5/1    1.5/3   2/4   3/6 ###Each one takes like 35 minutes
dir.create("data/Fig3B/slope_var_0p5")
for(i in 1:10){
  combo1 = simulate_fig3_FP_cells(slope=.5, sd=1)
  combo2 = simulate_fig3_FP_cells(slope=.75, sd=1.5)
  combo3 = simulate_fig3_FP_cells(slope=1, sd=2)
  combo4 = simulate_fig3_FP_cells(slope=2, sd=4)
  combo5 = simulate_fig3_FP_cells(slope=3, sd=6)
  slope_var_0p5 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3B/slope_var_0p5/", "/rep", i, sep=""))
  write.csv(slope_var_0p5, file_path)
}

##slope/var = 1:  .25/.25    .5/.5    1/1    2/2   3/3
dir.create("data/Fig3B/slope_var_1")
for(i in 1:10){
  combo1 = simulate_fig3_FP_cells(slope = .25, sd=.25)
  combo2 = simulate_fig3_FP_cells(slope = .5, sd=.5)
  combo3 = simulate_fig3_FP_cells(slope = 1, sd=1)
  combo4 = simulate_fig3_FP_cells(slope = 2, sd=2)
  combo5 = simulate_fig3_FP_cells(slope = 3, sd=3)
  slope_var_1 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3B/slope_var_1/", "/rep", i, sep=""))
  write.csv(slope_var_1, file_path)
}

##slope/var = 2:  .5/.25    1/.5    2/1    4/2   6/3
dir.create("data/Fig3B/slope_var_2")
for(i in 1:10){
  combo1 = simulate_fig3_FP_cells(slope = .5, sd=.25)
  combo2 = simulate_fig3_FP_cells(slope = 1, sd=.5)
  combo3 = simulate_fig3_FP_cells(slope = 2, sd=1)
  combo4 = simulate_fig3_FP_cells(slope = 4, sd=2)
  combo5 = simulate_fig3_FP_cells(slope = 6, sd=3)
  slope_var_2 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3B/slope_var_2/", "/rep", i, sep=""))
  write.csv(slope_var_2, file_path)
}

##slope/var = 4:  # 2/.5    3/.75   4/1   8/2   12/3
dir.create("data/Fig3B/slope_var_4")
  for(i in 1:10){
  combo1 = simulate_fig3_FP_cells(slope = 2, sd=.5)
  combo2 = simulate_fig3_FP_cells(slope = 3, sd=.75)
  combo3 = simulate_fig3_FP_cells(slope = 4, sd=1)
  combo4 = simulate_fig3_FP_cells(slope = 8, sd=2)
  combo5 = simulate_fig3_FP_cells(slope = 12, sd=3)
  slope_var_4 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3B/slope_var_4/", "/rep", i, sep=""))
  write.csv(slope_var_4, file_path)
}

##slope/var = 6:   3/.5    4.5/.75   6/1   12/2  18/3
dir.create("data/Fig3B/slope_var_6")
  for(i in 1:10){
  combo1 = simulate_fig3_FP_cells(slope = 3, sd=.5)
  combo2 = simulate_fig3_FP_cells(slope = 4.5, sd=.75)
  combo3 = simulate_fig3_FP_cells(slope = 6, sd=1)
  combo4 = simulate_fig3_FP_cells(slope = 12, sd=2)
  combo5 = simulate_fig3_FP_cells(slope = 18, sd=3)
  slope_var_6 = do.call("rbind", list(combo1,combo2,combo_3,combo_4,combo_5))
  file_path = as.character(paste("data/Fig3B/slope_var_6/", "/rep", i, sep=""))
  write.csv(slope_var_6, file_path)
  }

##slope/var = 10:   5/.5    7.5/.75   10/1  20/2   30/3
dir.create("data/Fig3B/slope_var_10")
  for(i in 1:10){
  combo1 = simulate_fig3_FP_cells(slope = 5, sd=.5)
  combo2 = simulate_fig3_FP_cells(slope = 7.5, sd=.75)
  combo3 = simulate_fig3_FP_cells(slope = 10, sd=1)
  combo4 = simulate_fig3_FP_cells(slope = 20, sd=2)
  combo5 = simulate_fig3_FP_cells(slope = 30, sd=3)
  slope_var_10 = do.call("rbind", list(combo1,combo2,combo3,combo4,combo5))
  file_path = as.character(paste("data/Fig3B/slope_var_10/", "/rep", i, sep=""))
  write.csv(slope_var_10, file_path)
}