clean up of summary scripts

gagneurlab · Apr 5, 2023 · ca9e588 · ca9e588
1 parent fa05486
commit ca9e588
Show file tree

Hide file tree

Showing 2 changed files with 59 additions and 40 deletions.
diff --git a/drop/modules/aberrant-splicing-pipeline/Counting/Summary.R b/drop/modules/aberrant-splicing-pipeline/Counting/Summary.R
@@ -44,54 +44,73 @@ devNull <- saveFraserDataSet(fdsMerge,dir=workingDir, name=paste0("raw-", datase
 #' Local: `r sum(!as.logical(fdsMerge@colData$isExternal))`  
 #' External: `r sum(as.logical(fdsMerge@colData$isExternal))`  
 #' 
+#' ## Number of introns:  
+#' Local (before filtering): `r length(rowRanges(fdsLocal, type = "j"))`  
+#' ```{asis, echo = has_external} 
+#' Merged with external counts (before filtering):
+#' ``` 
+#' ```{r, eval = has_external, echo=FALSE} 
+#' length(rowRanges(fdsMerge, type = "j"))
+#' ```
+#' After filtering: `r sum(mcols(fdsMerge, type="j")[,"passed"])`
+#' 
+#' ## Number of splice sites: 
+#' Local: `r length(rowRanges(fdsLocal, type = "theta"))`  
+#' ```{asis, echo = has_external}
+#' Merged with external counts:
+#' ``` 
+#' ```{r, eval = has_external, echo=FALSE} 
+#' length(rowRanges(fdsMerge, type = "theta"))
+#' ```
+#' 
+
+#' ```{asis, echo = has_external}
+#' ## Comparison of local and external counts  
 #' **Using external counts**  
 #' External counts introduce some complexity into the problem of counting junctions
 #' because it is unknown whether or not a junction is not counted (because there are no reads)
 #' compared to filtered and not present due to legal/personal sharing reasons. As a result,
 #' after merging the local (counted from BAM files) counts and the external counts, only the junctions that are 
 #' present in both remain. As a result it is likely that the number of junctions will decrease after merging.
 #' 
+#' ```
+#' ```{r, eval = has_external, echo=has_external}
+#' if(has_external){
+#'     externalCountIDs <- colData(fdsMerge)[as.logical(colData(fdsMerge)[,"isExternal"]),"sampleID"]
+#'     localCountIDs <- colData(fdsMerge)[!as.logical(colData(fdsMerge)[,"isExternal"]),"sampleID"]
 #' 
-#' ### Number of introns (psi5 or psi3) before and after merging:  
-#' Local: `r length(rowRanges(fdsLocal, type = "psi5"))`  
-#' Merged: `r length(rowRanges(fdsMerge, type = "psi5"))`  
+#'     cts <- K(fdsMerge,"psi5")
+#'     ctsLocal<- cts[,localCountIDs,drop=FALSE]
+#'     ctsExt<- cts[,externalCountIDs,drop=FALSE]
 #' 
-#' ### Number of splice sites (theta) before and after merging: 
-#' Local: `r length(rowRanges(fdsLocal, type = "theta"))`  
-#' Merged: `r length(rowRanges(fdsMerge, type = "theta"))`  
+#'     rowMeanLocal <- rowMeans(ctsLocal)
+#'     rowMeanExt <- rowMeans(ctsExt)
+#' 
+#'     dt <- data.table("Mean counts of local samples" = rowMeanLocal,
+#'                      "Mean counts of external samples" = rowMeanExt)
+#'                  
+#'     ggplot(dt,aes(x = `Mean counts of local samples`, y= `Mean counts of external samples`)) +
+#'        geom_hex() + theme_cowplot(font_size = 16) +
+#' 	   theme_bw() + scale_x_log10() + scale_y_log10() + 
+#'        geom_abline(slope = 1, intercept =0) +
+#'        scale_color_brewer(palette="Dark2") 
+#' }
+#' ```
 #' 
-
-#' ### Comparison of local and external counts  
-if(has_external){
-    externalCountIDs <- colData(fdsMerge)[as.logical(colData(fdsMerge)[,"isExternal"]),"sampleID"]
-    localCountIDs <- colData(fdsMerge)[!as.logical(colData(fdsMerge)[,"isExternal"]),"sampleID"]
-
-    cts <- K(fdsMerge,"psi5")
-    ctsLocal<- cts[,localCountIDs,drop=FALSE]
-    ctsExt<- cts[,externalCountIDs,drop=FALSE]
-
-    rowMeanLocal <- rowMeans(ctsLocal)
-    rowMeanExt <- rowMeans(ctsExt)
-
-    dt <- data.table("Mean counts of local samples" = rowMeanLocal,
-                     "Mean counts of external samples" = rowMeanExt)
-
-    ggplot(dt,aes(x = `Mean counts of local samples`, y= `Mean counts of external samples`)) +
-       geom_hex() + theme_cowplot(font_size = 16) +
-	   theme_bw() + scale_x_log10() + scale_y_log10() + 
-       geom_abline(slope = 1, intercept =0) +
-       scale_color_brewer(palette="Dark2") 
-}else{
-	print("No external counts, comparison is ommitted")
-}
 
 #' ## Expression filtering
-#' Min expression cutoff: `r snakemake@config$aberrantSplicing$minExpressionInOneSample`
-plotFilterExpression(fdsMerge) + theme_cowplot(font_size = 16)
+#' The expression filtering step removes introns that are lowly expressed. The requirements for an intron to pass this filter are:
+#' 
+#' * at least 1 sample needs to have a count (K) of `r snakemake@config$aberrantSplicing$minExpressionInOneSample` or more reads for the intron
+#' * at least `r 100*(1-snakemake@config$aberrantSplicing$quantileForFiltering)`% of the samples need to have a total of at least `r snakemake@config$aberrantSplicing$quantileMinExpression` reads for the splice metric denominator (N) of the intron
+plotFilterExpression(fdsMerge) + 
+    labs(title="", x="Mean Intron Expression", y="Introns") +
+    theme_cowplot(font_size = 16)
 
 #' ## Variability filtering
-#' Variability cutoff: `r snakemake@config$aberrantSplicing$minDeltaPsi`
-plotFilterVariability(fdsMerge) + theme_cowplot(font_size = 16)
-
-#' Introns that passed filter (after merging)
-table(mcols(fdsMerge, type="j")[,"passed"])
+#' The variability filtering step removes introns that have no or little variability in the splice metric values across samples. The requirement for an intron to pass this filter is:
+#' 
+#' * at least 1 sample has a difference of at least `r snakemake@config$aberrantSplicing$minDeltaPsi` in the splice metric compared to the mean splice metric of the intron
+plotFilterVariability(fdsMerge) + 
+    labs(title="", y="Introns") +
+    theme_cowplot(font_size = 16)
diff --git a/drop/modules/aberrant-splicing-pipeline/FRASER/Summary.R b/drop/modules/aberrant-splicing-pipeline/FRASER/Summary.R
@@ -39,9 +39,9 @@ hasExternal <- length(levels(colData(fds)$isExternal) > 1)
 
 #' Number of samples: `r nrow(colData(fds))`
 #' 
-#' Number of introns: `r length(rowRanges(fds, type = "psi5"))`
+#' Number of introns: `r length(rowRanges(fds, type = "j"))`
 #' 
-#' Number of splice sites: `r length(rowRanges(fds, type = "theta"))`
+#' Number of splice sites: `r length(rowRanges(fds, type = "ss"))`
 
 # used for most plots
 dataset_title <- paste0("Dataset: ", dataset, "--", annotation)
@@ -70,7 +70,7 @@ topJ <- 10000
 anno_color_scheme <- brewer.pal(n = 3, name = 'Dark2')[1:2]
 
 for(type in psiTypes){
-  for(normalized in c(T,F)){
+  for(normalized in c(F,T)){
     hm <- plotCountCorHeatmap(
       object=fds,
       type = type,