Merge pull request #20 from markowetzlab/dev

v1.2.0 development PR
markowetzlab · Dec 8, 2023 · 9cd18c7 · 9cd18c7
2 parents 46b9278 + d3af314
commit 9cd18c7
Show file tree

Hide file tree

Showing 66 changed files with 1,046 additions and 196 deletions.
diff --git a/.gitignore b/.gitignore
@@ -6,3 +6,4 @@ inst/4_*
 /doc/
 /Meta/
 inst/doc
+tests/testthat/Rplots.pdf
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: CINSignatureQuantification
-Title: Simple and quick measuring of copy number signatures in cancers
-Version: 1.1.2
+Title: Computing and measuring of copy number signatures in cancers
+Version: 1.2.0
 Authors@R: 
     c(person(given = "Philip S",
            family = "Smith",
@@ -27,13 +27,11 @@ biocViews:
 Depends:
     R (>= 3.6.0)
 Imports:
-    base,
     Biobase (>= 2.46.0),
     data.table (>= 1.14),
     graphics (>= 1.3.7),
     limSolve (>= 1.5.6),
     methods,
-    parallel,
     stats,
     stringr (>= 1.4),
     utils
@@ -47,6 +45,6 @@ License: ASL + file LICENSE
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.2
+RoxygenNote: 7.2.3
 VignetteBuilder: knitr
 Config/testthat/edition: 3
diff --git a/R/AllClasses.R b/R/AllClasses.R
@@ -36,7 +36,7 @@ ExpQuant <- setClass("ExpQuant",
                                   feature.method = "character"
                      ),
                      prototype = list(experimentName = NULL,
-                                      init.date = as.character(Sys.time()),
+                                      init.date = NULL,
                                       last.modified = "NA",
                                       samples.full = NULL,
                                       samples.current = NULL,

diff --git a/R/AllGenerics.R b/R/AllGenerics.R
@@ -154,6 +154,10 @@ setGeneric("getFeatures", function(object,feat=NULL) standardGeneric("getFeature
 #' @param cores Number of CPU threads/cores to utilise via doParallel. Default
 #'   is 1. Maximum number is equal to the number of features to extract (drews &
 #'   mac methods = 6 features).
+#' @param DCIN Threshold for required number of non-diploid segments to compute
+#'   copy number features (and subsequently copy number signatures) using method
+#'   "drews". Default is 20. This parameter should not need to be changed and
+#'   will affect feature values and signature activity.
 #' @return A CNQuant class object with extracted features stored in the
 #'   "featData" slot
 #' @examples
@@ -165,7 +169,7 @@ setGeneric("getFeatures", function(object,feat=NULL) standardGeneric("getFeature
 #' @docType methods
 #' @rdname calculateFeatures-methods
 #'
-setGeneric("calculateFeatures",function(object, method="drews",smooth.diploid=TRUE,cores=1)
+setGeneric("calculateFeatures",function(object, method="drews",smooth.diploid=TRUE,cores=1,DCIN = 20)
     standardGeneric("calculateFeatures"))
 
 #' calculateSampleByComponentMatrix

diff --git a/R/addSampleFeatures.R b/R/addSampleFeatures.R
@@ -16,7 +16,10 @@
 #'   cnobj <- addsampleFeatures(cnobj,sample.data=test.sample.features)
 #'   getSamplefeatures(cnobj)
 #' @export
-addsampleFeatures <- function(object,sample.data=NULL,id.col = "sample"){
+addsampleFeatures <- function(object=NULL,sample.data=NULL,id.col = "sample"){
+    if(is.null(object)){
+        stop("No object provided")
+    }
     if(!class(object) %in% c("CNQuant","SigQuant")){
         stop("this function requires a CNQuant or SigQuant class object")
     }
@@ -31,11 +34,11 @@ addsampleFeatures <- function(object,sample.data=NULL,id.col = "sample"){
     }
     sampFeat <- object@samplefeatData
     newDataSamples <- sample.data[,which(colnames(sample.data) == id.col)]
-    if(!all(newDataSamples %in% rownames(sampFeat))){
+    if(!any(newDataSamples %in% rownames(sampFeat))){
         stop("no overlapping samples in sample.data")
     }
     mergedsampfeats <- merge.data.frame(sampFeat,sample.data,by.x = "row.names",
-                                        by.y = id.col,all = T)
+                                        by.y = id.col,all.x = T)
     rownames(mergedsampfeats) <- mergedsampfeats$Row.names
     mergedsampfeats <- mergedsampfeats[,-1]
     if(!all(rownames(mergedsampfeats) == names(object@segments))){

diff --git a/R/calculateActivityDrews.R b/R/calculateActivityDrews.R
@@ -2,7 +2,7 @@ calculateActivityDrews = function(object,cancer.subset=NULL) {
 
     # Extract relevant information from object
     V = object@featFitting$sampleByComponent
-    nSamp = length(getSamples(object))
+    nSamp = nrow(object@featFitting$sampleByComponent)
     nFeat = ncol(object@featFitting$sampleByComponent)
 
     # Load signatures

diff --git a/R/calculateFeatures.R b/R/calculateFeatures.R
@@ -3,7 +3,7 @@
 #' @importFrom data.table rbindlist
 setMethod("calculateFeatures",
           signature=c(object="CNQuant"),
-          definition=function(object,method = NULL,smooth.diploid=TRUE,cores=1){
+          definition=function(object,method = NULL,smooth.diploid=TRUE,cores=1,DCIN = 20){
               methods <- c("mac","drews")
               if(is.null(method)){
                   stop("no method provided")
@@ -28,15 +28,33 @@ setMethod("calculateFeatures",
                                                       colNameChr = "chromosome",
                                                       IGNOREDELS = FALSE)
                     # Avoid measurement errors
-                    smoothed = avoidMeasurementErrors(smoothed)
+                    smoothed <- avoidMeasurementErrors(smoothed)
                     # Filter samples without CIN
-                    filtered = removeQuietSamples(smoothed, DCIN = 20)
+                    filtered <- removeQuietSamples(smoothed, DCIN = DCIN)
+
+                    kept <- ifelse(unique(smoothed$sample) %in% unique(filtered$sample),TRUE,FALSE)
+                    if(any(!kept)){
+                        message(cat("Low segment counts: Features not computed for",
+                                    length(unique(smoothed$sample)[!kept]),"of",
+                                    length(unique(smoothed$sample)),
+                                           "samples - see `getSampleFeatures()`"))
+                    }
+                    object <- addsampleFeatures(object = object,
+                                                sample.data = data.frame(sample=unique(smoothed$sample),
+                                                                         computeSigs=kept))
+                    if(nrow(filtered) == 0){
+                        warning("no samples with sufficient copy number alteration counts")
+                        methods::initialize(object,featData=list(),
+                                            ExpData = methods::initialize(object@ExpData,
+                                                                          last.modified = as.character(Sys.time()),
+                                                                          feature.method = method))
+                    } else {
                     # Extract
-                    featData = startCopynumberFeatureExtractionDrews(filtered, cores = cores, RMNORM = TRUE, build=object@ExpData@build)
-                    methods::initialize(object,featData=featData,
+                        featData = startCopynumberFeatureExtractionDrews(filtered, cores = cores, RMNORM = TRUE, build=object@ExpData@build)
+                        methods::initialize(object,featData=featData,
                                         ExpData = methods::initialize(object@ExpData,
                                                                       last.modified = as.character(Sys.time()),
                                                                       feature.method = method))
-
+                    }
                 })
           })
diff --git a/R/calculateSampleByComponentMatrix.R b/R/calculateSampleByComponentMatrix.R
@@ -4,7 +4,7 @@ setMethod("calculateSampleByComponentMatrix",
           signature=c(object="CNQuant"),
           definition=function(object, method=NULL){
               if(length(object@featData) == 0){
-                  stop("Copy number features unavailable - run 'calculateFeatures()'")
+                  stop("Copy number features unavailable (run 'calculateFeatures()) or could not be calculated due to too few segments'")
               }
               if(is.null(method)){
                 method <- getExperiment(object)@feature.method

diff --git a/R/checkInputTypes.R b/R/checkInputTypes.R
@@ -0,0 +1,12 @@
+# Set input types for segtable columns
+checkInputTypes <- function(x=NULL){
+    if(is.null(x)){
+        stop("no table")
+    }
+    x$chromsome <- as.character(x$chromosome)
+    x$start <- as.numeric(x$start)
+    x$end <- as.numeric(x$end)
+    x$segVal <- as.numeric(x$segVal)
+    x$sample <- as.factor(x$sample)
+    return(x)
+}
diff --git a/R/checkbinned.R b/R/checkbinned.R
@@ -1,5 +1,5 @@
 checkbinned <- function(segTable){
-    t.chr <- unique(segTable$chromosome)[1]
+    t.chr <- segTable$chromosome[1]
     t.end <- segTable$end[segTable$chromosome == t.chr]
     t.start <- segTable$start[segTable$chromosome == t.chr]
     startend.len <- length(unique(t.end - t.start))

diff --git a/R/clinPredictionDenovo.R b/R/clinPredictionDenovo.R
@@ -3,28 +3,29 @@
 setMethod("clinPredictionDenovo",
           signature=c(object="SigQuant"),
           definition=function(object, sampTrain, sigsTrain){
-              if(getExperiment(object)@feature.method != "drews"){
-                  stop("This function is only applicable to objects using drews method.")
-              }
-              # Load normalised signature activities
-              mNorm = object@activities$normAct1
-
-              # Extract samples for training
-              if(is.null(sampTrain)) { stop("No sample names supplied.")}
-              if(is.null(sigsTrain)) { stop("No signature names supplied.")}
-              if(length(sigsTrain) != 2) { stop("So far only two signatures can be used.")}
-              mTrain = mNorm[ rownames(mNorm) %in% sampTrain, colnames(mNorm) %in% sigsTrain]
-              mTest = mNorm[ ! rownames(mNorm) %in% sampTrain, colnames(mNorm) %in% sigsTrain]
-
-              # Scale training data and apply to test cohort
-              scaledTrain = scale(mTest)
-              lModel = list(mean = attr(scaledTrain, "scaled:center"),
-                            scale = attr(scaledTrain, "scaled:scale"))
-
-              scaledTest = scaleByModel(mTest, lModel)
-
-              # Do classification
-              vPred = ifelse(scaledTest[,sigsTrain[1]] >= scaledTest[,sigsTrain[2]], paste("Signature", sigsTrain[1], "higher"),
-                             paste("Signature", sigsTrain[2], "higher"))
-              return(vPred)
+              message("this function is depreciated")
+              # if(getExperiment(object)@feature.method != "drews"){
+              #     stop("This function is only applicable to objects using drews method.")
+              # }
+              # # Load normalised signature activities
+              # mNorm = object@activities$normAct1
+              #
+              # # Extract samples for training
+              # if(is.null(sampTrain)) { stop("No sample names supplied.")}
+              # if(is.null(sigsTrain)) { stop("No signature names supplied.")}
+              # if(length(sigsTrain) != 2) { stop("So far only two signatures can be used.")}
+              # mTrain = mNorm[ rownames(mNorm) %in% sampTrain, colnames(mNorm) %in% sigsTrain]
+              # mTest = mNorm[ ! rownames(mNorm) %in% sampTrain, colnames(mNorm) %in% sigsTrain]
+              #
+              # # Scale training data and apply to test cohort
+              # scaledTrain = scale(mTest)
+              # lModel = list(mean = attr(scaledTrain, "scaled:center"),
+              #               scale = attr(scaledTrain, "scaled:scale"))
+              #
+              # scaledTest = scaleByModel(mTest, lModel)
+              #
+              # # Do classification
+              # vPred = ifelse(scaledTest[,sigsTrain[1]] >= scaledTest[,sigsTrain[2]], paste("Signature", sigsTrain[1], "higher"),
+              #                paste("Signature", sigsTrain[2], "higher"))
+              # return(vPred)
 })
diff --git a/R/createCNQuant.R b/R/createCNQuant.R
@@ -47,94 +47,51 @@ createCNQuant <- function(data=NULL,experimentName = "defaultExperiment",build =
         if(file.exists(data)){
             header <- colnames(data.table::fread(input = data,
                                                  header = T,
-                                                 colClasses = c("character",
-                                                                "numeric","numeric",
-                                                                "numeric","character"),
                                                  nrows = 1))
-            if(!any(header == c("chromosome","start","end","segVal","sample"))){
-                stop("Header does not match the required naming")
+            if(!all(header %in% c("chromosome","start","end","segVal","sample"))){
+                stop("Header does not match the required naming ['chromosome','start','end','segVal','sample']")
             }
             segTable <- data.table::fread(input = data,
-                                          header = T,
-                                          colClasses = c("character","numeric",
-                                                         "numeric","numeric",
-                                                         "character"))
-            if(checkSegValRounding(segTable$segVal)){
-                warning("segVal appears to be rounded, copy number signatures
-                        were defined on unrounded absolute copy numbers, use
-                        caution when interpretting and comparing between rounded
-                        and unrounded inputs.")
-            }
-            ## Binned inputs (fixed width not supported yet)
-            # if(checkbinned(segTable)){
-            #     #segTable <- getSegTable()
-            #     #split(segTable,f = as.factor(segTable$sample))
-            # } else {
-            #     segTable <- split(segTable,f = as.factor(segTable$sample))
-            # }
-            ## Temp split until fixed bin input implemented
-            segTable$chromosome <- checkChromosomeFormat(segTable$chromosome)
-            segTable <- dropChromosomes(segTable)
-            segTable <- droplevels(segTable)
-            segTable <- split(segTable,f = as.factor(segTable$sample))
-
-            samplefeatData <- generateSampleFeatData(x = segTable)
-            methods::new("CNQuant",segments = segTable,samplefeatData = samplefeatData,
-                ExpData = methods::new("ExpQuant",
-                              build = build,
-                              samples.full = length(segTable),
-                              samples.current = length(segTable),
-                              experimentName = experimentName))
+                                          header = T)
+            segTable <- segTable[,c("chromosome","start","end","segVal","sample")]
         }
     } else if("QDNAseqCopyNumbers" %in% class(data)){
         segTable <- getSegTable(x = data)
-        if(checkSegValRounding(segTable$segVal)){
-            warning("segVal appears to be rounded, copy number signatures were
-                    defined on unrounded absolute copy numbers, use caution when
-                    interpretting and comparing between rounded and unrounded inputs.")
-        }
-        segTable$chromosome <- checkChromosomeFormat(segTable$chromosome)
-        segTable <- dropChromosomes(segTable)
-        segTable <- droplevels(segTable)
-        segTable <- split(segTable,f = as.factor(segTable$sample))
-        samplefeatData <- generateSampleFeatData(x = segTable)
-        methods::new("CNQuant",segments = segTable,samplefeatData = samplefeatData,
-                     ExpData = methods::new("ExpQuant",
-                          build = build,
-                          samples.full = length(segTable),
-                          samples.current = length(segTable),
-                          experimentName = experimentName))
     } else if(is.data.frame(data)){
         header <- colnames(data)
         if(!all(header %in% c("chromosome","start","end","segVal","sample"))){
             stop("Header does not match the required naming ['chromosome','start','end','segVal','sample']")
         }
         segTable <- data
-        if(checkSegValRounding(segTable$segVal)){
-            warning("segVal appears to be rounded, copy number signatures were
-                    defined on unrounded absolute copy numbers, use caution when
-                    interpretting and comparing between rounded and unrounded inputs.")
-        }
-        ## Binned inputs (fixed width not supported yet)
-        # if(checkbinned(segTable)){
-        #     #segTable <- getSegTable()
-        #     #split(segTable,f = as.factor(segTable$sample))
-        # } else {
-        #     segTable <- split(segTable,f = as.factor(segTable$sample))
-        # }
-        ## Temp split until fixed bin input implemented
-        segTable$chromosome <- checkChromosomeFormat(segTable$chromosome)
-        segTable <- dropChromosomes(segTable)
-        segTable <- droplevels(segTable)
-        segTable <- split(segTable,f = as.factor(segTable$sample))
-        samplefeatData <- generateSampleFeatData(x = segTable)
-        methods::new("CNQuant",segments=segTable,samplefeatData = samplefeatData,
-                     ExpData = methods::new("ExpQuant",
-                          build = build,
-                          samples.full = length(segTable),
-                          samples.current = length(segTable),
-                          experimentName = experimentName))
+        segTable <- segTable[,c("chromosome","start","end","segVal","sample")]
     } else {
         stop("Unknown input format\n")
     }
+    ## Binned inputs (fixed width not supported yet)
+    # if(checkbinned(segTable)){
+    #     #segTable <- getSegTable()
+    #     #split(segTable,f = as.factor(segTable$sample))
+    # } else {
+    #     segTable <- split(segTable,f = as.factor(segTable$sample))
+    # }
+    if(checkSegValRounding(segTable$segVal)){
+        warning("segVal appears to be rounded, copy number signatures
+                        were defined on unrounded absolute copy numbers, use
+                        caution when interpretting and comparing between rounded
+                        and unrounded inputs.")
+    }
+    ## Temp split until fixed bin input implemented
+    segTable <- checkInputTypes(segTable)
+    segTable$chromosome <- checkChromosomeFormat(segTable$chromosome)
+    segTable <- dropChromosomes(segTable)
+    segTable <- droplevels(segTable)
+    segTable <- split(segTable,f = as.factor(segTable$sample))
+    samplefeatData <- generateSampleFeatData(x = segTable)
+    methods::new("CNQuant",segments = segTable,samplefeatData = samplefeatData,
+                 ExpData = methods::new("ExpQuant",
+                                        build = build,
+                                        init.date = as.character(Sys.time()),
+                                        samples.full = length(segTable),
+                                        samples.current = length(segTable),
+                                        experimentName = experimentName))
 }
-Original file line number
+Diff line change
@@ Expand Up / @@ -6,3 +6,4 @@ inst/4_* @@
     /doc/
     /Meta/
     inst/doc
+    tests/testthat/Rplots.pdf