Merge pull request #123 from griffithlab/filter_and_BH

Filter and bh

Author: yang-yangfeng

scripts/compare_junctions_hist.R

@@ -139,7 +139,7 @@ a <- function(x, all_cse_identify_data){
 
   # Aggregate data across junction-samples
   if (nrow(variant_junctions_data)[[1]] > 0){
-    non_variant_junctions_aggr = non_variant_junctions_data[, list(mean_norm_score_non=mean(norm_score),sd_norm_score_non=sd(norm_score),norm_scores_non=sd(norm_score),total_score_non=sum(score)), 
+    non_variant_junctions_aggr = non_variant_junctions_data[, list(mean_norm_score_non=mean(norm_score),sd_norm_score_non=sd(norm_score),norm_scores_non=list(norm_score),total_score_non=sum(score)), 
                                                     by=list(chrom,start,end,strand,anchor,variant_info,info)]
   } else {
     return(data.table())

scripts/filter_and_BH.R

@@ -0,0 +1,41 @@
+# filter_and_BH.R
+library(data.table)
+library(stats)
+
+debug = F
+
+if (debug){
+  tag = "_default"
+} else {
+  # get options tag
+  argc = length(commandArgs())
+  tag = paste("_", commandArgs(trailingOnly = F)[argc], sep="")
+  
+  if ( substr(tag, 2, 3) == "--"){
+    stop("Please specify an option tag (e.g. \"default\", \"i20e5\")")
+  }
+}
+
+
+read_file=paste("compare_junctions/hist/", "junction_pvalues", tag, ".tsv", sep="")
+regtools_data = unique(data.table::fread(file=read_file, sep = '\t', header = TRUE, stringsAsFactors = FALSE))
+regtools_data_filtered = regtools_data[(regtools_data$total_score_variant > 5 & 
+											regtools_data$pvalue >= 0 & 
+											(regtools_data$anchor == "D" | 
+												regtools_data$anchor == "A" | 
+												regtools_data$anchor == "NDA"))]
+
+p = regtools_data_filtered$pvalue
+adjusted_p = p.adjust(p, method = "BH")
+regtools_data_filtered$adjusted_p = adjusted_p
+regtools_data_filtered_sorted = regtools_data_filtered[order(adjusted_p)]
+
+write_file = paste("compare_junctions/hist/", "junction_pvalues_filtered_BH", tag, ".tsv", sep="")
+write.table(regtools_data_filtered_sorted, file=write_file, quote=FALSE, sep='\t', row.names = FALSE)
+
+threshold = 0.05
+is_significant = regtools_data_filtered_sorted$adjusted_p < 0.05
+regtools_data_significant_filtered_sorted = regtools_data_filtered_sorted[is_significant] 
+
+write_file = paste("compare_junctions/hist/", "junction_pvalues_significant_",threshold,"_filtered_BH", tag, ".tsv", sep="")
+write.table(regtools_data_significant_filtered_sorted, file=write_file, quote=FALSE, sep='\t', row.names = FALSE)