Merged ddpinto's changes

leafcutter/R/differential_splicing.R

@@ -4,11 +4,18 @@
 #'
 #' @param results From \code{\link{differential_splicing}}
 #' @return Data.frame with columns status, log likelihood ratio, degrees of freedom, p-value
+#' @importFrom dplyr bind_rows
 #' @export
 cluster_results_table=function(results) {
-  as.data.frame(cbind(cluster=names(results), foreach(res=results, .combine=rbind) %do% 
-{ if ( is.character(res) | ("error" %in% class(res)) ) data.frame(status=as.character(res), loglr=NA, df=NA, p=NA) else 
-  data.frame(status="Success", loglr=res$loglr, df=res$df, p=res$lrtp) } ) )
+  cluster_table=as.data.frame(cbind(cluster=names(results), foreach(res=results, .combine=bind_rows) %do% { 
+    if ( is.character(res) | ("error" %in% class(res)) ) 
+      data.frame(status=as.character(res), loglr=NA, df=NA, p=NA) else 
+      data.frame(status="Success", loglr=res$loglr, df=res$df, p=res$lrtp) 
+    } ) )
+  # Make sure we strip newlines from any reported errors
+  cluster_table$status   = gsub("\n", "", cluster_table$status) 
+  cluster_table$p.adjust = p.adjust(cluster_table$p, method="fdr")
+  cluster_table
 }
 
 #' Convert K-1 representation of parameters to real

scripts/leafcutter_ds.R

@@ -10,7 +10,7 @@ arguments <- parse_args(OptionParser(usage = "%prog [options] counts_file groups
   make_option(c("-c","--min_coverage"), default=20, help="Require min_samples_per_group samples in each group to have at least this many reads [default %default]"), 
   make_option(c("-t","--timeout"), default=30, help="Maximum time (in seconds) allowed for a single optimization run [default %default]"),
   make_option(c("-p","--num_threads"), default=1, help="Number of threads to use [default %default]"),
-    make_option(c("-e","--exon_file"), default=NULL, help="File defining known exons, example in data/gencode19_exons.txt.gz. Columns should be chr, start, end, strand, gene_name. Optional, only just to label the clusters."))),
+  make_option(c("-e","--exon_file"), default=NULL, help="File defining known exons, example in data/gencode19_exons.txt.gz. Columns should be chr, start, end, strand, gene_name. Optional, only just to label the clusters."))),
   positional_arguments = 2)
 
 opt=arguments$opt
@@ -65,28 +65,31 @@ results <- differential_splicing(counts, numeric_x, confounders=confounders, max
 
 cat("Saving results...\n")
 
-cluster_table=cluster_results_table(results)
-cluster_table$cluster=add_chr(cluster_table$cluster)
+# Make cluster table
+cluster_table          = cluster_results_table(results)
+cluster_table$cluster  = add_chr(cluster_table$cluster)
 
+# Add gene names to clusters if an exon file is available
 if (!is.null(opt$exon_file)) {
   cat("Loading exons from",opt$exon_file,"\n")
   if (file.exists(opt$exon_file)) {
      tryCatch( {
-          exons_table=read.table(opt$exon_file, header=T, stringsAsFactors = F)
-          intron_meta=get_intron_meta(rownames(counts))
-          exons_table$chr=add_chr(exons_table$chr)
-          intron_meta$chr=add_chr(intron_meta$chr)
-          clu_gene_map=map_clusters_to_genes(intron_meta, exons_table)
-          rownames(clu_gene_map)=clu_gene_map$clu
-          cluster_table$genes=clu_gene_map[ cluster_table$clu, "genes" ]
+          exons_table     = read.table(opt$exon_file, header=T, stringsAsFactors = F)
+          intron_meta     = get_intron_meta(rownames(counts))
+          exons_table$chr = add_chr(exons_table$chr)
+          intron_meta$chr = add_chr(intron_meta$chr)
+          clu_gene_map    = map_clusters_to_genes(intron_meta, exons_table)
+          cluster_table   = merge(cluster_table, clu_gene_map, by.x="cluster", by.y="clu", all.x=TRUE)
      }, error=function(err) warning(as.character(err)) ) # ignore errors here
   } else warning("File ",opt$exon_file," does not exist")
 } else cat("No exon_file provided.\n")
-
 write.table( cluster_table, paste0(opt$output_prefix,"_cluster_significance.txt"), quote=F, sep="\t", row.names = F)
 
-effect_size_table=leaf_cutter_effect_sizes(results)
-colnames(effect_size_table)[3:4]=group_names
+# Write effect size table
+effect_size_table                = leaf_cutter_effect_sizes(results)
+colnames(effect_size_table)[3:4] = group_names
 write.table( effect_size_table, paste0(opt$output_prefix,"_effect_sizes.txt"), quote=F, col.names = T, row.names = F, sep="\t")
 
+# Save RData image
+# save.image(paste0(opt$output_prefix,"_cluster_significance.RData"));
 cat("All done, exiting\n")