Browse code

Keep allosome log ratio.

Markus Riester authored on 28/03/2024 18:50:25
Showing 2 changed files

... ...
@@ -197,6 +197,25 @@ normalDB.min.coverage, normalDB.max.missing) {
197 197
     intervalsUsed
198 198
 }
199 199
 
200
+.filterIntervalsAllosome <- function(intervalsUsed, tumor, sex) {
201
+    if (is.null(sex)) return(intervalsUsed)
202
+    sex.chr <- .getSexChr(seqlevels(tumor))
203
+    remove.chrs <- c()
204
+    if (sex == "M" || sex == "?") {
205
+        remove.chrs <- sex.chr
206
+    } else if (sex == "F") {
207
+        remove.chrs <- sex.chr[2]
208
+    }
209
+    nBefore <- sum(intervalsUsed)
210
+    intervalsUsed <- intervalsUsed & !seqnames(tumor) %in% remove.chrs
211
+    nAfter <- sum(intervalsUsed)
212
+    if (nAfter < nBefore) {
213
+        flog.info("Removing %i non-diploid allosome intervals.", nBefore - nAfter)
214
+    }
215
+    intervalsUsed
216
+}
217
+
218
+
200 219
 .filterIntervalsTotalNormalCoverage <- function(intervalsUsed, normal,
201 220
     min.targeted.base, min.coverage) {
202 221
     if (is.null(min.targeted.base) || is.null(min.coverage)) {
... ...
@@ -416,11 +416,10 @@ runAbsoluteCN <- function(normal.coverage.file = NULL,
416 416
     }
417 417
 
418 418
     sex <- .getSex(match.arg(sex), normal, tumor)
419
-    tumor <- .fixAllosomeCoverage(sex, tumor)
420
-
421 419
     if (!is.null(interval.file)) {
422 420
         tumor <- .addGCData(tumor, interval.file)
423 421
     }
422
+
424 423
     if (is.null(centromeres) && !missing(genome)) {
425 424
         centromeres <- .getCentromerePositions(centromeres, genome,
426 425
             if (is.null(tumor)) NULL else .getSeqlevelsStyle(tumor))
... ...
@@ -442,6 +441,16 @@ runAbsoluteCN <- function(normal.coverage.file = NULL,
442 441
     # chr.hash is an internal data structure, so we need to do this separately.
443 442
     intervalsUsed <- .filterIntervalsChrHash(intervalsUsed, tumor, chr.hash)
444 443
     intervalsUsed <- .filterIntervalsCentromeres(intervalsUsed, tumor, centromeres)
444
+    intervalsUsedAllosome <- intervalsUsed
445
+    intervalsUsed <- .filterIntervalsAllosome(intervalsUsed, tumor, sex)
446
+
447
+    if (!is.null(normalDB$sd$weights)) {
448
+        tumor$weights <- subsetByOverlaps(normalDB$sd$weights, tumor)$weights
449
+    }
450
+    tumorAllosome <- tumor
451
+    tumorAllosome$log.ratio <- log.ratio
452
+    tumorAllosome <- tumorAllosome[!intervalsUsed & intervalsUsedAllosome, ]
453
+    
445 454
     intervalsUsed <- which(intervalsUsed)
446 455
     if (length(tumor) != length(normal) ||
447 456
         length(tumor) != length(log.ratio)) {
... ...
@@ -468,9 +477,6 @@ runAbsoluteCN <- function(normal.coverage.file = NULL,
468 477
          flog.info("Mean off-target bin size: %.0f",
469 478
             mean(width(tumor[!tumor$on.target]), na.rm = TRUE))
470 479
     }
471
-    if (!is.null(normalDB$sd$weights)) {
472
-        tumor$weights <- subsetByOverlaps(normalDB$sd$weights, tumor)$weights
473
-    }
474 480
     # not needed anymore
475 481
     normalDB <- NULL
476 482
 
... ...
@@ -1145,7 +1151,7 @@ runAbsoluteCN <- function(normal.coverage.file = NULL,
1145 1151
     list(
1146 1152
         candidates = candidate.solutions,
1147 1153
         results = results,
1148
-        input = list(tumor = tumor.coverage.file, normal = normal.coverage.file,
1154
+        input = list(tumor = tumor.coverage.file, normal = normal.coverage.file, allosome = tumorAllosome,
1149 1155
             log.ratio = GRanges(normal[, 1], on.target = normal$on.target, log.ratio = log.ratio),
1150 1156
             log.ratio.sdev = sd.seg, mapd = mapd, vcf = vcf, sampleid = sampleid,
1151 1157
             test.num.copy = test.num.copy,