Merge pull request #410 from gagneurlab/summary_edits

Summary edits
gagneurlab · Dec 6, 2022 · 7fc6c43 · 7fc6c43
2 parents ffc786d + 0dce3ec
commit 7fc6c43
Show file tree

Hide file tree

Showing 3 changed files with 35 additions and 42 deletions.
diff --git a/drop/modules/aberrant-expression-pipeline/Counting/Summary.R b/drop/modules/aberrant-expression-pipeline/Counting/Summary.R
@@ -164,15 +164,16 @@ exp_genes_cols <- c(Rank = "expressedGenesRank",`Expressed\ngenes` = "expressedG
                     `Intersection of\nexpressed genes` = "intersectionExpressedGenes", 
                     `Genes passed\nfiltering` = "passedFilterGenes", `Is External` = "isExternal")
 
-expressed_genes <- as.data.table(colData(ods)[,exp_genes_cols])
-colnames(expressed_genes) <- names(exp_genes_cols)
+expressed_genes <- as.data.table(colData(ods)[,exp_genes_cols], keep.rownames = TRUE)
+colnames(expressed_genes) <- c('RNA_ID', names(exp_genes_cols))
 
 #+ expressedGenes, fig.height=6, fig.width=8
 plotExpressedGenes(ods) + 
   theme_cowplot() +
   background_grid(major = "y") +
-  geom_point(data =melt(expressed_genes,id.vars = c("Rank","Is External")),
-             aes(x = Rank, y = value, col = variable, shape = `Is External`),show.legend = has_external)
+  geom_point(data = melt(expressed_genes, id.vars = c("RNA_ID", "Rank", "Is External")),
+             aes(Rank, value, col = variable, shape = `Is External`), 
+             show.legend = has_external)
 
 if(has_external){
     DT::datatable(expressed_genes[order(Rank)],rownames = F)
@@ -238,5 +239,8 @@ if(isEmpty(sex_idx)){
 
     DT::datatable(sex_dt[match_sex == F], caption = 'Sex mismatches')
 
+    # Write table
+    fwrite(sex_dt, gsub('ods_unfitted.Rds', 'xist_uty.tsv', snakemake@input$ods), 
+           sep = '\t', quote = F)
   }
 }
diff --git a/drop/modules/aberrant-splicing-pipeline/FRASER/07_extract_results_FraseR.R b/drop/modules/aberrant-splicing-pipeline/FRASER/07_extract_results_FraseR.R
@@ -81,7 +81,7 @@ if(nrow(res_junc_dt) > 0){
 
   # add colData to the results
   res_junc_dt <- merge(res_junc_dt, as.data.table(colData(fds)), by = "sampleID")
-  res_junc_dt[, c("bamFile", "pairedEnd", "STRAND") := NULL]
+  res_junc_dt[, c("bamFile", "pairedEnd", "STRAND", "COUNT_MODE", "COUNT_OVERLAPS") := NULL]
 } else{
   warning("The aberrant splicing pipeline gave 0 results for the ", dataset, " dataset.")
 }
@@ -90,7 +90,7 @@ if(nrow(res_junc_dt) > 0){
 if(length(res_junc) > 0){
   res_genes_dt <- resultsByGenes(res_junc) %>% as.data.table
   res_genes_dt <- merge(res_genes_dt, as.data.table(colData(fds)), by = "sampleID")
-  res_genes_dt[, c("bamFile", "pairedEnd", "STRAND") := NULL]
+  res_genes_dt[, c("bamFile", "pairedEnd", "STRAND", "COUNT_MODE", "COUNT_OVERLAPS") := NULL]
 
   # add HPO overlap information
   sa <- fread(snakemake@config$sampleAnnotation, 

diff --git a/drop/modules/aberrant-splicing-pipeline/FRASER/Summary.R b/drop/modules/aberrant-splicing-pipeline/FRASER/Summary.R
@@ -68,47 +68,36 @@ plotAberrantPerSample(fds, padjCutoff = padj_cutoff, zScoreCutoff = zScore_cutof
 topN <- 30000
 topJ <- 10000
 anno_color_scheme <- brewer.pal(n = 3, name = 'Dark2')[1:2]
+
 for(type in psiTypes){
-  before <- plotCountCorHeatmap(
-    object=fds,
-    type = type,
-    logit = TRUE,
-    topN = topN,
-    topJ = topJ,
-    plotType = "sampleCorrelation",
-    normalized = FALSE,
-    annotation_col = "isExternal",
-    annotation_row = NA,
-    sampleCluster = NA,
-    minDeltaPsi = snakemake@config$aberrantSplicing$minDeltaPsi,
-    plotMeanPsi=FALSE,
-    plotCov = FALSE,
-    annotation_legend = TRUE,
-	annotation_colors = list(isExternal = c("FALSE" = anno_color_scheme[1],"TRUE" =  anno_color_scheme[2]))
-  )
-  before
-  after <- plotCountCorHeatmap(
-    object=fds,
-    type = type,
-    logit = TRUE,
-    topN = topN,
-    topJ = topJ,
-    plotType = "sampleCorrelation",
-    normalized = TRUE,
-    annotation_col = "isExternal",
-    annotation_row = NA,
-    sampleCluster = NA,
-    minDeltaPsi = snakemake@config$aberrantSplicing$minDeltaPsi,
-    plotMeanPsi=FALSE,
-    plotCov = FALSE,
-    annotation_legend = TRUE,
-	annotation_colors = list(isExternal = c("FALSE" = anno_color_scheme[1],"TRUE" =  anno_color_scheme[2]))
-  )
-  after
+  for(normalized in c(T,F)){
+    hm <- plotCountCorHeatmap(
+      object=fds,
+      type = type,
+      logit = TRUE,
+      topN = topN,
+      topJ = topJ,
+      plotType = "sampleCorrelation",
+      normalized = normalized,
+      annotation_col = "isExternal",
+      annotation_row = NA,
+      sampleCluster = NA,
+      minDeltaPsi = minDeltaPsi,
+      plotMeanPsi = FALSE,
+      plotCov = FALSE,
+      annotation_legend = TRUE,
+      annotation_colors = list(isExternal = c("FALSE" = anno_color_scheme[1], "TRUE" = anno_color_scheme[2]))
+    )
+    hm
+  }
 }
 
 #' ## Results
 res <- fread(snakemake@input$results)
+
+#' Total gene-level splicing outliers: `r nrow(res)`
+#' 
+
 file <- snakemake@output$res_html
 write_tsv(res, file = file)
 #+ echo=FALSE, results='asis'