#' @name compFGA
#' @title Comparison of fraction genome altered
#' @description This function calculates Fraction Genome Altered (FGA), Fraction Genome Gained (FGG), and Fraction Genome Lost (FGL) seperately, and compares them among curent subtypes identified from multi-omics integrative clustering algorithms.
#' @param moic.res An object returned by `getMOIC()` with one specified algorithm or `get\%algorithm_name\%` or `getConsensusMOIC()` with a list of multiple algorithms.
#' @param segment A data frame containing segmented copy number and columns must exactly include the following elements: c('sample','chrom','start','end','value'). Column of `value` should be segments value when \code{iscopynumber = FALSE} but copy-number value when \code{iscopynumber = TRUE}. Copy-number will be converted to segments by log2(copy-number/2).
#' @param iscopynumber A logical value to indicate if the fifth column of segment input is copy-number. If segment file derived from CNV calling provides copy number instead of segment_mean value, this argument must be switched to TRUE. FALSE by default.
#' @param cnathreshold A numeric value to indicate the cutoff for identifying copy-number gain or loss. 0.2 by default.
#' @param test.method A string value to indicate the method for statistical testing. Allowed values contain c('nonparametric', 'parametric'); nonparametric means two-sample wilcoxon rank sum test for two subtypes and Kruskal-Wallis rank sum test for multiple subtypes; parametric means two-sample t-test when only two subtypes are identified, and anova for multiple subtypes comparison; "nonparametric" by default.
#' @param barcolor A string vector to indicate the mapping color for bars of FGA, FGG and FGL.
#' @param clust.col A string vector storing colors for each subtype.
#' @param fig.path A string value to indicate the output path for storing the barplot.
#' @param fig.name A string value to indicate the name of the barplot.
#' @param width A numeric value to indicate the width of barplot.
#' @param height A numeric value to indicate the height of barplot.
#' @return A list contains the following components:
#'
#' \code{summary} a table summarizing the measurements of FGA, FGG, and FGL per sample
#'
#' \code{FGA.p.value} a nominal p value quantifying the difference of FGA among current subtypes
#'
#' \code{pairwise.FGA.test} a pairwise BH adjusted p value matrix for multiple comparisons of FGA if more than 2 subtypes were identified
#'
#' \code{FGG.p.value} a nominal p value quantifying the difference of FGG among current subtypes
#'
#' \code{pairwise.FGG.test} a pairwise BH adjusted p value matrix for multiple comparisons of FGG if more than 2 subtypes were identified
#'
#' \code{FGL.p.value} a nominal p value quantifying the difference of FGL among current subtypes
#'
#' \code{pairwise.FGL.test} a pairwise BH adjusted p value matrix for multiple comparisons of FGL if more than 2 subtypes were identified
#'
#' \code{test.method} a string value indicating the statistical testing method to calculate p values
#'
#' @export
#' @import ggplot2
#' @import patchwork
#' @importFrom dplyr group_by summarize %>%
#' @importFrom ggpubr stat_compare_means
#' @examples # There is no example and please refer to vignette.
#' @references Cerami E, Gao J, Dogrusoz U, et al. (2012). The cBio Cancer Genomics Portal: An Open Platform for Exploring Multidimensional Cancer Genomics Data. Cancer Discov, 2(5):401-404.
#'
#' Gao J, Aksoy B A, Dogrusoz U, et al. (2013). Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal, 6(269):pl1-pl1.
compFGA <- function(moic.res = NULL,
segment = NULL,
iscopynumber = FALSE,
cnathreshold = 0.2,
test.method = "nonparametric",
barcolor = c("#008B8A", "#F2042C", "#21498D"),
clust.col = c("#2EC4B6","#E71D36","#FF9F1C","#BDD5EA","#FFA5AB","#011627","#023E8A","#9D4EDD"),
fig.path = getwd(),
fig.name = NULL,
width = 8,
height = 4) {
# check arguments
if(!all(is.element(c("sample","chrom","start","end","value"), colnames(segment)))) {
stop("segment data must have the following columns: sample, chrom, start, end, value.")
}
if(iscopynumber) {
segment$value <- log2(segment$value/2) # convert copy-number to segment-mean value
}
comsam <- intersect(moic.res$clust.res$samID, unique(segment[,1]))
# check data
if(length(comsam) == nrow(moic.res$clust.res)) {
message("--all samples matched.")
} else {
message(paste0("--",(nrow(moic.res$clust.res)-length(comsam))," samples mismatched from current subtypes."))
}
if(!is.element(test.method, c("nonparametric","parametric"))) {
stop("test.method can be one of nonparametric or parametric.")
}
clust.res <- moic.res$clust.res[comsam, , drop = FALSE]
segment <- segment[which(segment$sample %in% comsam),]
n.moic <- length(unique(clust.res$clust))
# data process
segment$bases <- segment$end - segment$start
# calculate FGA, FGG and FGL
display.progress = function (index, totalN, breakN=20) {
if ( index %% ceiling(totalN/breakN) == 0 ) {
cat(paste(round(index*100/totalN), "% ", sep = ""))
}
}
std <- function(x, na.rm = TRUE) {
if(na.rm) {
x <- as.numeric(na.omit(x))
sd(x)/sqrt(length(x))
} else {sd(x)/sqrt(length(x))}
}
outTab <- data.frame()
for (i in 1:length(unique(segment$sample))) {
display.progress(index = i, totalN = length(unique(segment$sample)))
tmp <- segment[segment$sample == names(table(segment$sample))[i],]
if (length(tmp[abs(tmp$value) > cnathreshold,"bases"][6]) == 0) {
FGA = 0
} else {
FGA = sum(tmp[abs(tmp$value) > cnathreshold,"bases"]) / sum(tmp[,"bases"])
}
if (length(tmp[tmp$value > cnathreshold,"bases"][6]) == 0) {
FGG = 0
} else {
FGG = sum(tmp[tmp$value > cnathreshold,"bases"]) / sum(tmp[,"bases"])
}
if (length(tmp[tmp$value < (-cnathreshold),"bases"][6]) == 0) {
FGL = 0
} else {
FGL = sum(tmp[tmp$value < (-cnathreshold),"bases"]) / sum(tmp[,"bases"])
}
tmp <- data.frame(samID = names(table(segment$sample))[i],
FGA = FGA,
FGG = FGG,
FGL = FGL,
stringsAsFactors = FALSE)
outTab <- rbind.data.frame(outTab, tmp, stringsAsFactors = FALSE)
}
outTab$Subtype <- paste0("CS", clust.res[outTab$samID, "clust"])
# calculate mean and se
summaryFGA <- outTab %>% group_by(Subtype) %>% dplyr::summarize(mean = mean(FGA, na.rm = TRUE), se = std(FGA, na.rm = TRUE))
summaryFGG <- outTab %>% group_by(Subtype) %>% dplyr::summarize(mean = mean(FGG, na.rm = TRUE), se = std(FGG, na.rm = TRUE))
summaryFGL <- outTab %>% group_by(Subtype) %>% dplyr::summarize(mean = mean(FGL, na.rm = TRUE), se = std(FGL, na.rm = TRUE))
summaryFGGL <- data.frame(rbind.data.frame(summaryFGG,summaryFGL),class = rep(c("FGG","FGL"),c(nrow(summaryFGG),nrow(summaryFGL))),stringsAsFactors = FALSE)
# statistical testing
# generate boxviolin plot with statistical testing
if(n.moic == 2 & test.method == "nonparametric") {
statistic <- "wilcox.test"
FGA.test <- wilcox.test(outTab$FGA ~ outTab$Subtype)$p.value
FGG.test <- wilcox.test(outTab$FGG ~ outTab$Subtype)$p.value
FGL.test <- wilcox.test(outTab$FGL ~ outTab$Subtype)$p.value
}
if(n.moic == 2 & test.method == "parametric") {
statistic <- "t.test"
FGA.test <- t.test(outTab$FGA ~ outTab$Subtype)$p.value
FGG.test <- t.test(outTab$FGG ~ outTab$Subtype)$p.value
FGL.test <- t.test(outTab$FGL ~ outTab$Subtype)$p.value
}
if(n.moic > 2 & test.method == "nonparametric") {
statistic <- "kruskal.test"
FGA.test <- kruskal.test(outTab$FGA ~ outTab$Subtype)$p.value
FGG.test <- kruskal.test(outTab$FGG ~ outTab$Subtype)$p.value
FGL.test <- kruskal.test(outTab$FGL ~ outTab$Subtype)$p.value
pairwise.FGA.test <- pairwise.wilcox.test(outTab$FGA,outTab$Subtype,p.adjust.method = "BH")
pairwise.FGG.test <- pairwise.wilcox.test(outTab$FGG,outTab$Subtype,p.adjust.method = "BH")
pairwise.FGL.test <- pairwise.wilcox.test(outTab$FGL,outTab$Subtype,p.adjust.method = "BH")
}
if(n.moic > 2 & test.method == "parametric") {
statistic <- "anova"
FGA.test <- summary(aov(outTab$FGA ~ outTab$Subtype))[[1]][["Pr(>F)"]][1]
FGG.test <- summary(aov(outTab$FGG ~ outTab$Subtype))[[1]][["Pr(>F)"]][1]
FGL.test <- summary(aov(outTab$FGL ~ outTab$Subtype))[[1]][["Pr(>F)"]][1]
pairwise.FGA.test <- pairwise.t.test(outTab$FGA,outTab$Subtype,p.adjust.method = "BH")
pairwise.FGG.test <- pairwise.t.test(outTab$FGG,outTab$Subtype,p.adjust.method = "BH")
pairwise.FGL.test <- pairwise.t.test(outTab$FGL,outTab$Subtype,p.adjust.method = "BH")
}
# generate barplot
FGA.col <- barcolor[1]
FGG.col <- barcolor[2]
FGL.col <- barcolor[3]
p1 <- ggplot(summaryFGA, aes(x = Subtype, y = mean,fill=rep("0",nrow(summaryFGA)))) +
geom_bar(stat = 'identity') +
geom_errorbar(aes(ymax = mean+se, ymin = mean-se),position = position_dodge(0.9), width = 0.15) +
annotate(geom="text",
#hjust = ifelse(n.moic %% 2 == 0, 0.5, 0),
x = n.moic / 2 + 0.5,
#y = as.numeric(summaryFGA[which.max(summaryFGA$mean),"mean"] + summaryFGA[which.max(summaryFGA$mean),"se"]),
y = as.numeric(summaryFGA[which.max(summaryFGA$mean),"mean"]),
size = 8, angle = 90, fontface = "bold",
#label = paste0(statistic, " p = ", formatC(FGA.test, digits = 1, format = "e")),
label = cut(FGA.test,c(0,0.001,0.01,0.05,0.1,1),labels = c("****","***","**","*","."))) +
scale_x_discrete(name = "",position = "top") +
theme_bw() +
theme(axis.line.y = element_line(size = 0.8),
axis.ticks.y = element_line(size = 0.2),
axis.text.y = element_blank(),
axis.title.x = element_text(vjust = -0.3,size = 12),
axis.text.x = element_text(size = 10, color = "black"),
plot.margin = unit(c(0.3, -1.7, 0.3, 0.3), "lines"),
legend.title = element_blank()) +
coord_flip() +
scale_fill_manual(values = FGA.col, breaks = c("0"), labels = c("Copy number-altered genome")) +
scale_y_reverse(expand = c(0.01,0),
name = "FGA (Fraction of Genome Altered)", position = "left")
p2 <- ggplot(summaryFGGL, aes(x = Subtype, y = ifelse(class == 'FGG',mean,-mean),fill=class)) +
geom_bar(stat = 'identity') +
geom_errorbar(data=summaryFGGL[summaryFGGL$class=='FGG',],aes(ymax = mean+se, ymin =mean-se),position = position_dodge(0.9), width = 0.15) +
geom_errorbar(data=summaryFGGL[summaryFGGL$class=='FGL',],aes(ymax = -mean-se, ymin =-mean+se),position = position_dodge(0.9), width = 0.15) +
annotate(geom="text",
x = n.moic / 2 + 0.5,
y = -as.numeric(summaryFGL[which.max(summaryFGL$mean),"mean"]),
size = 8, angle = 90, fontface = "bold",
label = cut(FGL.test,c(0,0.001,0.01,0.05,0.1,1),labels = c("****","***","**","*","."))) +
annotate(geom="text",
x = n.moic / 2 + 0.5,
y = as.numeric(summaryFGG[which.max(summaryFGG$mean),"mean"]),
size = 8, angle = 90, fontface = "bold",
label = cut(FGG.test,c(0,0.001,0.01,0.05,0.1,1),labels = c("****","***","**","*","."))) +
scale_x_discrete(name = "") +
theme_bw() +
theme(axis.line.y = element_line(size = 0.8),
axis.ticks.y = element_line(size = 0.2),
axis.text.y = element_blank(),
axis.title.x = element_text(vjust = -0.3, size = 12),
axis.text.x = element_text(size = 10, color = "black"),
plot.margin = unit(c(0.3, 0.3, 0.3, -1), "lines"),
legend.title = element_blank()) +
coord_flip() +
scale_fill_manual(values = c(FGL.col, FGG.col), breaks = c("FGL","FGG"),
labels = c("Copy number-lost genome","Copy number-gained genome")) +
scale_y_continuous(expand = c(0.01,0),
name = "FGL or FGG (Fraction of Genome Lost or Gained)")
pp <- ggplot() +
# geom_text(data = summaryFGGL,
# aes(label = Subtype, x=Subtype), y = 0.5,
# size = 0.8*11/.pt, # match font size to theme
# hjust = 0.5, vjust = 0.5) +
geom_label(data = summaryFGGL,
aes(label = Subtype, x = Subtype, fill = Subtype),
y = 0.5,
color = "white",
size = 0.9*11/.pt, # match font size to theme
hjust = 0.4, vjust = 0.5) +
scale_fill_manual(values = clust.col) +
theme_minimal()+
theme(axis.line.y =element_blank(),
axis.ticks.y =element_blank(),
axis.text.y =element_blank(),
axis.title.y =element_blank(),
axis.title.x =element_blank(),
plot.margin = unit(c(0.3, 0, 0.3, 0), "lines")
) +
guides(fill = FALSE) +
coord_flip() +
scale_y_reverse()
pal <- p1 + pp + p2 +
plot_layout(widths = c(7,1,7), guides = 'collect') & theme(legend.position = 'top')
# save to pdf
if(is.null(fig.name)) {
outFig <- "barplot of FGA.pdf"
} else {
outFig <- paste0(fig.name,".pdf")
}
ggsave(file.path(fig.path, outFig), width = width, height = height)
# print to screen
print(pal)
if(n.moic > 2) {
return(list(summary = outTab,
FGA.p.value = FGA.test,
pairwise.FGA.test = pairwise.FGA.test,
FGG.p.value = FGG.test,
pairwise.FGG.test = pairwise.FGG.test,
FGL.p.value = FGL.test,
pairwise.FGL.test = pairwise.FGL.test,
test.method = statistic))
} else {
return(list(summary = outTab,
FGA.p.value = FGA.test,
FGG.p.value = FGG.test,
FGL.p.value = FGL.test,
test.method = statistic))
}
}
Datasets
Models
