In [1]:
library(tadaA)
setwd("/project/wgs/analysis/TADA/TADA-A")
Adjust mutation rate¶
In [2]:
TADA_A_adjust_mutrate(window_file = 'test_data/test_windows.txt',
mut_file = 'test_data/test_simulated_DNM.bed',
genes = "all",sample_size = 5000,
mutrate_mode = "regular",
scaling_file_name = "test_data/test_simulated_windows_mutrate_with_div_score_scaling_file_for_test_DNM.txt",
scale_features = c("GC_content", "div_score"))
Read in DNM data, annotation data and mutation rate data¶
In [3]:
system.time(compact_data <- TADA_A_read_info(
mut_files = c("test_data/test_simulated_DNM_with_allele_info.txt"),
window_file = "test_data/test_windows.txt",
mutrate_scaling_files = c("test_data/test_simulated_windows_mutrate_with_div_score_scaling_file_for_test_DNM.txt"),
sample_sizes = c(5000),
gene_prior_file = "test_data/test_gene_prior.txt",
nonAS_noncoding_annotations = "test_data/test_non_allele_specific_annotation_1.bed",
AS_noncoding_annotations = NA,
report_proportion = 2000/18665,
chunk_partition_num =1,
node_n = 1,
mutrate_ref_files = c("test_data/test_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_A.mutrate.bw",
"test_data/test_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_C.mutrate.bw",
"test_data/test_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_G.mutrate.bw" ,
"test_data/test_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_T.mutrate.bw")))
In [17]:
length(compact_data$base_info)
In [16]:
compact_data$base_info$ACTL8
Estimate relative risks¶
In [4]:
RR = TADA_A_RR_estimate(data = compact_data$base_info,
selected_annotations = c(1), gene_prior_file = "test_data/test_gene_prior.txt",
optimization_iteration = 2000)
In [5]:
RR
Bayesian FDR control¶
In [18]:
Bayesian_FDR(BF=c(seq(1,10)), pi0=0.05, alpha=0.05)
Bayes factor calculation¶
- This function is for calculating Bayes factors for based on noncoding annotations and then combine with coding Bayes factors in order to predict risk genes.
- calculate the Bayes factor for each gene based on informative noncodng annotations.
parameters:
- [data] is the [base_info] returned from [TADA_A_read_info], each element of the list has compact base-level information for each gene.
- [selected_annotations] determines which features in the [data] are going to be used for calculating bayes factors. e.g., if [selected_annotations] is c(1,3). Then the first and third element of feature_vector in the [base_info] of [data] will be selected.
- [rr] is the estimated relative risks of [selected_annotations]. RRs are estimated from [TADA_A_RR_estimate].
- RRs are in the log scale
- [additional_BF_file], a file that has additional Bayes factors that could be uesd along with noncoding-derived BFs to increase the power of risk gene mapping. e.g., "../data/Example_gene_coding_BF.txt"
- [TADA_p0] the proportion of genes being non-risk genes. Default = 0.94 for in the context of ASD risk gene mapping
Note:
%*% : non-conformable arguments
- There is a safer method for handling the "matrix-multiplication" of vectors, the crossprod function.
In [76]:
TADA_A_get_BFs(data = compact_data$base_info,
selected_annotations = c(1),
rr = RR$rr_report$relative_risk,
additional_BF_file = "tmp/tmp_gene_ad_bf",
TADA_p0 = 0.94)
TADA_A_DNM_generator¶
In [ ]:
window_file = 'test_data/test_windows.txt',
mut_file = 'test_data/test_simulated_DNM.bed',
genes = "all",sample_size = 5000,
mutrate_mode = "regular",
scaling_file_name = "test_data/test_simulated_windows_mutrate_with_div_score_scaling_file_for_test_DNM.txt",
scale_features = c("GC_content", "div_score")
- [data] is the [base_info] returned from [TADA_A_read_info], each element of the list has compact base-level information for each gene.
- [selected_annotations] determines which features in the [data] are going to be used for calculating bayes factors. e.g., if [selected_annotations] is c(1,3). Then the first and third element of feature_vector in the [base_info] of [data] will be selected.
- [rr] is the estimated relative risks of [selected_annotations]. RRs are estimated from [TADA_A_RR_estimate].
- RRs are in the log scale
- [additional_BF_file], a file that has additional Bayes factors that could be uesd along with noncoding-derived BFs to increase the power of risk gene mapping. e.g., "../data/Example_gene_coding_BF.txt"
- [TADA_p0] the proportion of genes being non-risk genes. Default = 0.94 for in the context of ASD risk gene mapping
In [103]:
TADA_A_DNM_generator(window_file = 'test_data/test_windows.txt',
mutrate_scaling_files = "test_data/test_simulated_windows_mutrate_with_div_score_scaling_file_for_test_DNM.txt",
sample_sizes = c(27),
gene_prior_file = "test_data/test_gene_prior.txt",
nonAS_noncoding_annotations = "test_data/test_non_allele_specific_annotation_1.bed",
AS_noncoding_annotations = NA,
report_proportion = 100/18665,
chunk_partition_num =1 ,
node_n = 6,
mutrate_ref_files = c("test_data/test_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_A.mutrate.bw",
"test_data/test_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_C.mutrate.bw",
"test_data/test_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_G.mutrate.bw",
"test_data/test_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_T.mutrate.bw"),
rr = RR$rr_report$relative_risk,
output_allele_info_files = "test_data/output_allele_info" ,
output_bed_files = "test_data/bed_file",
output_risk_genes_file = "test_data/riskgenes",
compact_mut_output = NA,
MPI = 1)
In [104]:
window_file = 'test_data/test_windows.txt'
mutrate_scaling_files = "test_data/test_simulated_windows_mutrate_with_div_score_scaling_file_for_test_DNM.txt"
sample_sizes = c(27)
gene_prior_file = "test_data/test_gene_prior.txt"
nonAS_noncoding_annotations = "test_data/test_non_allele_specific_annotation_1.bed"
AS_noncoding_annotations = NA
report_proportion = 100/18665
chunk_partition_num =1
node_n = 6
mutrate_ref_files = c("test_data/test_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_A.mutrate.bw",
"test_data/test_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_C.mutrate.bw",
"test_data/test_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_G.mutrate.bw",
"test_data/test_windows_extended_1bp_for_getting_base_level_mutrate.bed.fasta.tri.alt_T.mutrate.bw")
rr = c(1, 0 ,0.5)
output_allele_info_files = "test_data/output_allele_info"
output_bed_files = "test_data/bed_file"
output_risk_genes_file = "test_data/riskgenes"
compact_mut_output = NA
MPI = 1
In [109]:
base_ID = genename = mut_count = risk = NULL
prefix <- prefix <- as.integer((as.double(Sys.time())*1000+Sys.getpid()) %% 2^31) # prefix for temporary files that will be deleted at the end of the pipeline
prefix <- paste("tmp/", prefix, MPI, sep = "")
# make a tmp folder for tmp files
system("mkdir -p tmp")
#mut <- fread(mut_file, header = FALSE, sep = "\t", stringsAsFactors = FALSE)
windows <- fread(window_file, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
coverage <- windows
# the number of genomic windows in [mutrate_scaling] is less than the number of windows in [windows] because there are a few windows with mutration rate equal to 0, and thus removed.
for(i in 1:length(mutrate_scaling_files)){
mutrate_scaling <- fread(mutrate_scaling_files[i], header = TRUE, sep = "\t", stringsAsFactors = FALSE)
system(paste("echo \"Finished reading mutrate scaling file ", mutrate_scaling_files[i], ".\"", sep = ""))
system("date")
coverage <- coverage[mutrate_scaling, on = "site_index"]
coverage <- coverage[complete.cases(coverage)] # release memory
colnames(coverage)[length(colnames(coverage))] <- paste("scaling_factor_study_", i, sep = "")
rm(mutrate_scaling) # release memory
}
# get the piror probability of genes.
gene_prior <- fread(gene_prior_file, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
# merge gene prior info
coverage <- coverage[gene_prior, on = "genename"]
coverage <-coverage[complete.cases(coverage)]
In [110]:
head(coverage)
In [108]:
# select genes based on TADA prior probability and [report_proportion]
genes_for_report_with_prior <- gene_prior[order(gene_prior[,2],decreasing = TRUE),]
genes_for_report_with_prior <- genes_for_report_with_prior[1:floor(nrow(genes_for_report_with_prior)*report_proportion),]
# randomly assign risk genes based on prior probability
genes_for_report_with_prior$risk <- rbinom(nrow(genes_for_report_with_prior), 1, genes_for_report_with_prior$prior)
risk_genes <- genes_for_report_with_prior[risk == 1,]$genename
nonrisk_genes <- genes_for_report_with_prior[risk == 0,]$genename
coverage_risk <- coverage[is.element(genename, risk_genes)]
coverage_risk$risk <- 1
coverage_nonrisk <- coverage[is.element(genename, nonrisk_genes)]
coverage_nonrisk$risk <- 0
coverage <- rbind(coverage_risk, coverage_nonrisk)
# now need to extropolate window mutation file to base level mutation file
coverage <-coverage[,c("chr","start","end",paste("scaling_factor_study_", seq(1,length(mutrate_scaling_files)), sep = ""),"genename","risk"),with = FALSE]
coverage$ID <- paste(coverage$genename, coverage$start, sep = "_")
total_rows <- nrow(coverage)
interval <- floor(total_rows/chunk_partition_num)
data_bins <- c(rep(seq(1,chunk_partition_num), each = interval),rep(chunk_partition_num, total_rows -interval*chunk_partition_num))
# split into 20 different chunks, then for each chunk split by genes, then by feature configuration. If split by gene at the first level, implementation would take too long
coverage <- split(coverage, data_bins)
#funtion to expand windows to bases
window_expansion <- function(table_row){
start <- seq(as.numeric(table_row[2]),as.numeric(table_row[3])-1)
data.frame(table_row[1], start, start+1, paste(table_row["genename"],start,sep = "_"), table_row["ID"])
}
bed_extension <- function(table_row){
data.frame(chr = rep(table_row["chr"], table_row["mut_count"]), start = rep(table_row["start"], table_row["mut_count"]), end = rep(table_row["end"], table_row["mut_count"]))
}
options(warn=-1)
if(node_n != 1){
cl <- makeCluster(node_n)
}
# use parallel computing and rbindlist to increase the speed by thousands of times.
environment(window_expansion) <- .GlobalEnv
environment(bed_extension) <- .GlobalEnv
# get nonAS feature number
if(is.na(nonAS_noncoding_annotations[1])){
nonAS_feature_number <- 0
}else{
nonAS_feature_number <- length(nonAS_noncoding_annotations)
}
# get AS feature number
if(is.na(AS_noncoding_annotations[1])){
AS_feature_number <- 0
}else{
AS_feature_number <- length(AS_noncoding_annotations)
}
# get total feature number
feature_number = nonAS_feature_number + AS_feature_number
# function to get effective information of each element of partition_by_gene
# These information are those necessary to compute log-likelihood in the optimization function
partition_feature <- function(pbg){
# input is one element of the list of partition_by_gene
pbg_split <- split(pbg, pbg[,4:(4 + feature_number - 1)],drop = TRUE)
feature_combination_number <- length(pbg_split)
# this function below is different from the function used in dealing with dataset without reading by chunk. Here, prior is not incoporated at this step.
info_for_each_feature <- function(feature_set){
list(feature_vector = c(as.numeric(feature_set[1,4:(4 + feature_number - 1)])), sum_mut_rate_count = sum(feature_set$mut_count*log(feature_set$adjusted_base_mutrate)), sum_mut_rate = sum(feature_set$adjusted_base_mutrate), sum_mut_count = sum(feature_set$mut_count), log_fcount = sum(log(factorial(feature_set$mut_count))))
}
sapply(pbg_split, info_for_each_feature,simplify = FALSE)
}
alt_letters <- c("A","C","G","T")
mut_output_list <- list()
data_partition <- list()
for(m in 1:length(mutrate_scaling_files)){
mut_output_list[[m]] <- data.table(chr = character(), start = integer(), end = integer(), ref = character(), alt = character())
}
# here only deals with the noncoding parts that are within 10kb of TSSs of genes. Will deal with
for(i in 1:chunk_partition_num){
# split coverage_noncoding to 10 parts, and for each part, generate feature table (which will be used for )
if(node_n != 1){
coverage_noncoding_for_base_mutrate <- rbindlist(parApply(cl, coverage[[i]], 1, window_expansion))
}else{
coverage_noncoding_for_base_mutrate <- rbindlist(apply(coverage[[i]], 1, window_expansion))
}
system(paste("echo \"Finished partitioning base-level coordinates data at Round ", i, ".\"", sep = ""))
system("date")
colnames(coverage_noncoding_for_base_mutrate) <- c("chr","start","end","base_ID","ID")
coverage_noncoding_for_base_mutrate <- coverage_noncoding_for_base_mutrate[!duplicated(base_ID)]
coverage_noncoding_for_base_mutrate$start <- as.integer(coverage_noncoding_for_base_mutrate$start)
coverage_noncoding_for_base_mutrate$end <- as.integer(coverage_noncoding_for_base_mutrate$end)
# write out a bed file to get base-level mutation rates
fwrite(coverage_noncoding_for_base_mutrate[,1:4],paste(prefix, "_temp_for_mutrate.bed", sep = ""), col.names = FALSE, row.names = FALSE, sep = "\t", quote = FALSE)
# read in allele-specific base-level mutation rate
for(j in 1:length(mutrate_ref_files)){
command <- paste("../external_tools/bigWigAverageOverBed ", mutrate_ref_files[j], " ", paste(prefix, "_temp_for_mutrate.bed", sep = ""), " ", paste(prefix, "_temp_for_mutrate.bed.mutrate", sep = "" ), sep = "")
system(command)
command <- paste("awk {'print $1\"\t\"$4'} ", paste(prefix, "_temp_for_mutrate.bed.mutrate", sep = ""), " > ", paste(prefix, "_temp_for_mutrate.bed.mutrate.txt", sep = ""), sep = "")
system(command)
coverage_noncoding_base_mutrate <-fread(paste(prefix, "_temp_for_mutrate.bed.mutrate.txt", sep = ""), header = FALSE, stringsAsFactors = FALSE, sep = "\t")
colnames(coverage_noncoding_base_mutrate) <- c("base_ID",paste("base_mutrate_alt_", alt_letters[j], sep = ""))
coverage_noncoding_for_base_mutrate <- coverage_noncoding_for_base_mutrate[coverage_noncoding_base_mutrate, on = "base_ID"]
system(paste("echo \"Finished obtaining base-level mutation rate for alt allele ", alt_letters[j], ".\"", sep = ""))
system("date")
}
# read in non allele-specific epigenomic annotations
epi_ID = 1
if (!is.na(nonAS_noncoding_annotations)[1]){ # then epigenomic_marks must be a vector of epigenomic bed files that need to be compard with the mutation data
for(epi in nonAS_noncoding_annotations){
command <- paste("../external_tools/bedtools-2.17.0/bin/bedtools coverage -a ", epi, " -b ", paste(prefix, "_temp_for_mutrate.bed", sep = ""), " > ", paste(prefix,"_temp_for_mutrate_overlap_epi.bed", sep = ""),sep = "")
system(command)
base_in_epi <- fread(paste(prefix,"_temp_for_mutrate_overlap_epi.bed", sep = ""), header = FALSE, sep = "\t", stringsAsFactors = FALSE)
base_in_epi <- base_in_epi[,c("V4","V5"), with = FALSE]
colnames(base_in_epi) <- c("base_ID", paste("Anno",epi_ID, sep = "_"))
coverage_noncoding_for_base_mutrate <- coverage_noncoding_for_base_mutrate[base_in_epi, on = "base_ID"]
system(paste("echo \"Finished reading non-allele specific noncoding annotations ", epi_ID, ".\"", sep = ""))
system("date")
epi_ID <- epi_ID + 1
}
}
# read in allele-specific epigenomic annotations, now the base assignment for such noncoding annotations is based on the 50-bp windows.
# Should not be a big problem as distal introns are assigned based on refseq_ID for intron regions, should be consistent with assignment by spidex.
# though there is a small chance that a gene's distal intron, which is close to the promoter of another gene and within 10 kb from the TSSs of the latter gene, may be mis-assigned to the latter gene.
# To completely correct for this issue, need to allow epigenetic annotation to take its own gene assignment, which might be necessary in some situations under strict criteria, such as splicing annotaion.
if (!is.na(AS_noncoding_annotations)[1]){ # then epigenomic_marks must be a vector of epigenomic bed files that need to be compard with the mutation data
for(epi in AS_noncoding_annotations){
for(k in 1:length(epi)){
command <- paste("../external_tools/bedtools-2.17.0/bin/bedtools coverage -a ", epi[k], " -b ", paste(prefix, "_temp_for_mutrate.bed", sep = ""), " > ", paste(prefix,"_temp_for_mutrate_overlap_epi.bed", sep = ""),sep = "")
system(command)
base_in_epi <- fread(paste(prefix,"_temp_for_mutrate_overlap_epi.bed", sep = ""), header = FALSE, sep = "\t", stringsAsFactors = FALSE)
base_in_epi <- base_in_epi[,c("V4","V5"), with = FALSE]
colnames(base_in_epi) <- c("base_ID", paste("Anno",epi_ID, alt_letters[k], sep = "_"))
coverage_noncoding_for_base_mutrate <- coverage_noncoding_for_base_mutrate[base_in_epi, on = "base_ID"]
}
system(paste("echo \"Finished reading allele specific noncoding annotations ", epi_ID, ".\"", sep = ""))
system("date")
epi_ID <- epi_ID + 1
}
}
# now for each sample size, generate mut data.
for(m in 1:length(mutrate_scaling_files)){
mut_output_temp <- data.table(chr = character(), start = integer(), end = integer(), ref = character(), alt = character())
for(letter in alt_letters){
if(!is.na(nonAS_noncoding_annotations)[1] & !is.na(AS_noncoding_annotations)[1]){
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate[, c("base_ID", "ID",paste("base_mutrate_alt_", letter, sep = ""), paste("Anno_", seq(1, nonAS_feature_number), sep = ""), paste("Anno_", seq(nonAS_feature_number +1, nonAS_feature_number + AS_feature_number), "_", letter ,sep = ""), "chr", "start", "end"), with = FALSE]
}else if(!is.na(nonAS_noncoding_annotations)[1] & is.na(AS_noncoding_annotations)[1]){
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate[, c("base_ID", "ID",paste("base_mutrate_alt_", letter, sep = ""), paste("Anno_", seq(1, nonAS_feature_number), sep = ""), "chr", "start", "end"), with = FALSE]
}else if(is.na(nonAS_noncoding_annotations)[1] & !is.na(AS_noncoding_annotations)[1]){
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate[, c("base_ID", "ID",paste("base_mutrate_alt_", letter, sep = ""), paste("Anno_", seq(nonAS_feature_number +1, nonAS_feature_number + AS_feature_number), "_", letter ,sep = ""), "chr", "start", "end"), with = FALSE]
}
# a very important step here is removing bases that have adjusted_mutrate_base 0. This happens when the allele of the mutrate we are using is just the reference allele.
# By doing this, we make the computation of likelihood valid, also we automatically removed bases with nonAS annotations but with mutant allele the same with ref allele at this step
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate_temp[coverage_noncoding_for_base_mutrate_temp[[paste("base_mutrate_alt_", letter, sep = "")]] !=0]
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate_temp[coverage[[i]][,c(paste("scaling_factor_study_", m, sep = ""), "genename","ID","risk"), with = FALSE],on = "ID"]
# removing sites with allele-specific mutrate == 0 will result in incomplete rows after merging.
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate_temp[complete.cases(coverage_noncoding_for_base_mutrate_temp)]
# derive the adjusted mutation rate for bases of risk genes and nonrisk genes.
# first multiplies sample size and scaling factor, adjusted_base_mutrate will be used for generating data in a compact format.
coverage_noncoding_for_base_mutrate_temp$adjusted_base_mutrate <- coverage_noncoding_for_base_mutrate_temp[, paste("base_mutrate_alt_", letter, sep = ""), with = FALSE] * coverage_noncoding_for_base_mutrate_temp[, paste("scaling_factor_study_", m, sep = ""), with = FALSE] * 2 * sample_sizes[m]
# Then multiplies effect size after taking exponential, to generate a mutrate that has taken into account RRs, will be used for simulating mutations
#
coverage_noncoding_for_base_mutrate_temp$adjusted_base_mutrate_rr <- exp(as.matrix(coverage_noncoding_for_base_mutrate_temp[,4 : (3 + feature_number)]) %*% rr * coverage_noncoding_for_base_mutrate_temp$risk) * coverage_noncoding_for_base_mutrate_temp$adjusted_base_mutrate
# generate random allele-specific mutations based on adjusted mutation rate.
coverage_noncoding_for_base_mutrate_temp$mut_count <- rpois(nrow(coverage_noncoding_for_base_mutrate_temp), coverage_noncoding_for_base_mutrate_temp$adjusted_base_mutrate_rr)
if(!is.na(compact_mut_output)){
# have to collpase data at this point, otherwise, the I will run out of RAM if process 1000 genes every time.
anno_count <- rep(0, nrow(coverage_noncoding_for_base_mutrate_temp))
for(p in 1:feature_number){
anno_count <- anno_count + coverage_noncoding_for_base_mutrate_temp[[(4 + p -1)]]
}
# remove rows(bases) that don't have any non-coding features, this could save a lot of RAM, so I could use smaller partition number which would greatly accelerate speed when read in data for all genes.
coverage_noncoding_for_base_mutrate_temp1 <- coverage_noncoding_for_base_mutrate_temp[anno_count >0,]
# first partition by gene for the current chunk
coverage_noncoding_for_base_mutrate_temp1<- split(coverage_noncoding_for_base_mutrate_temp1, coverage_noncoding_for_base_mutrate_temp1$genename)
# then partition by feature configuration for each gene in the current chunk
coverage_noncoding_for_base_mutrate_temp1 <- sapply(coverage_noncoding_for_base_mutrate_temp1, partition_feature, simplify = FALSE)
# add compact data
data_partition <- append(data_partition, coverage_noncoding_for_base_mutrate_temp1)
rm(coverage_noncoding_for_base_mutrate_temp1) # release memory
system(paste("echo \"Finished read in mutation data and make them into the compact format for Study ", m, " and allele ", letter, ".\"", sep = ""))
system("date")
}
# then generate bed file format, this is a pseudobed file as I will put N at the 4th column as the ref allele is not relevant and mutant allele at the last column.
# remove rows(bases) that don't have any DNMs
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate_temp[mut_count >0,]
if(node_n != 1){
mut_bed <- rbindlist(parApply(cl, coverage_noncoding_for_base_mutrate_temp, 1, bed_extension))
}else{
mut_bed <- rbindlist(apply(coverage_noncoding_for_base_mutrate_temp, 1, bed_extension))
}
mut_bed$ref <- "N"
mut_bed$alt <- letter
mut_bed$start <- as.integer(as.character(mut_bed$start))
mut_bed$end <- as.integer(as.character(mut_bed$end))
mut_output_temp <- rbind(mut_output_temp, mut_bed)
}
mut_output_list[[m]] <- rbind(mut_output_list[[m]], mut_output_temp)
}
}
if(node_n != 1){
stopCluster(cl)
}
options(warn = 0)
# remove temporary files
system(paste("rm ", prefix, "_temp*", sep = ""))
print(paste("echo \"Temp files cleaned and mutation generation finished finished!\""))
system("date")
# now write mutation file and mutation with allele info file
for(m in 1:length(mutrate_scaling_files)){
fwrite(mut_output_list[[m]], output_allele_info_files[m], col.names = FALSE, row.names = FALSE, sep = "\t", quote = FALSE)
fwrite(mut_output_list[[m]][,1:3], output_bed_files[m], col.names = FALSE, row.names = FALSE, sep = "\t", quote = FALSE)
}
write.table(risk_genes, output_risk_genes_file, col.names = FALSE, row.names = FALSE, sep = "\t", quote = FALSE)
if(!is.na(compact_mut_output)){
saveRDS(data_partition, compact_mut_output)
}
In [107]:
library(parallel)
In [ ]:
TADA_A_DNM_generatora <- function(window_file,
mutrate_scaling_files,
sample_sizes,
gene_prior_file,
nonAS_noncoding_annotations,
AS_noncoding_annotations = NA,
report_proportion,
chunk_partition_num = 1,
node_n = 6,
mutrate_ref_files,
rr,
output_allele_info_files,
output_bed_files,
output_risk_genes_file,
compact_mut_output = NA,
MPI = 1){
base_ID = genename = mut_count = risk = NULL
prefix <- prefix <- as.integer((as.double(Sys.time())*1000+Sys.getpid()) %% 2^31) # prefix for temporary files that will be deleted at the end of the pipeline
prefix <- paste("tmp/", prefix, MPI, sep = "")
# make a tmp folder for tmp files
system("mkdir -p tmp")
#mut <- fread(mut_file, header = FALSE, sep = "\t", stringsAsFactors = FALSE)
windows <- fread(window_file, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
coverage <- windows
# the number of genomic windows in [mutrate_scaling] is less than the number of windows in [windows] because there are a few windows with mutration rate equal to 0, and thus removed.
for(i in 1:length(mutrate_scaling_files)){
mutrate_scaling <- fread(mutrate_scaling_files[i], header = TRUE, sep = "\t", stringsAsFactors = FALSE)
system(paste("echo \"Finished reading mutrate scaling file ", mutrate_scaling_files[i], ".\"", sep = ""))
system("date")
coverage <- coverage[mutrate_scaling, on = "site_index"]
coverage <- coverage[complete.cases(coverage)] # release memory
colnames(coverage)[length(colnames(coverage))] <- paste("scaling_factor_study_", i, sep = "")
rm(mutrate_scaling) # release memory
}
# get the piror probability of genes.
gene_prior <- fread(gene_prior_file, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
# merge gene prior info
coverage <- coverage[gene_prior, on = "genename"]
coverage <-coverage[complete.cases(coverage)]
# select genes based on TADA prior probability and [report_proportion]
genes_for_report_with_prior <- gene_prior[order(gene_prior[,2],decreasing = TRUE),]
genes_for_report_with_prior <- genes_for_report_with_prior[1:floor(nrow(genes_for_report_with_prior)*report_proportion),]
# randomly assign risk genes based on prior probability
genes_for_report_with_prior$risk <- rbinom(nrow(genes_for_report_with_prior), 1, genes_for_report_with_prior$prior)
risk_genes <- genes_for_report_with_prior[risk == 1,]$genename
nonrisk_genes <- genes_for_report_with_prior[risk == 0,]$genename
coverage_risk <- coverage[is.element(genename, risk_genes)]
coverage_risk$risk <- 1
coverage_nonrisk <- coverage[is.element(genename, nonrisk_genes)]
coverage_nonrisk$risk <- 0
coverage <- rbind(coverage_risk, coverage_nonrisk)
# now need to extropolate window mutation file to base level mutation file
coverage <-coverage[,c("chr","start","end",paste("scaling_factor_study_", seq(1,length(mutrate_scaling_files)), sep = ""),"genename","risk"),with = FALSE]
coverage$ID <- paste(coverage$genename, coverage$start, sep = "_")
total_rows <- nrow(coverage)
interval <- floor(total_rows/chunk_partition_num)
data_bins <- c(rep(seq(1,chunk_partition_num), each = interval),rep(chunk_partition_num, total_rows -interval*chunk_partition_num))
# split into 20 different chunks, then for each chunk split by genes, then by feature configuration. If split by gene at the first level, implementation would take too long
coverage <- split(coverage, data_bins)
#funtion to expand windows to bases
window_expansion <- function(table_row){
start <- seq(as.numeric(table_row[2]),as.numeric(table_row[3])-1)
data.frame(table_row[1], start, start+1, paste(table_row["genename"],start,sep = "_"), table_row["ID"])
}
bed_extension <- function(table_row){
data.frame(chr = rep(table_row["chr"], table_row["mut_count"]), start = rep(table_row["start"], table_row["mut_count"]), end = rep(table_row["end"], table_row["mut_count"]))
}
options(warn=-1)
if(node_n != 1){
cl <- makeCluster(node_n)
}
# use parallel computing and rbindlist to increase the speed by thousands of times.
environment(window_expansion) <- .GlobalEnv
environment(bed_extension) <- .GlobalEnv
# get nonAS feature number
if(is.na(nonAS_noncoding_annotations[1])){
nonAS_feature_number <- 0
}else{
nonAS_feature_number <- length(nonAS_noncoding_annotations)
}
# get AS feature number
if(is.na(AS_noncoding_annotations[1])){
AS_feature_number <- 0
}else{
AS_feature_number <- length(AS_noncoding_annotations)
}
# get total feature number
feature_number = nonAS_feature_number + AS_feature_number
# function to get effective information of each element of partition_by_gene
# These information are those necessary to compute log-likelihood in the optimization function
partition_feature <- function(pbg){
# input is one element of the list of partition_by_gene
pbg_split <- split(pbg, pbg[,4:(4 + feature_number - 1)],drop = TRUE)
feature_combination_number <- length(pbg_split)
# this function below is different from the function used in dealing with dataset without reading by chunk. Here, prior is not incoporated at this step.
info_for_each_feature <- function(feature_set){
list(feature_vector = c(as.numeric(feature_set[1,4:(4 + feature_number - 1)])), sum_mut_rate_count = sum(feature_set$mut_count*log(feature_set$adjusted_base_mutrate)), sum_mut_rate = sum(feature_set$adjusted_base_mutrate), sum_mut_count = sum(feature_set$mut_count), log_fcount = sum(log(factorial(feature_set$mut_count))))
}
sapply(pbg_split, info_for_each_feature,simplify = FALSE)
}
alt_letters <- c("A","C","G","T")
mut_output_list <- list()
data_partition <- list()
for(m in 1:length(mutrate_scaling_files)){
mut_output_list[[m]] <- data.table(chr = character(), start = integer(), end = integer(), ref = character(), alt = character())
}
# here only deals with the noncoding parts that are within 10kb of TSSs of genes. Will deal with
for(i in 1:chunk_partition_num){
# split coverage_noncoding to 10 parts, and for each part, generate feature table (which will be used for )
if(node_n != 1){
coverage_noncoding_for_base_mutrate <- rbindlist(parApply(cl, coverage[[i]], 1, window_expansion))
}else{
coverage_noncoding_for_base_mutrate <- rbindlist(apply(coverage[[i]], 1, window_expansion))
}
system(paste("echo \"Finished partitioning base-level coordinates data at Round ", i, ".\"", sep = ""))
system("date")
colnames(coverage_noncoding_for_base_mutrate) <- c("chr","start","end","base_ID","ID")
coverage_noncoding_for_base_mutrate <- coverage_noncoding_for_base_mutrate[!duplicated(base_ID)]
coverage_noncoding_for_base_mutrate$start <- as.integer(coverage_noncoding_for_base_mutrate$start)
coverage_noncoding_for_base_mutrate$end <- as.integer(coverage_noncoding_for_base_mutrate$end)
# write out a bed file to get base-level mutation rates
fwrite(coverage_noncoding_for_base_mutrate[,1:4],paste(prefix, "_temp_for_mutrate.bed", sep = ""), col.names = FALSE, row.names = FALSE, sep = "\t", quote = FALSE)
# read in allele-specific base-level mutation rate
for(j in 1:length(mutrate_ref_files)){
command <- paste("../external_tools/bigWigAverageOverBed ", mutrate_ref_files[j], " ", paste(prefix, "_temp_for_mutrate.bed", sep = ""), " ", paste(prefix, "_temp_for_mutrate.bed.mutrate", sep = "" ), sep = "")
system(command)
command <- paste("awk {'print $1\"\t\"$4'} ", paste(prefix, "_temp_for_mutrate.bed.mutrate", sep = ""), " > ", paste(prefix, "_temp_for_mutrate.bed.mutrate.txt", sep = ""), sep = "")
system(command)
coverage_noncoding_base_mutrate <-fread(paste(prefix, "_temp_for_mutrate.bed.mutrate.txt", sep = ""), header = FALSE, stringsAsFactors = FALSE, sep = "\t")
colnames(coverage_noncoding_base_mutrate) <- c("base_ID",paste("base_mutrate_alt_", alt_letters[j], sep = ""))
coverage_noncoding_for_base_mutrate <- coverage_noncoding_for_base_mutrate[coverage_noncoding_base_mutrate, on = "base_ID"]
system(paste("echo \"Finished obtaining base-level mutation rate for alt allele ", alt_letters[j], ".\"", sep = ""))
system("date")
}
# read in non allele-specific epigenomic annotations
epi_ID = 1
if (!is.na(nonAS_noncoding_annotations)[1]){ # then epigenomic_marks must be a vector of epigenomic bed files that need to be compard with the mutation data
for(epi in nonAS_noncoding_annotations){
command <- paste("../external_tools/bedtools-2.17.0/bin/bedtools coverage -a ", epi, " -b ", paste(prefix, "_temp_for_mutrate.bed", sep = ""), " > ", paste(prefix,"_temp_for_mutrate_overlap_epi.bed", sep = ""),sep = "")
system(command)
base_in_epi <- fread(paste(prefix,"_temp_for_mutrate_overlap_epi.bed", sep = ""), header = FALSE, sep = "\t", stringsAsFactors = FALSE)
base_in_epi <- base_in_epi[,c("V4","V5"), with = FALSE]
colnames(base_in_epi) <- c("base_ID", paste("Anno",epi_ID, sep = "_"))
coverage_noncoding_for_base_mutrate <- coverage_noncoding_for_base_mutrate[base_in_epi, on = "base_ID"]
system(paste("echo \"Finished reading non-allele specific noncoding annotations ", epi_ID, ".\"", sep = ""))
system("date")
epi_ID <- epi_ID + 1
}
}
# read in allele-specific epigenomic annotations, now the base assignment for such noncoding annotations is based on the 50-bp windows.
# Should not be a big problem as distal introns are assigned based on refseq_ID for intron regions, should be consistent with assignment by spidex.
# though there is a small chance that a gene's distal intron, which is close to the promoter of another gene and within 10 kb from the TSSs of the latter gene, may be mis-assigned to the latter gene.
# To completely correct for this issue, need to allow epigenetic annotation to take its own gene assignment, which might be necessary in some situations under strict criteria, such as splicing annotaion.
if (!is.na(AS_noncoding_annotations)[1]){ # then epigenomic_marks must be a vector of epigenomic bed files that need to be compard with the mutation data
for(epi in AS_noncoding_annotations){
for(k in 1:length(epi)){
command <- paste("../external_tools/bedtools-2.17.0/bin/bedtools coverage -a ", epi[k], " -b ", paste(prefix, "_temp_for_mutrate.bed", sep = ""), " > ", paste(prefix,"_temp_for_mutrate_overlap_epi.bed", sep = ""),sep = "")
system(command)
base_in_epi <- fread(paste(prefix,"_temp_for_mutrate_overlap_epi.bed", sep = ""), header = FALSE, sep = "\t", stringsAsFactors = FALSE)
base_in_epi <- base_in_epi[,c("V4","V5"), with = FALSE]
colnames(base_in_epi) <- c("base_ID", paste("Anno",epi_ID, alt_letters[k], sep = "_"))
coverage_noncoding_for_base_mutrate <- coverage_noncoding_for_base_mutrate[base_in_epi, on = "base_ID"]
}
system(paste("echo \"Finished reading allele specific noncoding annotations ", epi_ID, ".\"", sep = ""))
system("date")
epi_ID <- epi_ID + 1
}
}
# now for each sample size, generate mut data.
for(m in 1:length(mutrate_scaling_files)){
mut_output_temp <- data.table(chr = character(), start = integer(), end = integer(), ref = character(), alt = character())
for(letter in alt_letters){
if(!is.na(nonAS_noncoding_annotations)[1] & !is.na(AS_noncoding_annotations)[1]){
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate[, c("base_ID", "ID",paste("base_mutrate_alt_", letter, sep = ""), paste("Anno_", seq(1, nonAS_feature_number), sep = ""), paste("Anno_", seq(nonAS_feature_number +1, nonAS_feature_number + AS_feature_number), "_", letter ,sep = ""), "chr", "start", "end"), with = FALSE]
}else if(!is.na(nonAS_noncoding_annotations)[1] & is.na(AS_noncoding_annotations)[1]){
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate[, c("base_ID", "ID",paste("base_mutrate_alt_", letter, sep = ""), paste("Anno_", seq(1, nonAS_feature_number), sep = ""), "chr", "start", "end"), with = FALSE]
}else if(is.na(nonAS_noncoding_annotations)[1] & !is.na(AS_noncoding_annotations)[1]){
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate[, c("base_ID", "ID",paste("base_mutrate_alt_", letter, sep = ""), paste("Anno_", seq(nonAS_feature_number +1, nonAS_feature_number + AS_feature_number), "_", letter ,sep = ""), "chr", "start", "end"), with = FALSE]
}
# a very important step here is removing bases that have adjusted_mutrate_base 0. This happens when the allele of the mutrate we are using is just the reference allele.
# By doing this, we make the computation of likelihood valid, also we automatically removed bases with nonAS annotations but with mutant allele the same with ref allele at this step
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate_temp[coverage_noncoding_for_base_mutrate_temp[[paste("base_mutrate_alt_", letter, sep = "")]] !=0]
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate_temp[coverage[[i]][,c(paste("scaling_factor_study_", m, sep = ""), "genename","ID","risk"), with = FALSE],on = "ID"]
# removing sites with allele-specific mutrate == 0 will result in incomplete rows after merging.
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate_temp[complete.cases(coverage_noncoding_for_base_mutrate_temp)]
# derive the adjusted mutation rate for bases of risk genes and nonrisk genes.
# first multiplies sample size and scaling factor, adjusted_base_mutrate will be used for generating data in a compact format.
coverage_noncoding_for_base_mutrate_temp$adjusted_base_mutrate <- coverage_noncoding_for_base_mutrate_temp[, paste("base_mutrate_alt_", letter, sep = ""), with = FALSE] * coverage_noncoding_for_base_mutrate_temp[, paste("scaling_factor_study_", m, sep = ""), with = FALSE] * 2 * sample_sizes[m]
# Then multiplies effect size after taking exponential, to generate a mutrate that has taken into account RRs, will be used for simulating mutations
#
coverage_noncoding_for_base_mutrate_temp$adjusted_base_mutrate_rr <- exp(as.matrix(coverage_noncoding_for_base_mutrate_temp[,4 : (3 + feature_number)]) %*% rr * coverage_noncoding_for_base_mutrate_temp$risk) * coverage_noncoding_for_base_mutrate_temp$adjusted_base_mutrate
# generate random allele-specific mutations based on adjusted mutation rate.
coverage_noncoding_for_base_mutrate_temp$mut_count <- rpois(nrow(coverage_noncoding_for_base_mutrate_temp), coverage_noncoding_for_base_mutrate_temp$adjusted_base_mutrate_rr)
if(!is.na(compact_mut_output)){
# have to collpase data at this point, otherwise, the I will run out of RAM if process 1000 genes every time.
anno_count <- rep(0, nrow(coverage_noncoding_for_base_mutrate_temp))
for(p in 1:feature_number){
anno_count <- anno_count + coverage_noncoding_for_base_mutrate_temp[[(4 + p -1)]]
}
# remove rows(bases) that don't have any non-coding features, this could save a lot of RAM, so I could use smaller partition number which would greatly accelerate speed when read in data for all genes.
coverage_noncoding_for_base_mutrate_temp1 <- coverage_noncoding_for_base_mutrate_temp[anno_count >0,]
# first partition by gene for the current chunk
coverage_noncoding_for_base_mutrate_temp1<- split(coverage_noncoding_for_base_mutrate_temp1, coverage_noncoding_for_base_mutrate_temp1$genename)
# then partition by feature configuration for each gene in the current chunk
coverage_noncoding_for_base_mutrate_temp1 <- sapply(coverage_noncoding_for_base_mutrate_temp1, partition_feature, simplify = FALSE)
# add compact data
data_partition <- append(data_partition, coverage_noncoding_for_base_mutrate_temp1)
rm(coverage_noncoding_for_base_mutrate_temp1) # release memory
system(paste("echo \"Finished read in mutation data and make them into the compact format for Study ", m, " and allele ", letter, ".\"", sep = ""))
system("date")
}
# then generate bed file format, this is a pseudobed file as I will put N at the 4th column as the ref allele is not relevant and mutant allele at the last column.
# remove rows(bases) that don't have any DNMs
coverage_noncoding_for_base_mutrate_temp <- coverage_noncoding_for_base_mutrate_temp[mut_count >0,]
if(node_n != 1){
mut_bed <- rbindlist(parApply(cl, coverage_noncoding_for_base_mutrate_temp, 1, bed_extension))
}else{
mut_bed <- rbindlist(apply(coverage_noncoding_for_base_mutrate_temp, 1, bed_extension))
}
mut_bed$ref <- "N"
mut_bed$alt <- letter
mut_bed$start <- as.integer(as.character(mut_bed$start))
mut_bed$end <- as.integer(as.character(mut_bed$end))
mut_output_temp <- rbind(mut_output_temp, mut_bed)
}
mut_output_list[[m]] <- rbind(mut_output_list[[m]], mut_output_temp)
}
}
if(node_n != 1){
stopCluster(cl)
}
options(warn = 0)
# remove temporary files
system(paste("rm ", prefix, "_temp*", sep = ""))
print(paste("echo \"Temp files cleaned and mutation generation finished finished!\""))
system("date")
# now write mutation file and mutation with allele info file
for(m in 1:length(mutrate_scaling_files)){
fwrite(mut_output_list[[m]], output_allele_info_files[m], col.names = FALSE, row.names = FALSE, sep = "\t", quote = FALSE)
fwrite(mut_output_list[[m]][,1:3], output_bed_files[m], col.names = FALSE, row.names = FALSE, sep = "\t", quote = FALSE)
}
write.table(risk_genes, output_risk_genes_file, col.names = FALSE, row.names = FALSE, sep = "\t", quote = FALSE)
if(!is.na(compact_mut_output)){
saveRDS(data_partition, compact_mut_output)
}
}
retrieve_mutation_info¶
In [ ]:
retrieve_mutation_info(data_partition,
genename,
noncoding_annotation)