<commit_before>require(ggplot2)
require(dplyr)
require(tidyr)

##############################################################################
# Function checks if elements in a exist in b
# Output is binary vector of length same as b
##############################################################################
toBin <- function(a,b){
	as.numeric(is.element(b,a))
}

##############################################################################
# Read data
##############################################################################
# GoNL DNMs
cat("Reading data...\n")
dnms <- read.table("/net/bipolar/jedidiah/mutation/reference_data/GoNL_DNMs.txt", header=T)
chrdnms <- dnms[dnms$CHROM==4,]

# Read predicted mutation rate data
chrp <- read.table("/net/bipolar/jedidiah/mutation/output/predicted/chr4_full_mask.txt", header=F)
names(chrp) <- c("POS", "MU")

# Index DNMs in site data
chrp$OBS <- toBin(chrdnms$POS, chrp$POS)

# Sort by logit-estimated rate and add column of cumulate proportion
chrp <- chrp %>% arrange(MU)
chrp$prop <- cumsum(chrp$OBS)/sum(chrp$OBS)

##############################################################################
# Simulate mutations by randomly selecting sites and checking if
# Bernoulli(mu)==1; loop continues until we have simulated the same number
# of sites as in observed data
##############################################################################
nsites <- sum(chrp$OBS)

for(i in 1:5){
	cat("Running simulation ", i, "of 10...\n")
	nsample <- 50000
	mutated <- c()

	while(length(mutated) < nsites){
		rowind <- sample(nrow(chrp), nsample)
		row <- chrp[rowind,]
		mu <- row$MU

		batch <- sapply(mu, function(x) rbinom(1,1,x))
		mutated <- c(mutated, row$POS[which(as.logical(batch))])
	}

	# chrp$TMP <- toBin(mutated[1:nsites], chrp$POS)
	chrp$TMP <- toBin(mutated, chrp$POS)
	colnames(chrp)[ncol(chrp)] <- paste0("SIMOBS", i)
}

##############################################################################
# Get cumulative proportions of simulated data
##############################################################################
for(i in 5:9){
	chrp$tmpprop <- cumsum(chrp[,i])/sum(chrp[,i])
	colnames(chrp)[ncol(chrp)] <- paste0("simprop", i)
}

##############################################################################
# Subsample 10M sites, resort by mu, add percentile column, and coerce
# from wide to long format
#
# Can delete chrp after this step
##############################################################################
cat("Subsampling data...\n")
chrp2<-chrp[sample(nrow(chrp), 1000000),] %>%
  arrange(MU) %>%
  mutate(ntile=ntile(MU, 1000)) %>%
  gather(group, val, .dots=grep("prop", names(chrp)))

cat("Cleaning up memory...\n")
rm(chrp)
gc()

##############################################################################
# Further summarize data--get means by ntile
#
# Can delete chrp2 after this step
##############################################################################
chrp3<-chrp2 %>%
  group_by(group, ntile) %>%
  summarise(val=mean(val))

##############################################################################
# Get subsets of true and simulated data
##############################################################################
chrp3a<-chrp3[chrp3$group!="prop",]
chrp3m<-chrp3[chrp3$group=="prop",]

##############################################################################
# Calculate lower/upper 95% CIs from simulations
##############################################################################
chrp3lb<-chrp3a %>%
  group_by(ntile) %>%
  summarise(lb=t.test(val)$conf.int[1]) %>%
  mutate(group="lb") %>%
  select(group, ntile, lb)

chrp3ub<-chrp3a %>%
  group_by(ntile) %>%
  summarise(ub=t.test(val)$conf.int[2]) %>%
  mutate(group="ub") %>%<commit_msg>finalize roc script<commit_after>require(ggplot2)
require(dplyr)
require(tidyr)

##############################################################################
# Function checks if elements in a exist in b
# Output is binary vector of length same as b
##############################################################################
toBin <- function(a,b){
	as.numeric(is.element(b,a))
}

##############################################################################
# Read data
##############################################################################
# GoNL DNMs
cat("Reading data...\n")
dnms <- read.table("/net/bipolar/jedidiah/mutation/reference_data/GoNL_DNMs.txt", header=T)
chrdnms <- dnms[dnms$CHROM==4,]

# Read predicted mutation rate data
chrp <- read.table("/net/bipolar/jedidiah/mutation/output/predicted/chr4_full_mask.txt", header=F)
names(chrp) <- c("POS", "MU")

# Index DNMs in site data
chrp$OBS <- toBin(chrdnms$POS, chrp$POS)

# Sort by logit-estimated rate and add column of cumulate proportion
chrp <- chrp %>% arrange(MU)
chrp$prop <- cumsum(chrp$OBS)/sum(chrp$OBS)

##############################################################################
# Simulate mutations by randomly selecting sites and checking if
# Bernoulli(mu)==1; loop continues until we have simulated the same number
# of sites as in observed data
##############################################################################
nsites <- sum(chrp$OBS)

for(i in 1:5){
	cat("Running simulation ", i, "of 10...\n")
	nsample <- 20000
	mutated <- c()

	while(length(mutated) < 400){
		rowind <- sample(nrow(chrp), nsample)
		row <- chrp[rowind,]
		mu <- row$MU

		batch <- sapply(mu, function(x) rbinom(1,1,x))
		mutated <- c(mutated, row$POS[which(as.logical(batch))])
	}

	# chrp$TMP <- toBin(mutated[1:nsites], chrp$POS)
	chrp$TMP <- toBin(mutated, chrp$POS)
	colnames(chrp)[ncol(chrp)] <- paste0("SIMOBS", i)
}

##############################################################################
# Get cumulative proportions of simulated data
##############################################################################
for(i in 5:9){
	chrp$tmpprop <- cumsum(chrp[,i])/sum(chrp[,i])
	colnames(chrp)[ncol(chrp)] <- paste0("simprop", i)
}

##############################################################################
# Subsample 10M sites, resort by mu, add percentile column, and coerce
# from wide to long format
#
# Can delete chrp after this step
##############################################################################
cat("Subsampling data...\n")
chrp2<-chrp[sample(nrow(chrp), 1000000),] %>%
  arrange(MU) %>%
  mutate(ntile=ntile(MU, 1000)) %>%
  gather(group, val, .dots=grep("prop", names(chrp)))

cat("Cleaning up memory...\n")
rm(chrp)
gc()

##############################################################################
# Further summarize data--get means by ntile
#
# Can delete chrp2 after this step
##############################################################################
chrp3<-chrp2 %>%
  group_by(group, ntile) %>%
  summarise(val=mean(val))

##############################################################################
# Get subsets of true and simulated data
##############################################################################
chrp3a<-chrp3[chrp3$group!="prop",]
chrp3m<-chrp3[chrp3$group=="prop",]

##############################################################################
# Calculate lower/upper 95% CIs from simulations
##############################################################################
chrp3lb<-chrp3a %>%
  group_by(ntile) %>%
  # summarise(lb=t.test(val)$conf.int[1]) %>%
	summarise(lb=max(val)) %>%
  mutate(group="lb") %>%
  select(group, ntile, lb)

chrp3ub<-chrp3a %>%
  group_by(ntile) %>%
  # summarise(ub=t.test(val)$conf.int[2]) %>%
	summarise(ub=min(val)) %>%
  mutate(group="ub") %>%
{'fraction_non_alphanumeric': 0.4437181330892927, 'fraction_numerical': 0.01372728461236763, 'mean_word_length': 7.585858585858586, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '686756', 'n_tokens_mistral': 2459, 'n_tokens_neox': 2171, 'n_words': 458}