removed normalization for dipLogR to 0; fixed plotting issue with y-a…

…xis for icn; fixed issue with calculating discordant tcn
mskcc · Apr 27, 2020 · f2b3a46 · f2b3a46
1 parent f8d4faa
commit f2b3a46
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 17 deletions.
diff --git a/R/check-fit.R b/R/check-fit.R
@@ -175,10 +175,10 @@ check_fit = function(facets_output,
         discordant_lcn = (lcn.em.original != lcn.original | (is.na(lcn.em.original) | is.na(lcn.original))) & !(is.na(lcn.em.original) & is.na(lcn.original))
     )][(discordant_tcn == TRUE | discordant_lcn == TRUE) & tcn < 10]
     discordant_stats = discordant_segs[, list(
-        n_discordant_tcn = sum(discordant_tcn & !discordant_lcn),
-        length_discordant_tcn = sum(length[discordant_tcn & !discordant_lcn]),
-        n_discordant_lcn = sum(length[!discordant_tcn & discordant_lcn]),
-        length_discordant_lcn = sum(length[!discordant_tcn & discordant_lcn]),
+        n_discordant_tcn = sum(discordant_tcn),
+        length_discordant_tcn = sum(length[discordant_tcn ]),
+        n_discordant_lcn = sum(length[discordant_lcn]),
+        length_discordant_lcn = sum(length[discordant_lcn]),
         n_discordant_both = sum(discordant_tcn & discordant_lcn),
         length_discordant_both = sum(length[discordant_tcn & discordant_lcn])
     )]

diff --git a/R/plot-facets.R b/R/plot-facets.R
@@ -13,7 +13,6 @@
 #' @param plotX If \code{TRUE}, includes chromosome X.
 #' @param genome Genome build.
 #' @param highlight_gene Highlight gene(s), provide gene symbol or mapped positions (internally).
-#' @param adjust_dipLogR Normalize by sample dipLogR.
 #' @param method When available, choose between plotting solution from \code{em} or \code{cncf} algorithm.
 #' @param subset_snps Subset the SNP profile to reduce weight of plotting, supply a factor by which to reduce or \code{TRUE} for default.
 #' @param plot_chroms Chromosomes to plot when using \code{closeup_plot}.
@@ -37,7 +36,6 @@ cnlr_plot = function(facets_data,
                      plotX = FALSE,
                      genome = c('hg19', 'hg18', 'hg38'),
                      highlight_gene = NULL,
-                     adjust_dipLogR = TRUE,
                      subset_snps = NULL,
                      return_object = FALSE) {
 
@@ -67,14 +65,7 @@ cnlr_plot = function(facets_data,
     ymax = ceiling(max(segs$cnlr.median, na.rm = T))
     if (ymin > -3) ymin = -3
     if (ymax < 3) ymax = 3
-
-    if (adjust_dipLogR) {
-        snps$cnlr = snps$cnlr - dipLogR
-        my_starts$cnlr_median = my_starts$cnlr_median - dipLogR
-        my_ends$cnlr_median = my_ends$cnlr_median - dipLogR
-        dipLogR = Inf
-    }
-
+
     if (!is.null(subset_snps)) {
         if (subset_snps == TRUE) {
             snps = subset_snps(snps)
@@ -90,7 +81,7 @@ cnlr_plot = function(facets_data,
         geom_point(aes(y = cnlr, x = chr_maploc), pch = 19, col = pt_cols, size = .4) +
         scale_x_continuous(breaks = mid, labels = names(mid), expand = c(.01, 0)) +
         scale_y_continuous(breaks = scales::pretty_breaks(), limits = c(ymin, ymax)) +
-        geom_hline(yintercept = 0, color = 'sandybrown', size = .8) +
+        geom_hline(yintercept = dipLogR, color = 'sandybrown', size = .8) +
         geom_segment(data = segs, aes(x = my_starts$chr_maploc, xend = my_ends$chr_maploc,
                                       y = my_starts$cnlr_median, yend = my_ends$cnlr_median),
                      col = 'red3', size = 1, lineend = 'butt') +
@@ -297,11 +288,15 @@ icn_plot = function(facets_data,
         tcnscaled[segs$tcn.em > 5 & !is.na(segs$tcn.em)] = (5 + (tcnscaled[segs$tcn.em > 5 & !is.na(segs$tcn.em)] - 5) / 3)
         lcn = rep(segs$lcn.em, segs$num.mark)
         my_ylab = 'Integer copy number (EM)'        
+        axis_breaks = c(0:5, 5 + (sort(unique(segs[tcn.em >5, tcn.em])) - 5)/3)
+        axis_labels = c(0:5, sort(unique(segs[tcn.em >5, tcn.em])))
     } else if (method == 'cncf') {
         tcnscaled = segs$tcn
         tcnscaled[segs$tcn > 5 & !is.na(segs$tcn)] = (5 + (tcnscaled[segs$tcn > 5 & !is.na(segs$tcn)] - 5) / 3)
         lcn = rep(segs$lcn, segs$num.mark)
         my_ylab = 'Integer copy number (CNCF)'        
+        axis_breaks = c(0:5, 5 + (sort(unique(segs[tcn >5, tcn])) - 5)/3)
+        axis_labels = c(0:5, sort(unique(segs[tcn >5, tcn])))
     }
     tcn = rep(tcnscaled, segs$num.mark)
 
@@ -314,8 +309,6 @@ icn_plot = function(facets_data,
     my_lcn_starts = snps[starts, c('chr_maploc', 'lcn')]
     my_lcn_ends = snps[ends, c('chr_maploc', 'lcn')]
 
-    axis_breaks = c(0:5, 5 + (sort(unique(segs[tcn >5, tcn])) - 5)/3)
-    axis_labels = c(0:5, sort(unique(segs[tcn >5, tcn])))
 
     icn = ggplot(segs) +
         geom_segment(col = 'red', size = 1,