Mercurial > repos > artbio > cnv_facets

comparison plot_facets_enhanced-v22.R @ 3:d1914f4d9daf draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/main/tools/facets commit 64ac36125f04497dd51028f307e059fca9ec0503
author artbio
date Sun, 05 Oct 2025 00:55:34 +0000
parents
children
comparison
equal deleted inserted replaced
2:66a56502199d 3:d1914f4d9daf
1 plot_facets_enhanced <- function(x, emfit = NULL, clustered = FALSE,
2 plot.type = c("em", "naive", "both", "none"),
3 sname = NULL, add.legend = TRUE) {
4 def.par <- par(no.readonly = TRUE)
5 plot.type <- match.arg(plot.type)
6
7 # Setup layout
8 if (plot.type == "none") {
9 layout(matrix(1:2, ncol = 1))
10 }
11 if (plot.type == "em") {
12 layout(matrix(rep(1:4, c(9, 9, 6, 1)), ncol = 1))
13 }
14 if (plot.type == "naive") {
15 layout(matrix(rep(1:4, c(9, 9, 6, 1)), ncol = 1))
16 }
17 if (plot.type == "both") {
18 layout(matrix(rep(1:6, c(9, 9, 6, 1, 6, 1)), ncol = 1))
19 }
20
21 par(
22 mar = c(0.25, 3, 0.25, 1), mgp = c(1.75, 0.6, 0),
23 oma = c(3, 0, 1.25, 0)
24 )
25
26 # Prepare data
27 jseg <- x$jointseg
28 chrbdry <- which(diff(jseg$chrom) != 0)
29
30 if (missing(emfit)) {
31 out <- x$out
32 if (plot.type == "em" | plot.type == "both") {
33 warning("emfit is missing; plot.type set to naive")
34 plot.type <- "naive"
35 }
36 } else {
37 out <- emfit$cncf
38 out$tcn <- x$out$tcn
39 out$lcn <- x$out$lcn
40 out$cf <- x$out$cf
41 }
42
43 if (clustered) {
44 cnlr.median <- out$cnlr.median.clust
45 mafR <- out$mafR.clust
46 mafR[is.na(mafR)] <- out$mafR[is.na(mafR)]
47 } else {
48 cnlr.median <- out$cnlr.median
49 mafR <- out$mafR
50 }
51
52 mafR <- abs(mafR)
53 chrcol <- 1 + rep(out$chrom - 2 * floor(out$chrom / 2), out$num.mark)
54 nn <- cumsum(table(jseg$chrom[is.finite(jseg$cnlr)]))
55 segbdry <- cumsum(c(0, out$num.mark))
56 segstart <- segbdry[-length(segbdry)]
57 segend <- segbdry[-1]
58
59 # Plot 1: log-ratio
60 plot(jseg$cnlr[is.finite(jseg$cnlr)],
61 pch = ".", cex = 2,
62 col = c("grey", "lightblue", "azure4", "slateblue")[chrcol],
63 ylab = "log-ratio", xaxt = "n"
64 )
65 abline(v = chrbdry, lwd = 0.25)
66
67 # Calculate reference lines
68 global_median <- median(jseg$cnlr, na.rm = TRUE)
69 diploid_logr <- x$dipLogR
70 lines_diff <- abs(global_median - diploid_logr)
71
72 # Only show reference lines if they are sufficiently different (> 0.1)
73 legend_items <- c("Segments")
74 legend_cols <- c(2)
75 legend_lwd <- c(1.75)
76 legend_lty <- c(1)
77
78 if (lines_diff > 0.1) {
79 # Lines are different enough to show both
80 abline(h = global_median, col = "green2", lwd = 2.5)
81 abline(h = diploid_logr, col = "magenta4", lwd = 2.5, lty = 2)
82 legend_items <- c(
83 legend_items,
84 sprintf("Global median (%.3f)", global_median),
85 sprintf("Diploid LogR (%.3f)", diploid_logr)
86 )
87 legend_cols <- c(legend_cols, "green2", "magenta4")
88 legend_lwd <- c(legend_lwd, 2.5, 2.5)
89 legend_lty <- c(legend_lty, 1, 2)
90 } else {
91 # Lines are too close, show only one with combined label
92 abline(h = global_median, col = "green2", lwd = 2.5)
93 legend_items <- c(
94 legend_items,
95 sprintf("Median/Diploid (%.3f)", global_median)
96 )
97 legend_cols <- c(legend_cols, "green2")
98 legend_lwd <- c(legend_lwd, 2.5)
99 legend_lty <- c(legend_lty, 1)
100 }
101
102 segments(segstart, cnlr.median, segend, cnlr.median, lwd = 1.75, col = 2)
103
104 # Add legend for log-ratio plot
105 if (add.legend) {
106 legend("topright",
107 legend = legend_items,
108 col = legend_cols,
109 lwd = legend_lwd,
110 lty = legend_lty,
111 bty = "n",
112 cex = 0.7
113 )
114 }
115
116 # Plot 2: log-odds-ratio
117 plot(jseg$valor[is.finite(jseg$cnlr)],
118 pch = ".", cex = 2.5,
119 col = c("grey", "lightblue", "azure4", "slateblue")[chrcol],
120 ylab = "log-odds-ratio", ylim = c(-4, 4), xaxt = "n"
121 )
122 abline(v = chrbdry, lwd = 0.25)
123 segments(segstart, sqrt(mafR), segend, sqrt(mafR), lwd = 1.75, col = 2)
124 segments(segstart, -sqrt(mafR), segend, -sqrt(mafR), lwd = 1.75, col = 2)
125
126 # Add legend for log-odds-ratio plot
127 if (add.legend) {
128 legend("topright",
129 legend = c("BAF segments"),
130 col = c(2),
131 lwd = c(1.75),
132 bty = "n",
133 cex = 0.7
134 )
135 }
136
137 cfpalette <- c(colorRampPalette(c("white", "steelblue"))(10), "bisque2")
138
139 # Plot 3: copy number (naive)
140 if (plot.type == "naive" | plot.type == "both") {
141 out$tcn[out$tcn > 10] <- 9 + log10(out$tcn[out$tcn > 10])
142 ii <- which(out$lcn > 5)
143 if (length(ii) > 0) {
144 out$lcn[ii] <- 5 + log10(out$lcn[ii])
145 }
146
147 plot(c(0, length(jseg$cnlr)), c(0, max(out$tcn)),
148 type = "n",
149 ylab = "copy number (nv)", xaxt = "n"
150 )
151 abline(v = chrbdry, lwd = 0.25)
152 segments(segstart, out$lcn, segend, out$lcn, lwd = 1.75, col = 2)
153 segments(segstart, out$tcn, segend, out$tcn, lwd = 1.75, col = 1)
154
155 # Add legends for copy number plot (including CF legend)
156 if (add.legend) {
157 # CN legend on top left
158 legend("topleft",
159 legend = c("Total CN (TCN)", "Minor CN (LCN)"),
160 col = c(1, 2),
161 lwd = 1.75,
162 bty = "n",
163 cex = 0.7
164 )
165
166 # CF legend on top right
167 legend("topright",
168 legend = c(
169 "CF = 0 (normal)", "CF = 0.2 (low tumor)",
170 "CF = 0.8 (high tumor)", "CF = 1.0 (pure tumor)"
171 ),
172 fill = c("white", "lightblue", "steelblue", "tan"),
173 border = "black",
174 bty = "o",
175 box.col = "black",
176 box.lwd = 1.5,
177 cex = 0.65,
178 title = "Cellular Fraction (CF)",
179 bg = "white"
180 )
181 }
182
183 # CF bar (naive)
184 plot(c(0, length(jseg$cnlr)), 0:1,
185 type = "n", ylab = "",
186 xaxt = "n", yaxt = "n"
187 )
188 mtext("cf-nv", side = 2, at = 0.5, line = 0.3, las = 2, cex = 0.75)
189 cfcol <- cfpalette[round(10 * out$cf + 0.501)]
190 rect(segstart, 0, segend, 1, col = cfcol, border = NA)
191 }
192
193 # Plot 4: copy number (EM)
194 if (plot.type == "em" | plot.type == "both") {
195 out$tcn.em[out$tcn.em > 10] <- 9 + log10(out$tcn.em[out$tcn.em > 10])
196 ii <- which(out$lcn.em > 5)
197 if (length(ii) > 0) {
198 out$lcn.em[ii] <- 5 + log10(out$lcn.em[ii])
199 }
200
201 plot(c(0, length(jseg$cnlr)), c(0, max(out$tcn.em)),
202 type = "n", ylab = "copy number (em)", xaxt = "n"
203 )
204 abline(v = chrbdry, lwd = 0.25)
205 segments(segstart, out$lcn.em, segend, out$lcn.em, lwd = 1.75, col = 2)
206 segments(segstart, out$tcn.em, segend, out$tcn.em, lwd = 1.75, col = 1)
207
208 # Add legends for EM copy number plot (including CF legend)
209 if (add.legend) {
210 # CN legend on top left
211 legend("topleft",
212 legend = c("Total CN (TCN)", "Minor CN (LCN)"),
213 col = c(1, 2),
214 lwd = 1.75,
215 bty = "n",
216 cex = 0.7
217 )
218
219 # CF legend on top right
220 legend("topright",
221 legend = c(
222 "CF = 0 (normal)", "CF = 0.2 (low tumor)",
223 "CF = 0.8 (high tumor)", "CF = 1.0 (pure tumor)"
224 ),
225 fill = c("white", "lightblue", "steelblue", "tan"),
226 border = "black",
227 bty = "o",
228 box.col = "black",
229 box.lwd = 1.5,
230 cex = 0.65,
231 title = "Cellular Fraction (CF)",
232 bg = "white"
233 )
234 }
235
236 # CF bar (EM)
237 plot(c(0, length(jseg$cnlr)), 0:1,
238 type = "n", ylab = "",
239 xaxt = "n", yaxt = "n"
240 )
241 mtext("cf-em", side = 2, at = 0.5, line = 0.2, las = 2, cex = 0.75)
242 cfcol <- cfpalette[round(10 * out$cf.em + 0.501)]
243 rect(segstart, 0, segend, 1, col = cfcol, border = NA)
244 }
245
246 # X-axis with chromosome labels
247 chromlevels <- x$chromlevels
248 if (is.null(chromlevels)) {
249 chromlevels <- 1:length(nn)
250 }
251 axis(
252 labels = chromlevels, side = 1,
253 at = (nn + c(0, nn[-length(nn)])) / 2, cex = 0.65
254 )
255 mtext(side = 1, line = 1.75, "Chromosome", cex = 0.8)
256
257 # Sample name title
258 if (!missing(sname)) {
259 mtext(sname, side = 3, line = 0, outer = TRUE, cex = 0.8)
260 }
261
262 par(def.par)
263 }