Nothing
Introduction"
In ROKET: Optimal Transport-Based Kernel Regression
knitr::opts_chunk$set(
collapse = TRUE,comment = "#>",
echo = TRUE,cache = TRUE,
dev = "png")
ot_eval = TRUE
\def\T{\text{T}}
\newcommand{\bf}[1]{\mathbf{#1}}
\newcommand{\bigPar}[1]{\left(#1\right)}
\newcommand{\bigCur}[1]{\left{#1\right}}
\newcommand{\bcSqu}[2]{\left[#1 \middle| #2\right]}
\newcommand{\cE}[2]{E\bcSqu{#1}{#2}}
\newcommand{\nexp}[1]{\exp\bigCur{#1}}
\newcommand{\ind}[1]{1\bigCur{#1}}
\newcommand{\w}[1]{\widehat{#1}}
Overview
Assuming all software dependencies and ROKET [@little2023associating] installation are installed, we can begin.
# Load libraries
req_packs = c("devtools","smarter","ggplot2","reshape2",
"survival","ggdendro","MiRKAT","ROKET")
for(pack in req_packs){
library(package = pack,character.only = TRUE)
}
# List package's exported functions
ls("package:ROKET")
# Fix seed
set.seed(2)
Distance Motivation
For $i = 1,\ldots,N$ and $g = 1,\ldots,G$, let
- $Z_{ig} = 1$ indicate the $i$th sample's $g$th gene is mutated and $Z_{ig} = 0$ otherwise
- $\bf{Z}i \equiv \bigPar{Z{i1},\ldots,Z_{iG}}^\T$
We would like to calculate the distance between the $i$th and $j$th samples in terms of mutated genes $\bf{Z}_i$ and $\bf{Z}_j$, denoted by $d\bigPar{\bf{Z}_i,\bf{Z}_j}$.
Optimal Transport
Background papers
Several optimal transport (OT) references to get familiarized with the framework, objective functions, updating equations, entropic regularizations, types of OT are provided:
- @cuturi2013sinkhorn,
- @zhang2021review,
- @jagarlapudi2020statistical,
- @fatras2021unbalanced,
- @chapel2021unbalanced
Basic Idea
For the $i$th and $j$th samples, let
- $\bf{P}^{(ij)} =$ an unknown $G$ by $G$ transport matrix to be estimated
- $P_{gh}^{(ij)} =$ the value of the $g$th row and $h$th column of $\bf{P}$, denotes the mass transported between the $g$th gene of the $i$th sample and the $h$th gene of the $j$th sample
- $\bf{W} =$ a known $G$ by $G$ cost matrix, independent of samples
- $W_{gh} =$ the value of the $g$th row and $h$th column of $\bf{W}$, denotes the cost to transport one unit of mass between the $g$th gene and $h$th gene
- The space of transport matrices $\bf{P}^{(ij)}$ to search over is governed through the penalized divergence between $\bf{Z}_i$ and the row sums of $\bf{P}^{(ij)}$ as well as between $\bf{Z}_j$ and the column sums of $\bf{P}^{(ij)}$
- $\w{\bf{P}^{(ij)}} = \text{argmin}{\bf{P}^{(ij)}} \bigPar{\sum_g \sum_h P{gh}^{(ij)} W_{gh}}$, the optimal transport matrix
- $d\bigPar{\bf{Z}i,\bf{Z}_j} = \sum_g \sum_h \w{P{gh}^{(ij)}} W_{gh}$
Balanced OT
If the mass of the two vectors are equal or normalized such that $\sum_g Z_{ig} = \sum_g Z_{jg}$, we could use balanced optimal transport.
Unbalanced OT
If the mass of the two vectors are not equal, $\sum_g Z_{ig} \neq \sum_g Z_{jg}$, we could use unbalanced optimal transport with penalty parameters.
Simulated Example
The code below will simulate samples and mutated genes. We welcome the reader to manipulate the following input arguments:
NN for sample size,PP for number of pathways,GG for number of genes, and bnd_same the lower bound gene similarity for genes sharing the same pathway and 1 - bnd_same, an upper bound on gene similarity for genes not sharing the same pathway
Ideally, PP should be much less than GG to allow multiple genes to be allocated to each pathway.
# number of samples
NN = 30
NN_nms = sprintf("S%s",seq(NN))
# number of pathways
PP = 4
PP_nms = sprintf("P%s",seq(PP))
# number of genes
GG = 30
GG_nms = sprintf("G%s",seq(GG))
# bound for gene similarity of two genes on same or different pathway
bnd_same = 0.75
# Gene and pathway relationship
GP = smart_df(PATH = sample(seq(PP),GG,replace = TRUE),
GENE = seq(GG))
table(GP$PATH)
# gene-gene similarity matrix
GS = matrix(NA,GG,GG)
dimnames(GS) = list(GG_nms,GG_nms)
diag(GS) = 1
tmp_mat = t(combn(seq(GG),2))
for(ii in seq(nrow(tmp_mat))){
G1 = tmp_mat[ii,1]
G2 = tmp_mat[ii,2]
same = GP$PATH[GP$GENE == G1] == GP$PATH[GP$GENE == G2]
if( same )
GS[G1,G2] = runif(1,bnd_same,1)
else
GS[G1,G2] = runif(1,0,1 - bnd_same)
}
GS[lower.tri(GS)] = t(GS)[lower.tri(GS)]
# round(GS,3)
Gene-Gene Similarity
Let's take a look at the gene similarity matrix.
show_tile = function(MAT,LABEL,TYPE = NULL,
LABx = NULL,LABy = NULL,DIGITS = 1){
min_val = min(MAT)
max_val = max(MAT)
med_val = (min_val + max_val) / 2
if( is.null(TYPE) ) stop("Specify TYPE")
if( isSymmetric(MAT) )
MAT[upper.tri(MAT,diag = !TRUE)] = NA
if( TYPE == "GSIM" ){
max_val = 1
med_val = 0.5
}
# else if( TYPE == "DIST" ){
# max_val = max(c(1,max(MAT)))
#}
dat = melt(MAT,na.rm = TRUE)
# class(dat); dim(dat); dat[1:5,]
gg = ggplot(data = dat,aes(x = Var1,y = Var2,fill = value)) +
geom_tile(color = "black") + ggtitle(LABEL) +
labs(fill = "Value")
if( is.null(LABx) ){
gg = gg + xlab("")
} else {
gg = gg + xlab(LABx)
}
if( is.null(LABy) ){
gg = gg + ylab("")
} else {
gg = gg + ylab(LABy)
}
if( max(dim(MAT)) <= 10 && DIGITS >= 0 )
gg = gg + geom_text(mapping = aes(label = smart_digits(value,DIGITS)))
if( TYPE %in% c("GSIM","DIST") ){
gg = gg + scale_fill_gradient2(midpoint = med_val,low = "deepskyblue",
mid = "white",high = "red",limit = c(min_val,max_val))
} else if( TYPE == "MUT" ){
gg = gg + scale_fill_gradient2(low = "black",
high = "red",limit = c(min_val,max_val))
}
gg = gg + guides(fill = guide_colorbar(frame.colour = "black")) +
scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
scale_y_discrete(guide = guide_axis(n.dodge = 2))
gg = gg + theme(legend.position = "right",
# legend.key.width = unit(1.5,'cm'),
legend.key.height = unit(1,'cm'),
legend.text = element_text(size = 12),
legend.title = element_text(size = 12,hjust = 0.5),
text = element_text(size = 12),
panel.background = element_blank(),
panel.grid.major = element_line(colour = "grey50",
size = 0.5,linetype = "dotted"),
axis.text = element_text(face = "italic"),
axis.text.x = element_text(angle = 0),
plot.title = element_text(hjust = 0.5))
return(gg)
}
hout = hclust(as.dist(1 - GS))
ord = hout$labels[hout$order]
show_tile(MAT = GS[ord,ord],
LABEL = "Simulated Gene Similarity",
TYPE = "GSIM",DIGITS = 2)
show_tile(MAT = 1 - GS[ord,ord],
LABEL = "Simulated Gene Dissimilarity",
TYPE = "GSIM",DIGITS = 2)
The function show_tile() is inherent to this vignette and not part of the ROKET package.
Simulate mutated gene statuses
# Mutated gene statuses
prob_mut = 0.2
prob_muts = c(1 - prob_mut,prob_mut)
while(TRUE){
ZZ = matrix(sample(c(0,1),NN*GG,replace = TRUE,prob = prob_muts),NN,GG)
# Ensure each sample has at least one mutated gene
if( min(rowSums(ZZ)) > 0 ) break
}
dimnames(ZZ) = list(NN_nms,GG_nms)
show_tile(MAT = ZZ,
LABEL = "Mutation Status: Gene by Sample",
TYPE = "MUT",DIGITS = 0)
# Store all distances
DD = array(data = NA,dim = c(NN,NN,5))
dimnames(DD)[1:2] = list(NN_nms,NN_nms)
dimnames(DD)[[3]] = c("EUC","OT_Balanced",sprintf("OT_LAM%s",c(0.5,1.0,5.0)))
We can look at the distribution of gene mutation frequencies.
freq = colSums(ZZ); # freq
dat = smart_df(GENE = names(freq),FREQ = as.integer(freq))
dat$GENE = factor(dat$GENE,levels = names(sort(freq,decreasing = TRUE)))
# dat
ggplot(data = dat,mapping = aes(x = GENE,y = FREQ)) +
geom_bar(stat = "identity") +
xlab("Gene") + ylab("Mutation Frequency") +
scale_x_discrete(guide = guide_axis(n.dodge = 2))
Euclidean distance
We can calculate the Euclidean distance, which does not incorporate relationships between pairs of genes.
DD[,,"EUC"] = as.matrix(dist(ZZ,diag = TRUE,upper = TRUE))
hout = hclust(as.dist(DD[,,"EUC"]))
ord = hout$labels[hout$order]
show_tile(MAT = DD[,,"EUC"][ord,ord],
LABEL = "Euclidean Pairwise Distances",
TYPE = "DIST",DIGITS = 2)
OT distance
Code between two samples
To demonstrate ROKET's functionality, the code below will run balanced OT (pre-normalizing input vectors) between two samples. Regardless of the values specified by LAMBDA1 and LAMBDA2 arguments, we need to set balance = TRUE. The OT cost matrix argument corresponds to 1 - GS, one minus the gene similarity matrix.
# Pick two samples
ii = 1
jj = 2
ZZ[c(ii,jj),colSums(ZZ[c(ii,jj),]) > 0]
outOT = run_myOT(XX = ZZ[ii,],YY = ZZ[jj,],
COST = 1 - GS,EPS = 1e-3,LAMBDA1 = 1,
LAMBDA2 = 1,balance = TRUE,verbose = FALSE)
# str(outOT)
# Optimal transport matrix
tmpOT = outOT$OT
tmpOT = tmpOT[rowSums(tmpOT) > 0,colSums(tmpOT) > 0]
show_tile(MAT = tmpOT,LABEL = "Balanced OT (showing only genes mutated in each sample)",
TYPE = "DIST",LABx = sprintf("Sample %s",ii),
LABy = sprintf("Sample %s",jj),
DIGITS = 2)
# Pairwise distance
outOT$DIST
Let's try again but with unbalanced OT and $\lambda = 0.5$. We need to set balance = FALSE and specify LAMBDA1 = 0.5 and LAMBDA2 = 0.5.
ZZ[c(ii,jj),colSums(ZZ[c(ii,jj),]) > 0]
LAM = 0.5
outOT = run_myOT(XX = ZZ[ii,],YY = ZZ[jj,],
COST = 1 - GS,EPS = 1e-3,LAMBDA1 = LAM,
LAMBDA2 = LAM,balance = FALSE,verbose = FALSE)
# str(outOT)
# Optimal transport matrix
tmpOT = outOT$OT
tmpOT = tmpOT[rowSums(tmpOT) > 0,colSums(tmpOT) > 0]
show_tile(MAT = tmpOT,
LABEL = "Unbalanced OT (showing only genes mutated in each sample)",TYPE = "DIST",
LABx = sprintf("Sample %s",ii),
LABy = sprintf("Sample %s",jj),
DIGITS = 2)
# Pairwise distance
outOT$DIST
Code between all sample pairs
ROKET can run optimal transport calculations across all $N$ choose 2 samples. Below is code to run balanced OT.
outOTs = run_myOTs(ZZ = t(ZZ),COST = 1 - GS,
EPS = 1e-3,balance = TRUE,LAMBDA1 = 1,
LAMBDA2 = 1,conv = 1e-5,max_iter = 3e3,
ncores = 1,verbose = FALSE)
hout = hclust(as.dist(outOTs))
ord = hout$labels[hout$order]
show_tile(MAT = outOTs[ord,ord],
LABEL = "Balanced OT Distances",
TYPE = "DIST",DIGITS = 1)
We can run calculations for various $\lambda$ values. For shorthand, $\lambda$ = Inf corresponds to balanced OT.
LAMs = c(0,0.5,1.0,5.0)
for(LAM in LAMs){
# LAM = LAMs[2]
BAL = ifelse(LAM == 0,TRUE,FALSE)
LAM2 = ifelse(BAL,1,LAM)
outOTs = run_myOTs(ZZ = t(ZZ),COST = 1 - GS,
EPS = 1e-3,balance = BAL,LAMBDA1 = LAM2,
LAMBDA2 = LAM2,conv = 1e-5,max_iter = 3e3,
ncores = 1,verbose = FALSE)
hout = hclust(as.dist(outOTs))
ord = hout$labels[hout$order]
LABEL = ifelse(BAL,"OT Distances (Balanced)",
sprintf("OT Distances (Lambda = %s)",LAM))
LABEL2 = ifelse(BAL,"OT_Balanced",sprintf("OT_LAM%s",LAM))
gg = show_tile(MAT = outOTs[ord,ord],
LABEL = LABEL,TYPE = "DIST",DIGITS = 2)
print(gg)
DD[,,LABEL2] = outOTs
rm(outOTs)
}
Dendrograms
We can see that Euclidean distance calculations on gene mutation statuses alone does not lead to strong evidence of sample clusters. However optimal transport-based distance calculations with integrated gene-gene similarities provide stronger evidence of sample clusters.
nms = dimnames(DD)[[3]]; nms
for(nm in nms){
# nm = nms[2]
hout = hclust(as.dist(DD[,,nm]))
hdend = as.dendrogram(hout)
dend_data = dendro_data(hdend,type = "rectangle")
gg = ggplot(dend_data$segments) +
geom_segment(aes(x = x,y = y,xend = xend,yend = yend)) +
ggtitle(nm) + xlab("") + ylab("") +
geom_text(data = dend_data$labels,
mapping = aes(x,y,label = label),vjust = 0.5,hjust = 1) +
theme_dendro() + coord_flip() +
theme(plot.title = element_text(hjust = 0.5))
print(gg)
rm(hout,hdend,dend_data,gg)
}
Kernel Regression and Association
The next step is to transform distance matrices into centered kernel matrices to perform hypothesis testing on our constructed kernels.
Models
For a binary or continuous outcome, $Y_i$, fitted with a generalized linear model, we have
$$
g\bigPar{E\bcSqu{Y_i}{\bf{X}_i,\bf{Z}_i}} = \bf{X}_i^\T \bf{\beta} + f(\bf{Z}_i) = \bf{X}_i^\T \bf{\beta} + \bf{K}_i^\T \bf{\alpha}
$$
and for Cox proportional hazards for survival outcomes, we have
$$
h(t;\bf{X}_i,\bf{Z}_i) = h_0(t) \nexp{\bf{X}_i^\T \bf{\beta} + f(\bf{Z}_i)} = h_0(t) \nexp{\bf{X}_i^\T \bf{\beta} + \bf{K}_i^\T \bf{\alpha}}
$$
where
- $h_0(t)$ is the baseline hazards function,
- $g\bigPar{\cdot}$ is a known link function,
- $f\bigPar{\cdot}$ is assumed to be generated by a positive semi-definite function $K\bigPar{\cdot,\cdot}$ such that $f(\cdot)$ lies in the reproducing kernel Hilbert space $\mathcal{H}_K$,
- $\bf{X}_i$ denotes $p$ baseline covariates to adjust for,
- $\bf{K}i = (K{i1},\ldots,K_{iN})^\T$, the $i$th column of kernel matrix $\bf{K}$, and
- $\bf{\alpha} = (\alpha_1,\ldots,\alpha_N)^\T$.
Distance to Kernel
The matrix $\bf{K}$ is constructed from the distance matrix (Euclidean or optimal transport), $\bf{D}$, by double centering the rows and columns of $\bf{D}$ with the formula
$$
\tilde{\bf{K}} = -\frac{1}{2} \bigPar{\bf{I}_N - \bf{J}_N \bf{J}_N^\T} \bf{D}^2 \bigPar{\bf{I}_N - \bf{J}_N \bf{J}_N^\T}
$$
where
- $\bf{D}^2$ is the element-wise squared distance matrix,
- $\bf{I}_N$ is an $N \times N$ identity matrix, and
- $\bf{J}_N$ is a column vector of $N$ ones.
Since $\tilde{\bf{K}}$ is not guaranteed to be positive semi-definite, we perform spectral decomposition and replace negative eigenvalues with zero and re-calculate the kernel to arrive at $\bf{K}$.
For a single candidate distance matrix, transform the distance matrix to a kernel matrix and convert the object into a list, denoted by the object KK, code is provided below from the R package, MiRKAT [@zhao2015testing;@plantinga2017mirkat].
# For example, with Euclidean distance
KK = MiRKAT::D2K(D = DD[,,"EUC"])
KK = list(EUC = KK)
If there are multiple distance matrices, in our case, contained in an array DD, store them into a list of kernel matrices KK. Example code is below.
KK = list()
for(nm in dimnames(DD)[[3]]){
KK[[nm]] = MiRKAT::D2K(D = DD[,,nm])
}
Hypothesis Testing
The user needs to pre-specify and fit the null model,
$$
H_0: \bf{\alpha} = 0,
$$
of baseline covariates, $\bf{X}_i$, to one of the three outcome models to obtain continuous, logistic, or martingale residuals (RESI). Some example code is provided below.
# Continuous
out_LM = lm(Y ~ .,data = data.frame(X))
RESI = residuals(out_LM)
# Logistic
out_LOG = glm(Y ~ .,data = data.frame(X),family = "binomial")
RESI = residuals(out_LOG)
# Survival
out_CX = coxph(Surv(TIME,CENS) ~ .,data = data.frame(X))
RESI = residuals(out_CX)
With one or multiple candidate kernels, ROKET will take
- a R array of kernel matrices (
aKK) (defined below),
- the residual vector (
RESI), and
- user-supplied number of permutations (
nPERMS)
to calculate individual kernel p-values as well as omnibus tests. An omnibus test assesses the significance of the minimum p-value kernel among kernels considered. The function kernTEST requires names(RESI) and dimension names per object as well as row/column names per kernel within aKK are named for sample order consistency. Sample code is provided below. Setting verbose = TRUE allows the user to track the permutation's progress, especially when requesting hundreds of thousands of permutations.
nPERMS = 1e5
nKK = length(KK)
# Array of kernels
aKK = array(data = NA,dim = dim(DD),
dimnames = dimnames(DD))
for(nm in dimnames(DD)[[3]]){
aKK[,,nm] = KK[[nm]]
}
# Create OMNI matrix
OMNI = matrix(0,nrow = 2,ncol = dim(aKK)[3],
dimnames = list(c("EUC","OT"),dimnames(aKK)[[3]]))
OMNI["EUC","EUC"] = 1
OMNI["OT",grepl("^OT",colnames(OMNI))] = 1
OMNI
# Hypothesis Testing
ROKET::kernTEST(RESI = RESI,
KK = aKK,
OMNI = OMNI,
nPERMS = nPERMS,
ncores = 1)
The final output contains individual kernel p-values and the omnibus p-value.
Have fun with utilizing kernel regression and optimal transport frameworks with ROKET!
Session Info
sessionInfo()
References
Try the ROKET package in your browser
Any scripts or data that you put into this service are public.
ROKET documentation built on April 4, 2025, 2:28 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE,comment = "#>", echo = TRUE,cache = TRUE, dev = "png") ot_eval = TRUE
\def\T{\text{T}} \newcommand{\bf}[1]{\mathbf{#1}} \newcommand{\bigPar}[1]{\left(#1\right)} \newcommand{\bigCur}[1]{\left{#1\right}} \newcommand{\bcSqu}[2]{\left[#1 \middle| #2\right]} \newcommand{\cE}[2]{E\bcSqu{#1}{#2}} \newcommand{\nexp}[1]{\exp\bigCur{#1}} \newcommand{\ind}[1]{1\bigCur{#1}} \newcommand{\w}[1]{\widehat{#1}}
Overview
Assuming all software dependencies and ROKET [@little2023associating] installation are installed, we can begin.
# Load libraries req_packs = c("devtools","smarter","ggplot2","reshape2", "survival","ggdendro","MiRKAT","ROKET") for(pack in req_packs){ library(package = pack,character.only = TRUE) } # List package's exported functions ls("package:ROKET") # Fix seed set.seed(2)
Distance Motivation
For $i = 1,\ldots,N$ and $g = 1,\ldots,G$, let
- $Z_{ig} = 1$ indicate the $i$th sample's $g$th gene is mutated and $Z_{ig} = 0$ otherwise
- $\bf{Z}i \equiv \bigPar{Z{i1},\ldots,Z_{iG}}^\T$
We would like to calculate the distance between the $i$th and $j$th samples in terms of mutated genes $\bf{Z}_i$ and $\bf{Z}_j$, denoted by $d\bigPar{\bf{Z}_i,\bf{Z}_j}$.
Optimal Transport
Background papers
Several optimal transport (OT) references to get familiarized with the framework, objective functions, updating equations, entropic regularizations, types of OT are provided:
- @cuturi2013sinkhorn,
- @zhang2021review,
- @jagarlapudi2020statistical,
- @fatras2021unbalanced,
- @chapel2021unbalanced
Basic Idea
For the $i$th and $j$th samples, let
- $\bf{P}^{(ij)} =$ an unknown $G$ by $G$ transport matrix to be estimated
- $P_{gh}^{(ij)} =$ the value of the $g$th row and $h$th column of $\bf{P}$, denotes the mass transported between the $g$th gene of the $i$th sample and the $h$th gene of the $j$th sample
- $\bf{W} =$ a known $G$ by $G$ cost matrix, independent of samples
- $W_{gh} =$ the value of the $g$th row and $h$th column of $\bf{W}$, denotes the cost to transport one unit of mass between the $g$th gene and $h$th gene
- The space of transport matrices $\bf{P}^{(ij)}$ to search over is governed through the penalized divergence between $\bf{Z}_i$ and the row sums of $\bf{P}^{(ij)}$ as well as between $\bf{Z}_j$ and the column sums of $\bf{P}^{(ij)}$
- $\w{\bf{P}^{(ij)}} = \text{argmin}{\bf{P}^{(ij)}} \bigPar{\sum_g \sum_h P{gh}^{(ij)} W_{gh}}$, the optimal transport matrix
- $d\bigPar{\bf{Z}i,\bf{Z}_j} = \sum_g \sum_h \w{P{gh}^{(ij)}} W_{gh}$
Balanced OT
If the mass of the two vectors are equal or normalized such that $\sum_g Z_{ig} = \sum_g Z_{jg}$, we could use balanced optimal transport.
Unbalanced OT
If the mass of the two vectors are not equal, $\sum_g Z_{ig} \neq \sum_g Z_{jg}$, we could use unbalanced optimal transport with penalty parameters.
Simulated Example
The code below will simulate samples and mutated genes. We welcome the reader to manipulate the following input arguments:
NNfor sample size,PPfor number of pathways,GGfor number of genes, andbnd_samethe lower bound gene similarity for genes sharing the same pathway and1 - bnd_same, an upper bound on gene similarity for genes not sharing the same pathway
Ideally, PP should be much less than GG to allow multiple genes to be allocated to each pathway.
# number of samples NN = 30 NN_nms = sprintf("S%s",seq(NN)) # number of pathways PP = 4 PP_nms = sprintf("P%s",seq(PP)) # number of genes GG = 30 GG_nms = sprintf("G%s",seq(GG)) # bound for gene similarity of two genes on same or different pathway bnd_same = 0.75 # Gene and pathway relationship GP = smart_df(PATH = sample(seq(PP),GG,replace = TRUE), GENE = seq(GG)) table(GP$PATH) # gene-gene similarity matrix GS = matrix(NA,GG,GG) dimnames(GS) = list(GG_nms,GG_nms) diag(GS) = 1 tmp_mat = t(combn(seq(GG),2)) for(ii in seq(nrow(tmp_mat))){ G1 = tmp_mat[ii,1] G2 = tmp_mat[ii,2] same = GP$PATH[GP$GENE == G1] == GP$PATH[GP$GENE == G2] if( same ) GS[G1,G2] = runif(1,bnd_same,1) else GS[G1,G2] = runif(1,0,1 - bnd_same) } GS[lower.tri(GS)] = t(GS)[lower.tri(GS)] # round(GS,3)
Gene-Gene Similarity
Let's take a look at the gene similarity matrix.
show_tile = function(MAT,LABEL,TYPE = NULL, LABx = NULL,LABy = NULL,DIGITS = 1){ min_val = min(MAT) max_val = max(MAT) med_val = (min_val + max_val) / 2 if( is.null(TYPE) ) stop("Specify TYPE") if( isSymmetric(MAT) ) MAT[upper.tri(MAT,diag = !TRUE)] = NA if( TYPE == "GSIM" ){ max_val = 1 med_val = 0.5 } # else if( TYPE == "DIST" ){ # max_val = max(c(1,max(MAT))) #} dat = melt(MAT,na.rm = TRUE) # class(dat); dim(dat); dat[1:5,] gg = ggplot(data = dat,aes(x = Var1,y = Var2,fill = value)) + geom_tile(color = "black") + ggtitle(LABEL) + labs(fill = "Value") if( is.null(LABx) ){ gg = gg + xlab("") } else { gg = gg + xlab(LABx) } if( is.null(LABy) ){ gg = gg + ylab("") } else { gg = gg + ylab(LABy) } if( max(dim(MAT)) <= 10 && DIGITS >= 0 ) gg = gg + geom_text(mapping = aes(label = smart_digits(value,DIGITS))) if( TYPE %in% c("GSIM","DIST") ){ gg = gg + scale_fill_gradient2(midpoint = med_val,low = "deepskyblue", mid = "white",high = "red",limit = c(min_val,max_val)) } else if( TYPE == "MUT" ){ gg = gg + scale_fill_gradient2(low = "black", high = "red",limit = c(min_val,max_val)) } gg = gg + guides(fill = guide_colorbar(frame.colour = "black")) + scale_x_discrete(guide = guide_axis(n.dodge = 2)) + scale_y_discrete(guide = guide_axis(n.dodge = 2)) gg = gg + theme(legend.position = "right", # legend.key.width = unit(1.5,'cm'), legend.key.height = unit(1,'cm'), legend.text = element_text(size = 12), legend.title = element_text(size = 12,hjust = 0.5), text = element_text(size = 12), panel.background = element_blank(), panel.grid.major = element_line(colour = "grey50", size = 0.5,linetype = "dotted"), axis.text = element_text(face = "italic"), axis.text.x = element_text(angle = 0), plot.title = element_text(hjust = 0.5)) return(gg) } hout = hclust(as.dist(1 - GS)) ord = hout$labels[hout$order] show_tile(MAT = GS[ord,ord], LABEL = "Simulated Gene Similarity", TYPE = "GSIM",DIGITS = 2) show_tile(MAT = 1 - GS[ord,ord], LABEL = "Simulated Gene Dissimilarity", TYPE = "GSIM",DIGITS = 2)
The function show_tile() is inherent to this vignette and not part of the ROKET package.
Simulate mutated gene statuses
# Mutated gene statuses prob_mut = 0.2 prob_muts = c(1 - prob_mut,prob_mut) while(TRUE){ ZZ = matrix(sample(c(0,1),NN*GG,replace = TRUE,prob = prob_muts),NN,GG) # Ensure each sample has at least one mutated gene if( min(rowSums(ZZ)) > 0 ) break } dimnames(ZZ) = list(NN_nms,GG_nms) show_tile(MAT = ZZ, LABEL = "Mutation Status: Gene by Sample", TYPE = "MUT",DIGITS = 0) # Store all distances DD = array(data = NA,dim = c(NN,NN,5)) dimnames(DD)[1:2] = list(NN_nms,NN_nms) dimnames(DD)[[3]] = c("EUC","OT_Balanced",sprintf("OT_LAM%s",c(0.5,1.0,5.0)))
We can look at the distribution of gene mutation frequencies.
freq = colSums(ZZ); # freq dat = smart_df(GENE = names(freq),FREQ = as.integer(freq)) dat$GENE = factor(dat$GENE,levels = names(sort(freq,decreasing = TRUE))) # dat ggplot(data = dat,mapping = aes(x = GENE,y = FREQ)) + geom_bar(stat = "identity") + xlab("Gene") + ylab("Mutation Frequency") + scale_x_discrete(guide = guide_axis(n.dodge = 2))
Euclidean distance
We can calculate the Euclidean distance, which does not incorporate relationships between pairs of genes.
DD[,,"EUC"] = as.matrix(dist(ZZ,diag = TRUE,upper = TRUE)) hout = hclust(as.dist(DD[,,"EUC"])) ord = hout$labels[hout$order] show_tile(MAT = DD[,,"EUC"][ord,ord], LABEL = "Euclidean Pairwise Distances", TYPE = "DIST",DIGITS = 2)
OT distance
Code between two samples
To demonstrate ROKET's functionality, the code below will run balanced OT (pre-normalizing input vectors) between two samples. Regardless of the values specified by LAMBDA1 and LAMBDA2 arguments, we need to set balance = TRUE. The OT cost matrix argument corresponds to 1 - GS, one minus the gene similarity matrix.
# Pick two samples ii = 1 jj = 2 ZZ[c(ii,jj),colSums(ZZ[c(ii,jj),]) > 0] outOT = run_myOT(XX = ZZ[ii,],YY = ZZ[jj,], COST = 1 - GS,EPS = 1e-3,LAMBDA1 = 1, LAMBDA2 = 1,balance = TRUE,verbose = FALSE) # str(outOT) # Optimal transport matrix tmpOT = outOT$OT tmpOT = tmpOT[rowSums(tmpOT) > 0,colSums(tmpOT) > 0] show_tile(MAT = tmpOT,LABEL = "Balanced OT (showing only genes mutated in each sample)", TYPE = "DIST",LABx = sprintf("Sample %s",ii), LABy = sprintf("Sample %s",jj), DIGITS = 2) # Pairwise distance outOT$DIST
Let's try again but with unbalanced OT and $\lambda = 0.5$. We need to set balance = FALSE and specify LAMBDA1 = 0.5 and LAMBDA2 = 0.5.
ZZ[c(ii,jj),colSums(ZZ[c(ii,jj),]) > 0] LAM = 0.5 outOT = run_myOT(XX = ZZ[ii,],YY = ZZ[jj,], COST = 1 - GS,EPS = 1e-3,LAMBDA1 = LAM, LAMBDA2 = LAM,balance = FALSE,verbose = FALSE) # str(outOT) # Optimal transport matrix tmpOT = outOT$OT tmpOT = tmpOT[rowSums(tmpOT) > 0,colSums(tmpOT) > 0] show_tile(MAT = tmpOT, LABEL = "Unbalanced OT (showing only genes mutated in each sample)",TYPE = "DIST", LABx = sprintf("Sample %s",ii), LABy = sprintf("Sample %s",jj), DIGITS = 2) # Pairwise distance outOT$DIST
Code between all sample pairs
ROKET can run optimal transport calculations across all $N$ choose 2 samples. Below is code to run balanced OT.
outOTs = run_myOTs(ZZ = t(ZZ),COST = 1 - GS, EPS = 1e-3,balance = TRUE,LAMBDA1 = 1, LAMBDA2 = 1,conv = 1e-5,max_iter = 3e3, ncores = 1,verbose = FALSE) hout = hclust(as.dist(outOTs)) ord = hout$labels[hout$order] show_tile(MAT = outOTs[ord,ord], LABEL = "Balanced OT Distances", TYPE = "DIST",DIGITS = 1)
We can run calculations for various $\lambda$ values. For shorthand, $\lambda$ = Inf corresponds to balanced OT.
LAMs = c(0,0.5,1.0,5.0) for(LAM in LAMs){ # LAM = LAMs[2] BAL = ifelse(LAM == 0,TRUE,FALSE) LAM2 = ifelse(BAL,1,LAM) outOTs = run_myOTs(ZZ = t(ZZ),COST = 1 - GS, EPS = 1e-3,balance = BAL,LAMBDA1 = LAM2, LAMBDA2 = LAM2,conv = 1e-5,max_iter = 3e3, ncores = 1,verbose = FALSE) hout = hclust(as.dist(outOTs)) ord = hout$labels[hout$order] LABEL = ifelse(BAL,"OT Distances (Balanced)", sprintf("OT Distances (Lambda = %s)",LAM)) LABEL2 = ifelse(BAL,"OT_Balanced",sprintf("OT_LAM%s",LAM)) gg = show_tile(MAT = outOTs[ord,ord], LABEL = LABEL,TYPE = "DIST",DIGITS = 2) print(gg) DD[,,LABEL2] = outOTs rm(outOTs) }
Dendrograms
We can see that Euclidean distance calculations on gene mutation statuses alone does not lead to strong evidence of sample clusters. However optimal transport-based distance calculations with integrated gene-gene similarities provide stronger evidence of sample clusters.
nms = dimnames(DD)[[3]]; nms for(nm in nms){ # nm = nms[2] hout = hclust(as.dist(DD[,,nm])) hdend = as.dendrogram(hout) dend_data = dendro_data(hdend,type = "rectangle") gg = ggplot(dend_data$segments) + geom_segment(aes(x = x,y = y,xend = xend,yend = yend)) + ggtitle(nm) + xlab("") + ylab("") + geom_text(data = dend_data$labels, mapping = aes(x,y,label = label),vjust = 0.5,hjust = 1) + theme_dendro() + coord_flip() + theme(plot.title = element_text(hjust = 0.5)) print(gg) rm(hout,hdend,dend_data,gg) }
Kernel Regression and Association
The next step is to transform distance matrices into centered kernel matrices to perform hypothesis testing on our constructed kernels.
Models
For a binary or continuous outcome, $Y_i$, fitted with a generalized linear model, we have
$$ g\bigPar{E\bcSqu{Y_i}{\bf{X}_i,\bf{Z}_i}} = \bf{X}_i^\T \bf{\beta} + f(\bf{Z}_i) = \bf{X}_i^\T \bf{\beta} + \bf{K}_i^\T \bf{\alpha} $$
and for Cox proportional hazards for survival outcomes, we have
$$ h(t;\bf{X}_i,\bf{Z}_i) = h_0(t) \nexp{\bf{X}_i^\T \bf{\beta} + f(\bf{Z}_i)} = h_0(t) \nexp{\bf{X}_i^\T \bf{\beta} + \bf{K}_i^\T \bf{\alpha}} $$
where
- $h_0(t)$ is the baseline hazards function,
- $g\bigPar{\cdot}$ is a known link function,
- $f\bigPar{\cdot}$ is assumed to be generated by a positive semi-definite function $K\bigPar{\cdot,\cdot}$ such that $f(\cdot)$ lies in the reproducing kernel Hilbert space $\mathcal{H}_K$,
- $\bf{X}_i$ denotes $p$ baseline covariates to adjust for,
- $\bf{K}i = (K{i1},\ldots,K_{iN})^\T$, the $i$th column of kernel matrix $\bf{K}$, and
- $\bf{\alpha} = (\alpha_1,\ldots,\alpha_N)^\T$.
Distance to Kernel
The matrix $\bf{K}$ is constructed from the distance matrix (Euclidean or optimal transport), $\bf{D}$, by double centering the rows and columns of $\bf{D}$ with the formula
$$ \tilde{\bf{K}} = -\frac{1}{2} \bigPar{\bf{I}_N - \bf{J}_N \bf{J}_N^\T} \bf{D}^2 \bigPar{\bf{I}_N - \bf{J}_N \bf{J}_N^\T} $$
where
- $\bf{D}^2$ is the element-wise squared distance matrix,
- $\bf{I}_N$ is an $N \times N$ identity matrix, and
- $\bf{J}_N$ is a column vector of $N$ ones.
Since $\tilde{\bf{K}}$ is not guaranteed to be positive semi-definite, we perform spectral decomposition and replace negative eigenvalues with zero and re-calculate the kernel to arrive at $\bf{K}$.
For a single candidate distance matrix, transform the distance matrix to a kernel matrix and convert the object into a list, denoted by the object KK, code is provided below from the R package, MiRKAT [@zhao2015testing;@plantinga2017mirkat].
# For example, with Euclidean distance KK = MiRKAT::D2K(D = DD[,,"EUC"]) KK = list(EUC = KK)
If there are multiple distance matrices, in our case, contained in an array DD, store them into a list of kernel matrices KK. Example code is below.
KK = list() for(nm in dimnames(DD)[[3]]){ KK[[nm]] = MiRKAT::D2K(D = DD[,,nm]) }
Hypothesis Testing
The user needs to pre-specify and fit the null model,
$$ H_0: \bf{\alpha} = 0, $$
of baseline covariates, $\bf{X}_i$, to one of the three outcome models to obtain continuous, logistic, or martingale residuals (RESI). Some example code is provided below.
# Continuous out_LM = lm(Y ~ .,data = data.frame(X)) RESI = residuals(out_LM) # Logistic out_LOG = glm(Y ~ .,data = data.frame(X),family = "binomial") RESI = residuals(out_LOG) # Survival out_CX = coxph(Surv(TIME,CENS) ~ .,data = data.frame(X)) RESI = residuals(out_CX)
With one or multiple candidate kernels, ROKET will take
- a R array of kernel matrices (
aKK) (defined below), - the residual vector (
RESI), and - user-supplied number of permutations (
nPERMS)
to calculate individual kernel p-values as well as omnibus tests. An omnibus test assesses the significance of the minimum p-value kernel among kernels considered. The function kernTEST requires names(RESI) and dimension names per object as well as row/column names per kernel within aKK are named for sample order consistency. Sample code is provided below. Setting verbose = TRUE allows the user to track the permutation's progress, especially when requesting hundreds of thousands of permutations.
nPERMS = 1e5 nKK = length(KK) # Array of kernels aKK = array(data = NA,dim = dim(DD), dimnames = dimnames(DD)) for(nm in dimnames(DD)[[3]]){ aKK[,,nm] = KK[[nm]] } # Create OMNI matrix OMNI = matrix(0,nrow = 2,ncol = dim(aKK)[3], dimnames = list(c("EUC","OT"),dimnames(aKK)[[3]])) OMNI["EUC","EUC"] = 1 OMNI["OT",grepl("^OT",colnames(OMNI))] = 1 OMNI # Hypothesis Testing ROKET::kernTEST(RESI = RESI, KK = aKK, OMNI = OMNI, nPERMS = nPERMS, ncores = 1)
The final output contains individual kernel p-values and the omnibus p-value.
Have fun with utilizing kernel regression and optimal transport frameworks with ROKET!
Session Info
sessionInfo()
References
Try the ROKET package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.