R/getCNfeatures.R

In iC10: A Copy Number and Expression-Based Classifier for Breast Tumours

getCNfeatures <- function (CN, Probes, Map, by.feat, ref, Synonyms)
{
    if (by.feat == "gene") {
        res <- matrix(NA, nrow = length(Probes), ncol = ncol(CN))
        counter <- 1
        no.na <- 0
        for (i in Probes) {
            ids <- which(tolower(rownames(CN)) == tolower(i))
            if (length(ids) > 0) {
                vars <- apply(CN[ids, , drop = F], 1, IQR, na.rm = TRUE)
                selected <- which.max(vars)
                res[counter, ] <- unlist(CN[ids, , drop = F][selected,
                  ])
                no.na <- no.na + 1
            }
            counter <- counter + 1
        }
        if (no.na < length(Probes)) {
            na.ids <- which(apply(res, 1, function(x) mean(is.na(x))) ==
                1)
            for (i in na.ids) {
                one.synon.found <- FALSE
                for (j in Synonyms[which(names(Synonyms) == Probes[i])]) {
                  if (!one.synon.found) {
                    ids <- which(tolower(rownames(CN)) == tolower(j))
                    if (length(ids) > 0) {
                      vars <- apply(CN[ids, , drop = F], 1, IQR,
                        na.rm = TRUE)
                      selected <- which.max(vars)
                      res[i, ] <- unlist(CN[ids, , drop = F][selected,
                        ])
                      no.na <- no.na + 1
                    }
                    one.synon.found <- TRUE
                  }
                }
            }
        }
        rownames(res) <- Probes
        colnames(res) <- colnames(CN)
        new.D1 <- res
        message(paste0("Found ", no.na, " out of ", length(Probes), " copy number features\n"))
    }
    else {
        Map <- Map[which(Map$Probe_ID %in% Probes), ]
        new.D1 <- data.frame(ID = unique(CN$ID))
        for (i in Probes) {
            id <- which(Map$Probe_ID == i)[1]

            sub.CN <- CN[which(CN$chromosome_name == Map[id,
                paste("chromosome_name", ref, sep = "_")] & ((CN$loc.start >=
                Map[id, paste("start_position", ref, sep = "_")] &
                CN$loc.start <= Map[id, paste("end_position",
                  ref, sep = "_")]) | (CN$loc.start <= Map[id,
                paste("start_position", ref, sep = "_")] & CN$loc.end >=
                Map[id, paste("start_position", ref, sep = "_")]))),
                ]
                if (nrow(sub.CN) > 0) {
            feo <- sapply(split(sub.CN, sub.CN$ID), function(x) {
                if (nrow(x) > 0)
                  res <- x$seg.mean[which.max(abs(x$seg.mean))]
                else res <- NA
                res
            })
            }
			else {
				feo <- rep(NA, length(unique(CN$ID)))
				names(feo) <- unique(CN$ID)
			}
            feo <- data.frame(ID = names(feo), feo)
            colnames(feo)[2] <- as.character(Map[id, "Probe_ID"])
            new.D1 <- merge(new.D1, feo, all.x = TRUE)
        }

        if (sum(is.na(new.D1)) > 0) {
            for (i in 1:nrow(new.D1)) {
                for (j in 1:ncol(new.D1)) {
                  if (is.na(new.D1[i, j])) {
                    id <- which(Map$Probe_ID == colnames(new.D1)[j])[1]
                    sub.CN <- CN[which(CN$ID == new.D1[i, 1] &
                      CN$chromosome_name == Map[id, paste("chromosome_name",
                        ref, sep = "_")] & (CN$loc.start >= Map[id,
                      paste("start_position", ref, sep = "_")])),
                      ][1, ]
                    new.D1[i, j] <- sub.CN$seg.mean
                  }
                }
            }
        }
        rownames(new.D1) <- new.D1$ID
        new.D1$ID <- NULL
        new.D1 <- t(new.D1)
    }
    CN <- new.D1
    CN

}