pairwise_comparisons_conditional_selection_simulated_cohorts

Here, we will perform pairwise comparisons of conditional seleciton on simulated cohorts

For a detailed description of the resampling/bootstrapping algorithm for pairwise comparisons, please refer to pseudo code

For a detailed description of the algorithm for simulating genes, please refer to methods

Plugging Data into simresults generator #2

These simulation have been optimized to used parallel computing. Please note that true_OR refers to the odds ratio of the gene of interest with a positive control
Please note that or_pair2 refers to teh odds ratio between the two positive controls (generally low)
Note to self: just convert all the “true OR” variables to or_pair1

##GenePair1
cores=detectCores()
cl= makeCluster(cores[1]-1)
registerDoParallel(cl)
tic()
or_pair1=c(.01,.05,.1,.2,.3,.4,.5,.6,.7,.8,.9,1)
or_pair2=seq(.01,.1,by=.02)
incidence=seq(4,36,by=1)
cohort_size=seq(100,1000,by=200)
# unlist(mut_excl_genes_generator(500,25,.4,.01)[1])

# or_pair1=c(.01,.05,.1,.2,.3,.4,.5,.6,.7,.8,.9,1)
# or_pair2=seq(.01,.01)
# incidence=seq(4,50,by=1)
# cohort_size=seq(500,500)

# or_pair1=.8
# or_pair2=.01
# incidence=25
# cohort_size=500

true_or_vals=or_pair1
cohort_size_vals=cohort_size
gene1_total_vals=incidence

# simresults_compiled=as.data.frame(matrix(nrow=1,ncol=12)) #10 because there is an OR vals column being added at the end of the simulations
simresults_compiled=as.list(seq_len(length(or_pair2)*length(or_pair1)*length(incidence)*length(cohort_size)))

ct=1
# for(j in 1:length(gene1_total_vals)){
parfunc=foreach(j =1:length(gene1_total_vals),.combine = rbind)%dopar%{
for(l in 1:length(or_pair2)){
for(k in 1:length(cohort_size_vals)){
    for(i in 1:length(true_or_vals)){
gene_pair_1=unlist(mut_excl_genes_generator(cohort_size[k],incidence[j],or_pair1[i],or_pair2[l])[1])
gene_pair_2=unlist(mut_excl_genes_generator(cohort_size[k],incidence[j],or_pair1[i],or_pair2[l])[2])
alldata_1=mut_excl_genes_datapoints(gene_pair_1)
alldata_2=mut_excl_genes_datapoints(gene_pair_2)
alldata_comp_1=alldata_compiler(alldata_1,"gene2","gene3","gene1",'N',"N/A","N/A")[[2]]
genex_replication_prop_1=alldata_compiler(alldata_1,"gene2","gene3","gene1",'N',"N/A","N/A")[[1]]
alldata_comp_2=alldata_compiler(alldata_2,"gene2","gene3","gene1",'N',"N/A","N/A")[[2]]
genex_replication_prop_2=alldata_compiler(alldata_2,"gene2","gene3","gene1",'N',"N/A","N/A")[[1]]
#Number of simulations are the minimum of: subsampling with replacement combinations of incidence of gene pair 1, pair 2, and 10000. 10000 because I don't want to do more than 1k simulations per experiment.
nsims=min(c(choose(sum(alldata_comp_1$genex)+sum(alldata_comp_1$genex)-1,sum(alldata_comp_1$genex)),
  choose(sum(alldata_comp_2$genex)+sum(alldata_comp_2$genex)-1,sum(alldata_comp_2$genex)),1000))
# nsubsamples=sum(alldata_comp_2$genex) #Verify that this is what it is.    

# Max subsample size is whichever integer is the minimum of:
# number of genex in sample vs. number of NOTgenex or PC1 or PC2 in sample*genexreplication proportion. For example we had 165 PC1s, therefore 340-175=165 PC1s in our data. Since we want to sample 7.35% positive hits and 92.65% negative hits, a sample containing too many positive hits can impose an upper bound on our subsampling size.
alldata_comp=alldata_comp_1
genex_replication_prop=genex_replication_prop_1
maxsubsamplesize_1=min(c(genex_replication_prop*length(alldata_comp$Positive_Ctrl1),
                       floor((length(alldata_comp$Positive_Ctrl1)-sum(alldata_comp$Positive_Ctrl1))*genex_replication_prop),
                       floor((length(alldata_comp$Positive_Ctrl1)-sum(alldata_comp$Positive_Ctrl2))*genex_replication_prop),
                       floor((length(alldata_comp$Positive_Ctrl1)-sum(alldata_comp$genex))*genex_replication_prop)))
alldata_comp=alldata_comp_2
genex_replication_prop=genex_replication_prop_2
maxsubsamplesize_2=min(c(genex_replication_prop*length(alldata_comp$Positive_Ctrl1),
                       floor((length(alldata_comp$Positive_Ctrl1)-sum(alldata_comp$Positive_Ctrl1))*genex_replication_prop),
                       floor((length(alldata_comp$Positive_Ctrl1)-sum(alldata_comp$Positive_Ctrl2))*genex_replication_prop),
                       floor((length(alldata_comp$Positive_Ctrl1)-sum(alldata_comp$genex))*genex_replication_prop)))

nsubsamples=min(c(maxsubsamplesize_1,maxsubsamplesize_2)) #minimum of either of the subsamples.
nexperiments=7

#####If Cohort size is greater than gene1total, run simresults_generator for that simulation####

    if(cohort_size_vals[k]>gene1_total_vals[j]){

      simresults_pair1=simresults_generator(alldata_comp_1,7)
      simresults_pair1$gene_pair=1
      
      simresults_pair2=simresults_generator(alldata_comp_2,7)
      simresults_pair2$gene_pair=2
      
      simresults=rbind(simresults_pair1,simresults_pair2)
      simresults=as.data.frame(simresults, stringsAsFactors = F)

simresults$true_OR=true_or_vals[i]
simresults$gene1_total=gene1_total_vals[j]
simresults$cohort_size=cohort_size_vals[k]
simresults$or_pair2=or_pair2[l]

# colnames(simresults_compiled)=names(simresults)
# simresults_compiled[[i+j+k+l-3]]=simresults
simresults_compiled[[ct]]=simresults
# simresults_compiled=rbind(simresults_compiled,simresults)

    }
ct=ct+1
}
}
}
simresults_compiled
}
simresults_compiled=parfunc
simresults_compiled2=do.call(rbind,simresults_compiled)
stopCluster(cl)
toc()
# tic()
# pair1=simresults_generator(alldata_comp_1,7)
# pair1$gene_pair=1
# pair2=simresults_generator(alldata_comp_1,7)
# pair2$gene_pair=2
# aaa=rbind(pair1,pair2)
# a=c()
# a=rbind(a,aaa)
# toc()
# write.csv(simresults_compiled2,"alkati_simulations_compiled_12219.csv")
write.csv(simresults_compiled2,"data/alkati_simulations_compiled_1000_3719.csv")
# simresults_compiled[[5]]
#Looking at which conditions didn't yield a maxsimnumber of 1000 revealed that at a subsample size of 1, the max simnum was 35.
# a=simresults_compiled2%>%
#   group_by(subsample_size,exp_num,gene_pair,gene1_total,cohort_size,or_pair2)%>%
#   summarize(maxsimnum=max(sim_num))