# Unit tests
library(COCOA)
library(data.table)
library(GenomicRanges)
context("Unit tests for COCOA.")
# generate synthetic data
# coordinateDT: positions of cytosines
chr1 <- seq(from=1, to = 4000, by = 200)
chr2 <- seq(from=4001, to = 6400, by = 200)
start <- c(chr1, chr2)
coordinateDT <- data.table(chr=c(rep("chr1", length(chr1)),
rep("chr2", length(chr2))),
start = start)
# loading values
loadingMat <- matrix(data = rep(1, (length(start) * 2)), ncol = 2)
colnames(loadingMat) <- c("PC1", "PC3")
dataDT <- data.table(coordinateDT, loadingMat)
# arbitrarily make some != 1
dataDT$PC1[3] <- 0
dataDT$PC3[5] <- 2
# fake region data (for testing ATAC-seq)
regionCoordDT <- data.table(chr=c(rep("chr1",4), rep("chr2", 2)),
start=c(300, 650, 3000, 10000, 5075, 6150),
end=c(600, 1000, 4000, 11000, 5125, 6300))
regionDataDT <- data.table(PC1=c(1:6), PC2=-1:4)
# region sets that will be tested
regionSet1 <- data.table(chr = c("chr1", "chr1", "chr1", "chr2", "chr2"),
start = c(1, 500, 3100, 5100, 6000),
end = c(400, 700, 3400, 5150, 6450))
regionSet1 <- MIRA:::dtToGr(regionSet1)
regionSet2 <- data.table(chr = c("chr1", "chr1", "chr1", "chr2", "chr2"),
start = c(1000, 1500, 3700, 4000, 4300),
end = c(1400, 1700, 3800, 4100, 4700))
regionSet2 <- MIRA:::dtToGr(regionSet2)
# don't use built in package data for tests because it will take longer
# running the tests
# aggregateSignal(signal = loadingMat,
# coordinateDT = coordinates,
# regionSet = regionSet)
# signalOLMetrics
dataDT <- cbind(coordinateDT, as.data.frame(loadingMat))
# making test data for Wilcoxon rank sum test
chr3 <- seq(from=100, to = 700, by = 100)
coordinateDTW <- data.table(chr=rep("chr3", length(chr3)),
start = chr3, end = chr3)
regionSetW <- data.table(chr = c("chr3", "chr3"),
start = c(50, 250),
end = c(150, 450))
regionSetW <- COCOA:::dtToGr(regionSetW)
loadingMatW <- matrix(data = c(-2:4, seq(from=8, to=2, by=-1)), ncol = 2)
colnames(loadingMatW) <- c("PC2", "PC3")
loadingMatW[7, "PC3"] <- 10
dataDTW <- cbind(coordinateDTW, as.data.frame(loadingMatW))
test_that("aggregateSignal, scoring metrics, and aggregateSignalGRList", {
# # test wilcoxon rank sum scoring metric
rsWResults <- COCOA:::rsWilcox(dataDT = dataDTW, regionSet = regionSetW)
PC2W <- wilcox.test(x = c(-2, 0, 1), y=c(-1, 2:4))$p.value
PC3W <- wilcox.test(x = c(8, 6, 5), y=c(7, 4, 3, 10))$p.value
expect_equal(c(PC2W, PC3W, 3, 2, 2, mean(width(regionSetW))),
c(rsWResults$PC2, rsWResults$PC3, rsWResults$signalCoverage,
rsWResults$regionSetCoverage,
rsWResults$totalRegionNumber,
rsWResults$meanRegionSize))
# # same test for Wilcoxon but with aggregateSignal (absolute value
# # of loadings will be taken), "greater" alternate hypothesis is used in
# # aggregateSignal
# rsWResults <- COCOA:::aggregateSignal(signal = loadingMatW,
# signalCoord = coordinateDTW,
# regionSet = regionSetW,
# signalCol = c("PC2", "PC3"),
# scoringMetric = "rankSum")
# PC2W <- wilcox.test(x = c(2, 0, 1), y=c(1, 2:4), alternative = "greater")$p.value
# PC3W <- wilcox.test(x = c(8, 6, 5), y=c(7, 4, 3, 10), alternative = "greater")$p.value
# expect_equal(c(PC2W, PC3W, 3, 2, 2, mean(width(regionSetW))),
# c(rsWResults$PC2, rsWResults$PC3, rsWResults$signalCoverage,
# rsWResults$regionSetCoverage,
# rsWResults$totalRegionNumber,
# rsWResults$meanRegionSize))
# test regionMean scoring method, average first within region, then
# between regions
regionMeanRes <- COCOA:::aggregateSignal(signal = loadingMatW,
signalCoord = coordinateDTW,
regionSet = regionSetW,
signalCol = c("PC2", "PC3"),
scoringMetric = "regionMean")
PC2R <- mean(c(2, mean(c(0, 1))))
PC3R <- mean(c(8, mean(c(6, 5))))
expect_equal(c(PC2R, PC3R, 3, 2, 2, mean(width(regionSetW))),
c(regionMeanRes$PC2, regionMeanRes$PC3, regionMeanRes$signalCoverage,
regionMeanRes$regionSetCoverage,
regionMeanRes$totalRegionNumber,
regionMeanRes$meanRegionSize))
########### test simpleMean scoring method ################
# this is a mean of CpG loading values just like "raw" but instead
# of averaging within regions then averaging regions together, this
# method does the simple average of all CpGs within the region set
simpleMeanRes <- COCOA:::aggregateSignal(signal = loadingMatW,
signalCoord = coordinateDTW,
regionSet = regionSetW,
signalCol = c("PC2", "PC3"),
scoringMetric = "simpleMean")
PC2RC <- mean(c(2, 0, 1))
PC3RC <- mean(c(8, 6, 5))
expect_equal(c(PC2RC, PC3RC, 3, 2, 2, mean(width(regionSetW))),
c(simpleMeanRes$PC2, simpleMeanRes$PC3, simpleMeanRes$signalCoverage,
simpleMeanRes$regionSetCoverage,
simpleMeanRes$totalRegionNumber,
simpleMeanRes$meanRegionSize))
###########################################################
# test simpleMedian scoring method
simpleMedRes <- COCOA:::aggregateSignal(signal = loadingMatW,
signalCoord = coordinateDTW,
regionSet = regionSetW,
signalCol = c("PC2", "PC3"),
scoringMetric = "simpleMedian")
PC2RC <- median(c(2, 0, 1))
PC3RC <- median(c(8, 6, 5))
expect_equal(c(PC2RC, PC3RC, 3, 2, 2, mean(width(regionSetW))),
c(simpleMedRes$PC2, simpleMedRes$PC3, simpleMedRes$signalCoverage,
simpleMedRes$regionSetCoverage,
simpleMedRes$totalRegionNumber,
simpleMedRes$meanRegionSize))
##########################################################
# test regionMedian scoring method
# median first within region, then
# between regions
tmpCoordDTW = data.frame(chr=c(coordinateDTW$chr, "chr3"),
start = c(coordinateDTW$start, 350),
end = c(coordinateDTW$end, 350))
tmpLoadMatW = data.frame(PC2= c(loadingMatW[, "PC2"], 20),
PC3=c(loadingMatW[, "PC3"], 85))
regionMedRes <- COCOA:::aggregateSignal(signal = tmpLoadMatW,
signalCoord = tmpCoordDTW,
regionSet = regionSetW,
signalCol = c("PC2", "PC3"),
scoringMetric = "regionMedian")
PC2R <- median(c(2, median(c(0, 1, 20))))
PC3R <- median(c(8, median(c(6, 5, 85))))
expect_equal(c(PC2R, PC3R, 4, 2, 2, mean(width(regionSetW))),
c(regionMedRes$PC2, regionMedRes$PC3, regionMedRes$signalCoverage,
regionMedRes$regionSetCoverage,
regionMedRes$totalRegionNumber,
regionMedRes$meanRegionSize))
#########################################################
# # test mean difference scoring method
# PC2Num <- mean(c(2, 0, 1)) - mean(c(1, 2, 3, 4))
# PC3Num <- mean(c(8, 6, 5)) - mean(c(7, 4, 3, 10))
# # see doi: 10.1186/s13040-015-0059-z
# # pooled variance times normalization factor for size of each set
# PC2Denom <- sqrt((sd(c(2, 0, 1))^2 + sd(c(1, 2, 3, 4))^2) / 2) * sqrt(1/3 - 1/4)
# PC3Denom <- sqrt((sd(c(8, 6, 5))^2 + sd(c(7, 4, 3, 10))^2) / 2) * sqrt(1/3 - 1/4)
#
# PC2MD <- PC2Num / PC2Denom
# PC3MD <- PC3Num / PC3Denom
# mdRes <- COCOA:::aggregateSignal(signal = loadingMatW,
# signalCoord = coordinateDTW,
# regionSet = regionSetW,
# signalCol = c("PC2", "PC3"),
# scoringMetric = "meanDiff")
# expect_equal(c(PC2MD, PC3MD, 3, 2, 2, mean(width(regionSetW))),
# c(mdRes$PC2, mdRes$PC3, mdRes$signalCoverage,
# mdRes$regionSetCoverage,
# mdRes$totalRegionNumber,
# mdRes$meanRegionSize))
################## test aggregateSignalGRList with meanDiff test data
coordinateDTW2 <- copy(coordinateDTW)
coordinateDTW2$end <- coordinateDTW$end
# make sure it works if there is metadata in GRanges object
# extraCol <- rep(0, nrow(coordinateDTW2))
# signalCoordW2 <- COCOA:::dtToGr(coordinateDTW2)
# mcols(signalCoordW2) <- data.frame(extraCol)
# twoResults <- aggregateSignalGRList(signal = loadingMatW,
# signalCoord = signalCoordW2,
# GRList = GRangesList(regionSetW, regionSetW),
# signalCol = c("PC2", "PC3"),
# scoringMetric = "meanDiff")
# expect_equal(rep(c(PC2MD, PC3MD, 3, 2, 2, mean(width(regionSetW))), each=2),
# c(twoResults$PC2, twoResults$PC3, twoResults$signalCoverage,
# twoResults$regionSetCoverage,
# twoResults$totalRegionNumber,
# twoResults$meanRegionSize))
########## testing signalOLMetrics, used for meanDiff scoring method #######
olMetrics <- COCOA:::signalOLMetrics(dataDT=dataDTW,
dataGR = dtToGr(dataDTW),
regionSet = regionSetW,
metrics = c("mean", "sd"),
alsoNonOLMet = TRUE)
PC2Met <- data.table(t(c(mean(c(-2, 0, 1)), sd(c(-2, 0, 1)),
mean(c(-1, 2:4)), sd(c(-1, 2:4)), 3, 2, 2, 151)))
PC2Met <- cbind("PC2", PC2Met)
colnames(PC2Met) <- colnames(olMetrics)
expect_equal(olMetrics[1, ], PC2Met)
PC3Met <- data.table(t(c(mean(c(8, 6, 5)), sd(c(8, 6, 5)),
mean(c(7, 4, 3, 10)), sd(c(7, 4, 3, 10)), 3, 2, 2, 151)))
PC3Met <- cbind("PC3", PC3Met)
colnames(PC3Met) <- colnames(olMetrics)
expect_equal(olMetrics[2, ], PC3Met)
})
test_that("ATAC-seq scoring methods", {
# test "regionOLWeightedMean"
weightedAve <- regionOLWeightedMean(signalMat = regionDataDT, signalGR = COCOA:::dtToGr(regionCoordDT),
regionSet = regionSet1, calcCols = c("PC1", "PC2"))
# proportion overlap is first then PC
correctAve <- data.frame(PC1=((101/400)*1+(101/201)*1+(51/201)*2+1*3+(26/51)*5+(151/451)*6) /
((101/400)+(101/201)+(51/201)+1+(26/51)+(151/451)),
PC2=((101/400)*-1+(101/201)*-1+(51/201)*0+1*1+(26/51)*3+(151/451)*4) /
((101/400)+(101/201)+(51/201)+1+(26/51)+(151/451)),
signalCoverage=5, regionSetCoverage=5,
sumProportionOverlap=(101/400)+(101/201)+(51/201)+1+(26/51)+(151/451))
expect_equal(weightedAve, correctAve)
# test "regionOLMean"
signalAve <- regionOLMean(signalDT = regionDataDT,
signalGR = COCOA:::dtToGr(regionCoordDT),
regionSet = regionSet1,
calcCols = c("PC1", "PC2"))
# proportion overlap is first then PC
correctAve <- data.frame(PC1=(1+1+2+3+5+6)/6,
PC2=(-1+-1+0+1+3+4)/6,
signalCoverage=5,
regionSetCoverage=5)
expect_equal(signalAve, correctAve)
# test "weightedAvePerRegion"
avePerRegion <- weightedAvePerRegion(signalDT= regionDataDT,
signalCoord=COCOA:::dtToGr(regionCoordDT),
regionSet=regionSet1,
calcCols = c("PC1", "PC2"))
correctAve <- data.table(PC1 = c(1*1, (101/201*1 + 51/201*2)/ (101/201 + 51/201), 1*3, 1*5, 1*6),
PC2 = c(1*-1, (101/201*-1 + 51/201*0)/ (101/201 + 51/201), 1*1, 1*3, 1*4))
expect_equal(avePerRegion$PC1, correctAve$PC1)
expect_equal(avePerRegion$PC2, correctAve$PC2)
})
test_that("averagePerRegion", {
loadingMatABR <- loadingMat
loadingMatABR[1, "PC1"] <- 3
loadingMatABR[32, "PC1"] <- 2
abr <- COCOA:::averagePerRegion(signal = loadingMatABR, signalCoord = COCOA:::dtToGr(coordinateDT),
regionSet = regionSet1, signalCol = "PC1",
returnQuantile = FALSE)
expect_equal(abr$PC1, c(2, 1, 1, 1.5))
# test quantile
abrq <- COCOA:::averagePerRegion(signal = loadingMatABR, signalCoord = COCOA:::dtToGr(coordinateDT),
regionSet = regionSet1, signalCol = "PC1",
returnQuantile = TRUE)
# the mean values, abr$PC1, converted to quantiles
# converting properly to quantiles?
expect_equal(abrq$PC1, ecdf(loadingMatABR[, "PC1"])(abr$PC1))
})
test_that("getMetaRegionProfile", {
loadingMatP <- matrix(data = c(1:8, 10, 10, rep(1, 20)), nrow = 30)
colnames(loadingMatP) <- "PC1"
# two CpGs per bin of regionSetP (with 5 bins)
.start <- c(seq(from=50, to=950, by = 100),
seq(from=2050, to=2950, by = 100),
seq(from=4050, to=4950, by = 100))
coordinateDTP <- data.frame(chr = rep("chr3", nrow(loadingMatP)),
start = .start,
end = .start,
extraCol = rep(1, length(.start)))
regionSetP <- data.table(chr = rep("chr3", 3),
start = c(1, 2001, 4001),
end = c(1010, 3010, 5010))
regionSetP <- COCOA:::dtToGr(regionSetP)
# COCOA:::dtToGr(coordinateDTP)
binnedP <- lapply(GRangesList(regionSetP),
function(x) getMetaRegionProfile(signal = loadingMatP,
signalCoord = coordinateDTP,
regionSet = x,
signalCol = "PC1",
binNum = 5))
meanPerBin <- (c(seq(from=1.5, to=7.5, by=2), 10) + rep(1, 5) + rep(1, 5)) / 3
# BSBinAggregate averages the profile around the center
symmetricalBin <- (meanPerBin + rev(meanPerBin)) / 2
expect_equal(binnedP[[1]]$PC1, symmetricalBin)
# making sure that different input formats still work
#########
alterOut <- lapply(GRangesList(regionSetP),
function(x) getMetaRegionProfile(signal = data.frame(loadingMatP),
signalCoord = COCOA:::dtToGr(coordinateDTP),
regionSet = x,
signalCol = "PC1",
binNum = 5))
expect_equal(binnedP, alterOut)
##########################################################
# # test multiBase data like ATAC-seq
# coordinateDTP <- data.frame(chr = rep("chr3", nrow(loadingMatP)),
# start = .start,
# end = .start+1,
# extraCol = rep(1, length(.start)))
# # make 1 region overlap multiple bins ( % first bin, % second bin)
# # 150 to
# coordinateDTP$end[2] = 225
# meanPerBin <- (c(seq(from=1.5, to=7.5, by=2), 10) * (x, x, 1, 1, 1) +
# rep(1, 5) + rep(1, 5)) / 3
#
})
################### test function inputs
# make sure various input formats work
#### take out these high level tests in order to save time
test_that("Input formats", {
normalOut <- aggregateSignalGRList(signal = loadingMatW,
signalCoord = coordinateDTW,
GRList = GRangesList(regionSetW),
signalCol = "PC2",
scoringMetric = "regionMean")
alterOut <- aggregateSignalGRList(signal = data.frame(loadingMatW),
signalCoord = COCOA:::dtToGr(coordinateDTW),
GRList = regionSetW,
signalCol = "PC2",
scoringMetric = "regionMean")
expect_equal(normalOut, alterOut)
})
