R Under development (unstable) (2024-02-16 r85931 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library("PSCBS") PSCBS v0.67.0 successfully loaded. See ?PSCBS for help. > > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > # Simulating copy-number data > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > set.seed(0xBEEF) > > # Number of loci > J <- 1000 > > x <- sort(runif(J, max=J)) * 1e5 > > mu <- double(J) > mu[200:300] <- mu[200:300] + 1 > mu[350:400] <- NA # centromere > mu[650:800] <- mu[650:800] - 1 > eps <- rnorm(J, sd=1/2) > y <- mu + eps > > outliers <- seq(from=1L, to=J, length.out=0.2*J) > y[outliers] <- y[outliers] + 1.5 > > w <- rep(1.0, times=J) > w[outliers] <- 0.01 > > data <- data.frame(chromosome=1L, x=x, y=y) > dataW <- cbind(data, w=w) > > > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > # Single-chromosome segmentation > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > par(mar=c(2,3,0.2,1)+0.1) > # Segment without weights > fit <- segmentByCBS(data) > sampleName(fit) <- "CBS_Example" > print(fit) sampleName chromosome start end nbrOfLoci mean 1 CBS_Example 1 136857.7 19138391 199 0.2712 2 CBS_Example 1 19138391.4 28682180 101 1.2168 3 CBS_Example 1 28682180.1 64690253 298 0.3027 4 CBS_Example 1 64690253.3 80738828 151 -0.7101 5 CBS_Example 1 80738828.3 99932904 200 0.3655 > plotTracks(fit) Warning message: In plotTracks.CBS(fit) : Setting default 'Clim' assuming the signal type is 'ratio' because signalType(fit) is unknown ('NA'). Use signalType(fit) <- 'ratio' to avoid this warning. > ## Highlight outliers (they pull up the mean levels) > points(x[outliers]/1e6, y[outliers], col="purple") > > # Segment without weights but with median > fitM <- segmentByCBS(data, avg="median") > sampleName(fitM) <- "CBS_Example (median)" > print(fitM) sampleName chromosome start end nbrOfLoci mean 1 CBS_Example (median) 1 136857.7 19138391 199 0.1203255 2 CBS_Example (median) 1 19138391.4 28682180 101 0.9949202 3 CBS_Example (median) 1 28682180.1 64690253 298 0.1471793 4 CBS_Example (median) 1 64690253.3 80738828 151 -0.8770443 5 CBS_Example (median) 1 80738828.3 99932904 200 0.2211061 > drawLevels(fitM, col="magenta", lty=3) NULL > > # Segment with weights > fitW <- segmentByCBS(dataW, avg="median") > sampleName(fitW) <- "CBS_Example (weighted)" > print(fitW) sampleName chromosome start end nbrOfLoci mean 1 CBS_Example (weighted) 1 136857.7 19138391 199 -0.02220950 2 CBS_Example (weighted) 1 19138391.4 28682180 101 0.92421628 3 CBS_Example (weighted) 1 28682180.1 64690253 298 -0.02364830 4 CBS_Example (weighted) 1 64690253.3 80738828 151 -1.04750872 5 CBS_Example (weighted) 1 80738828.3 99932904 200 0.08961195 > drawLevels(fitW, col="red") NULL > > # Segment with weights and median > fitWM <- segmentByCBS(dataW, avg="median") > sampleName(fitWM) <- "CBS_Example (weighted median)" > print(fitWM) sampleName chromosome start end nbrOfLoci 1 CBS_Example (weighted median) 1 136857.7 19138391 199 2 CBS_Example (weighted median) 1 19138391.4 28682180 101 3 CBS_Example (weighted median) 1 28682180.1 64690253 298 4 CBS_Example (weighted median) 1 64690253.3 80738828 151 5 CBS_Example (weighted median) 1 80738828.3 99932904 200 mean 1 -0.02220950 2 0.92421628 3 -0.02364830 4 -1.04750872 5 0.08961195 > drawLevels(fitWM, col="orange", lty=3) NULL > > legend("topright", bg="white", legend=c("outliers", "non-weighted CBS (mean)", "non-weighted CBS (median)", "weighted CBS (mean)", "weighted CBS (median)"), col=c("purple", "purple", "magenta", "red", "orange"), lwd=c(NA,3,3,3,3), lty=c(NA,1,3,1,3), pch=c(1,NA,NA,NA,NA)) > > ## Assert that weighted segment means are less biased > dmean <- getSegments(fit)$mean - getSegments(fitW)$mean > cat("Segment mean differences:\n") Segment mean differences: > print(dmean) [1] 0.2934095 0.2925837 0.3263483 0.3374087 0.2758881 > stopifnot(all(dmean > 0, na.rm=TRUE)) > > dmean <- getSegments(fitM)$mean - getSegments(fitWM)$mean > cat("Segment median differences:\n") Segment median differences: > print(dmean) [1] 0.14253502 0.07070392 0.17082758 0.17046439 0.13149418 > stopifnot(all(dmean > 0, na.rm=TRUE)) > > > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > # Multi-chromosome segmentation > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > data2 <- data > data2$chromosome <- 2L > data <- rbind(data, data2) > dataW <- cbind(data, w=w) > > par(mar=c(2,3,0.2,1)+0.1) > # Segment without weights > fit <- segmentByCBS(data) > sampleName(fit) <- "CBS_Example" > print(fit) sampleName chromosome start end nbrOfLoci mean 1 CBS_Example 1 136857.7 19138391 199 0.2712 2 CBS_Example 1 19138391.4 28682180 101 1.2168 3 CBS_Example 1 28682180.1 64690253 298 0.3027 4 CBS_Example 1 64690253.3 80738828 151 -0.7101 5 CBS_Example 1 80738828.3 99932904 200 0.3655 6 NA NA NA NA NA 7 CBS_Example 2 136857.7 19138391 199 0.2712 8 CBS_Example 2 19138391.4 28682180 101 1.2168 9 CBS_Example 2 28682180.1 64690253 298 0.3027 10 CBS_Example 2 64690253.3 80738828 151 -0.7101 11 CBS_Example 2 80738828.3 99932904 200 0.3655 > plotTracks(fit, Clim=c(-3,3)) > > # Segment without weights but with median > fitM <- segmentByCBS(data, avg="median") > sampleName(fitM) <- "CBS_Example (median)" > print(fitM) sampleName chromosome start end nbrOfLoci mean 1 CBS_Example (median) 1 136857.7 19138391 199 0.1203255 2 CBS_Example (median) 1 19138391.4 28682180 101 0.9949202 3 CBS_Example (median) 1 28682180.1 64690253 298 0.1471793 4 CBS_Example (median) 1 64690253.3 80738828 151 -0.8770443 5 CBS_Example (median) 1 80738828.3 99932904 200 0.2211061 6 NA NA NA NA NA 7 CBS_Example (median) 2 136857.7 19138391 199 0.1203255 8 CBS_Example (median) 2 19138391.4 28682180 101 0.9949202 9 CBS_Example (median) 2 28682180.1 64690253 298 0.1471793 10 CBS_Example (median) 2 64690253.3 80738828 151 -0.8770443 11 CBS_Example (median) 2 80738828.3 99932904 200 0.2211061 > drawLevels(fitM, col="magenta", lty=3) NULL > > # Segment with weights > fitW <- segmentByCBS(dataW, avg="median") > sampleName(fitW) <- "CBS_Example (weighted)" > print(fitW) sampleName chromosome start end nbrOfLoci mean 1 CBS_Example (weighted) 1 136857.7 19138391 199 -0.02220950 2 CBS_Example (weighted) 1 19138391.4 28682180 101 0.92421628 3 CBS_Example (weighted) 1 28682180.1 64690253 298 -0.02364830 4 CBS_Example (weighted) 1 64690253.3 80738828 151 -1.04750872 5 CBS_Example (weighted) 1 80738828.3 99932904 200 0.08961195 6 NA NA NA NA NA 7 CBS_Example (weighted) 2 136857.7 19138391 199 -0.02220950 8 CBS_Example (weighted) 2 19138391.4 28682180 101 0.92421628 9 CBS_Example (weighted) 2 28682180.1 64690253 298 -0.02364830 10 CBS_Example (weighted) 2 64690253.3 80738828 151 -1.04750872 11 CBS_Example (weighted) 2 80738828.3 99932904 200 0.08961195 > drawLevels(fitW, col="red") NULL > > # Segment with weights and median > fitWM <- segmentByCBS(dataW, avg="median") > sampleName(fitWM) <- "CBS_Example (weighted median)" > print(fitWM) sampleName chromosome start end nbrOfLoci 1 CBS_Example (weighted median) 1 136857.7 19138391 199 2 CBS_Example (weighted median) 1 19138391.4 28682180 101 3 CBS_Example (weighted median) 1 28682180.1 64690253 298 4 CBS_Example (weighted median) 1 64690253.3 80738828 151 5 CBS_Example (weighted median) 1 80738828.3 99932904 200 6 NA NA NA NA 7 CBS_Example (weighted median) 2 136857.7 19138391 199 8 CBS_Example (weighted median) 2 19138391.4 28682180 101 9 CBS_Example (weighted median) 2 28682180.1 64690253 298 10 CBS_Example (weighted median) 2 64690253.3 80738828 151 11 CBS_Example (weighted median) 2 80738828.3 99932904 200 mean 1 -0.02220950 2 0.92421628 3 -0.02364830 4 -1.04750872 5 0.08961195 6 NA 7 -0.02220950 8 0.92421628 9 -0.02364830 10 -1.04750872 11 0.08961195 > drawLevels(fitWM, col="orange", lty=3) NULL > > legend("topright", bg="white", legend=c("outliers", "non-weighted CBS (mean)", "non-weighted CBS (median)", "weighted CBS (mean)", "weighted CBS (median)"), col=c("purple", "purple", "magenta", "red", "orange"), lwd=c(NA,3,3,3,3), lty=c(NA,1,3,1,3), pch=c(1,NA,NA,NA,NA)) > > ## Assert that weighted segment means are less biased > dmean <- getSegments(fit)$mean - getSegments(fitW)$mean > cat("Segment mean differences:\n") Segment mean differences: > print(dmean) [1] 0.2934095 0.2925837 0.3263483 0.3374087 0.2758881 NA 0.2934095 [8] 0.2925837 0.3263483 0.3374087 0.2758881 > stopifnot(all(dmean > 0, na.rm=TRUE)) > > dmean <- getSegments(fitM)$mean - getSegments(fitWM)$mean > cat("Segment median differences:\n") Segment median differences: > print(dmean) [1] 0.14253502 0.07070392 0.17082758 0.17046439 0.13149418 NA [7] 0.14253502 0.07070392 0.17082758 0.17046439 0.13149418 > stopifnot(all(dmean > 0, na.rm=TRUE)) > > proc.time() user system elapsed 1.43 0.09 1.75