---
title: "A Vignette for DeMixT"
author: "Zeya Wang, Fan Gao, Shaolong Cao and Peng Yang"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{DeMixT.Rmd}
%\VignetteEngine{knitr::rmarkdown}
\usepackage[utf8]{inputenc}
---
```{r setup, include=FALSE, message=FALSE}
library(ggplot2)
library(DeMixT)
# knitr::opts_chunk$set(echo = TRUE)
```
# 1. Introduction
Transcriptomic deconvolution in cancer and other heterogeneous tissues remains
challenging. Available methods lack the ability to estimate both
component-specific proportions and expression profiles for individual samples.
We develop a three-component deconvolution model, DeMixT, for expression data
from a mixture of cancerous tissues, infiltrating immune cells and tumor
microenvironment. DeMixT is a software package that performs deconvolution on
transcriptome data from a mixture of two or three components.
DeMixT is a frequentist-based method and fast in yielding accurate estimates of
cell proportions and compart-ment-specific expression profiles for
two-component \and three-component deconvolution problem. Our method promises
to provide deeper insight into cancer biomarkers and assist in the development
of novel prognostic markers and therapeutic strategies.
The function DeMixT is designed to finish the whole pipeline of deconvolution
for two or three components. The newly added DeMixT_GS function is designed to
estimates the proportions of mixed samples for each mixing component based on a
new approach to select genes more effectively that utilizes profile likelihood.
DeMixT_S1 function is designed to estimate the proportions of all mixed samples
for each mixing component based on the gene differential expressions to select
genes. DeMixT_S2 function is designed to estimate the component-specific
deconvolved expressions of individual mixed samples for a given set of genes.
# 2 Feature Description
The DeMixT R-package builds the transcriptomic deconvolution with a couple of
novel features into R-based standard analysis pipeline through Bioconductor.
DeMixT showed high accuracy and efficiency from our designed experiment.
Hence, DeMixT can be considered as an important step towards linking tumor
transcriptomic data with clinical outcomes.
Different from most previous computational deconvolution methods, DeMixT has
integrated new features for the deconvolution with more than 2 components.
**Joint estimation**: jointly estimate component proportions and expression
profiles for individual samples by requiring reference samples instead of
reference genes; For the three-component deconvolution considering immune
infiltration, it provides a comprehensive view of tumor-stroma-immune
transcriptional dynamics, as compared to methods that address only immune
subtypes within the immune component, in each tumor sample.
**Efficient estimation**: DeMixT adopts an approach of iterated conditional
modes (ICM) to guarantee a rapid convergence to a local maximum. We also design
a novel gene-set-based component merging approach to reduce the bias of
proportion estimation for three-component deconvolutionthe.
**Parallel computing**: OpenMP enables parallel computing on single computer by
taking advantage of the multiple cores shipped on modern CPUs. The ICM
framework further enables parallel computing, which helps compensate for the
expensive computing time used in the repeated numerical double integrations.
# 3. Installation
## 3.1 Source file
DeMixT source files are compatible with Windows, Linux and macOS.
DeMixT_1.2.5 is the latest version, which is for a computer that has OpenMP.
It can be downloaded from Bioconductor website.
https://bioconductor.org/packages/release/bioc/html/DeMixT.html
To install DeMixT_1.2.5 from GitHub, start R and enter:
```{r install_DeMixT}
# devtools::install_github("wwylab/DeMixT")
```
For more information, please visit:
<http://bioinformatics.mdanderson.org/main/DeMixT>
## 3.2 Functions
The following table shows the functions included in DeMixT.
Table Header | Second Header
------- | ----------------------------------
DeMixT | Deconvolution of tumor samples with two or three components
DeMixT_GS | Estimates the proportions of mixed samples for each mixing component based on a new approach to select genes that utilizes profile likelihood
DeMixT_S1 | Estimates the proportions of mixed samples for each mixing component
DeMixT_S2 |Deconvolves expressions of each sample for unknown component
Optimum_KernelC | Call the C function used for parameter estimation in DeMixT
# 4. Methods
## 4.1 Model
Let \(Y_{ig}\) be the observed expression levels of the raw measured data from
clinically derived malignant tumor samples for gene \(g, g = 1, \cdots, G\) and
sample \(i, i = 1, \cdots, My\). \(G\) denotes the total number of probes/genes
and \(My\) denotes the number of samples. The observed expression levels for
solid tumors can be modeled as a linear combination of raw expression levels
from three components:
\[ {Y_{ig}} = \pi _{1,i}N_{1,ig} + \pi _{2,i}N_{2,ig} +
(1 - \pi_{1,i} - \pi _{2,i}){T_{ig}} \label{eq:1} \]
Here \(N_{1,ig}\), \(N_{2,ig}\) and \({T_{ig}}\) are the unobserved raw
expression levels from each of the three components. We call the two components
for which we require reference samples the \(N_1\)-component and the
\(N_2\)-component. We call the unknown component the T-component. We let
\(\pi_{1,i}\) denote the proportion of the \(N_1\)-component, \(\pi_{2,i}\)
denote the proportion of the \(N_2\)-component, and \(1 - \pi_{1,i}-\pi_{2,i}\)
denote the proportion of the T-component. We assume that the mixing proportions
of one specific sample remain the same across all genes.
Our model allows for one component to be unknown, and therefore does not
require reference profiles from all components. A set of samples for
\(N_{1,ig}\) and \(N_{2,ig}\), respectively, needs to be provided as input
data. This three-component deconvolution model is applicable to the linear
combination of any three components in any type of material. It can also be
simplified to a two-component model, assuming there is just one
\(N\)-component. For application in this paper, we consider tumor (\(T\)),
stromal (\(N_1\)) and immune components (\(N_2\)) in an admixed sample (\(Y\)).
Following the convention that \(\log_2\)-transformed microarray gene expression
data follow a normal distribution, we assume that the raw measures \(N_{1,ig}
\sim LN({\mu _{{N_1}g}},\sigma _{{N_1}g}^2)\), \(N_{2,ig}
\sim LN({\mu _{{N_2}g}},\sigma _{{N_2}g}^2)\) and \({T_{ig}}
\sim LN({\mu _{Tg}}, \sigma _{Tg}^2)\), where LN denotes a \(\log_2\)-normal
distribution and \(\sigma _{{N_1}g}^2\),\(\sigma _{{N_2}g}^2\),
\(\sigma _{Tg}^2\) reflect the variations under \(\log_2\)-transformed data.
Consequently, our model can be expressed as the convolution of the density
function for three \(\log_2\)-normal distributions. Because there is no closed
form of this convolution, we use numerical integration to evaluate the complete
likelihood function (see the full likelihood in the Supplementary Materials).
## 4.2 The DeMixT algorithm for deconvolution
DeMixT estimates all distribution parameters and cellular proportions and
reconstitutes the expression profiles for all three components for each gene
and each sample. The estimation procedure (summarized in Figure 1b) has two
main steps as follows.
1. Obtain a set of parameters \(\{\pi_{1,i}, \pi_{2,i}\}_{i=1}^{My}\), \(\{\mu_T,
\sigma_T\}_{g=1}^G\) to maximize the complete likelihood function, for which
\(\{\mu_{N_{1,g}}, \sigma_{N_{1,g}}, \mu_{N_{2,g}},
\sigma_{N_{2,g}}\}_{g=1}^G\) were already estimated from the available
unmatched samples of the \(N_1\) and \(N_2\) component tissues.
(See further details in our paper.)
2. Reconstitute the expression profiles by searching each set of \(\{n_{1,ig},
n_{2,ig}\}\) that maximizes the joint density of \(N_{1,ig}\), \(N_{2,ig}\) and
\(T_{ig}\). The value of \(t_{ig}\) is solved as \({y_{ig}} -
{{\hat \pi }_{1,i}}{n_{1,ig}} - {{\hat \pi }_{2,i}}{n_{2,ig}}\).
These two steps can be separately implemented using the function DeMixT_S1 or
DeMixT_GS for the first step and DeMixT_S2 for the second, which are combined
in the function DeMixT(Note: DeMixT_GS is the default function for first step).
In version 1.2.5, DeMixT added simulated normal reference samples, i.e., spike-in,
based on the observed normal reference samples. It has been shown to improve
accuracy in proportion estimation for the scenario where a dataset consists of
samples where true tumor proportions are skewed to the high end.
```{r Algorithm, echo=FALSE, out.width='100%'}
knitr::include_graphics(path = paste0("Algorithm.png"))
```
# 5. Examples
## 5.1 Simulated two-component data
```{r, sim_2comp_GS, results="hide", message=FALSE}
data("test.data.2comp")
# res.GS = DeMixT_GS(data.Y = test.data.2comp$data.Y,
# data.N1 = test.data.2comp$data.N1,
# niter = 30, nbin = 50, nspikein = 50,
# if.filter = TRUE, ngene.Profile.selected = 150,
# mean.diff.in.CM = 0.25, ngene.selected.for.pi = 150,
# tol = 10^(-5))
load('Res_2comp/res.GS.RData')
```
```{r sim_2comp_GS_res}
head(t(res.GS$pi))
head(res.GS$gene.name)
```
```{r, sim_2comp_S2, results="hide", message=FALSE}
data("test.data.2comp")
# res.S2 <- DeMixT_S2(data.Y = test.data.2comp$data.Y,
# data.N1 = test.data.2comp$data.N1,
# data.N2 = NULL,
# givenpi = c(t(res.S1$pi[-nrow(res.GS$pi),])), nbin = 50)
load('Res_2comp/res.S2.RData')
```
```{r sim_2comp_S2_res}
head(res.S2$decovExprT[,1:5],3)
head(res.S2$decovExprN1[,1:5],3)
head(res.S2$decovMu,3)
head(res.S2$decovSigma,3)
```
## 5.2 Simulated two-component data with DeMixT Spike-in Normal
In the simulation,
````markdown
## Simulate MuN and MuT for each gene
MuN <- rnorm(G, 7, 1.5)
MuT <- rnorm(G, 7, 1.5)
Mu <- cbind(MuN, MuT)
## Simulate SigmaN and SigmaT for each gene
SigmaN <- runif(n = G, min = 0.1, max = 0.8)
SigmaT <- runif(n = G, min = 0.1, max = 0.8)
## Simulate Tumor Proportion
PiT = truncdist::rtrunc(n = My,
spec = 'norm',
mean = 0.55,
sd = 0.2,
a = 0.25,
b = 0.95)
## Simulate Data
for(k in 1:G){
data.N1[k,] <- 2^rnorm(M1, MuN[k], SigmaN[k]); # normal reference
True.data.T[k,] <- 2^rnorm(My, MuT[k], SigmaT[k]); # True Tumor
True.data.N1[k,] <- 2^rnorm(My, MuN[k], SigmaN[k]); # True Normal
data.Y[k,] <- pi[1,]*True.data.N1[k,] + pi[2,]*True.data.T[k,] # Mixture Tumor
}
````
where $\pi_i \in (0.25, 0.95)$ is from truncated normal distribution. In general, the true distribution of tumor proportion does not follow a uniform distribution between $[0,1]$, but instead skewed to the upper part of the interval.
```{r read_data_2comp, warning=FALSE, message=FALSE}
# ## DeMixT_S1 without Spike-in Normal
# res.S1 = DeMixT_S1(data.Y = test.data.2comp$data.Y,
# data.N1 = test.data.2comp$data.N1,
# niter = 30, nbin = 50, nspikein = 0,
# if.filter = TRUE,
# mean.diff.in.CM = 0.25, ngene.selected.for.pi = 150,
# tol = 10^(-5))
# ## DeMixT_S1 with Spike-in Normal
# res.S1.SP = DeMixT_S1(data.Y = test.data.2comp$data.Y,
# data.N1 = test.data.2comp$data.N1,
# niter = 30, nbin = 50, nspikein = 50,
# if.filter = TRUE,
# mean.diff.in.CM = 0.25, ngene.selected.for.pi = 150,
# tol = 10^(-5))
# ## DeMixT_GS with Spike-in Normal
# res.GS.SP = DeMixT_GS(data.Y = test.data.2comp$data.Y,
# data.N1 = test.data.2comp$data.N1,
# niter = 30, nbin = 50, nspikein = 50,
# if.filter = TRUE, ngene.Profile.selected = 150,
# mean.diff.in.CM = 0.25, ngene.selected.for.pi = 150,
# tol = 10^(-5))
load('Res_2comp/res.S1.RData'); load('Res_2comp/res.S1.SP.RData');
load('Res_2comp/res.GS.RData'); load('Res_2comp/res.GS.SP.RData');
```
This simulation was designed to compare previous DeMixT resutls with DeMixT spike-in results under both gene selection method.
```{r sim_2comp, fig.height = 4, fig.width = 6, fig.align='center', warning=FALSE}
res.2comp = as.data.frame(cbind(round(rep(t(test.data.2comp$pi[2,]),3),2),
round(c(t(res.S1$pi[2,]),t(res.S1.SP$pi[2,]), t(res.GS.SP$pi[2,])),2),
rep(c('S1','S1-SP','GS-SP'), each = 100)), num = 1:2)
res.2comp$V1 <- as.numeric(as.character(res.2comp$V1))
res.2comp$V2 <- as.numeric(as.character(res.2comp$V2))
res.2comp$V3 = as.factor(res.2comp$V3)
names(res.2comp) = c('True.Proportion', 'Estimated.Proportion', 'Method')
## Plot
ggplot(res.2comp, aes(x=True.Proportion, y=Estimated.Proportion, group = Method, color=Method, shape=Method)) +
geom_point() +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "black", lwd = 0.5) +
xlim(0,1) + ylim(0,1) +
scale_shape_manual(values=c(seq(1:3))) +
labs(x = 'True Proportion', y = 'Estimated Proportion')
```
## 5.3 Simulated three-component data
In this simulation,
````markdown
G <- G1 + G2
## Simulate MuN1, MuN2 and MuT for each gene
MuN1 <- rnorm(G, 7, 1.5)
MuN2_1st <- MuN1[1:G1] + truncdist::rtrunc(n = 1,
spec = 'norm',
mean = 0,
sd = 1.5,
a = -0.1,
b = 0.1)
MuN2_2nd <- c()
for(l in (G1+1):G){
tmp <- MuN1[l] + truncdist::rtrunc(n = 1,
spec = 'norm',
mean = 0,
sd = 1.5,
a = 0.1,
b = 3)^rbinom(1, size=1, prob=0.5)
while(tmp <= 0) tmp <- MuN1[l] + truncdist::rtrunc(n = 1,
spec = 'norm',
mean = 0,
sd = 1.5,
a = 0.1,
b = 3)^rbinom(1, size=1, prob=0.5)
MuN2_2nd <- c(MuN2_2nd, tmp)
}
## Simulate SigmaN1, SigmaN2 and SigmaT for each gene
SigmaN1 <- runif(n = G, min = 0.1, max = 0.8)
SigmaN2 <- runif(n = G, min = 0.1, max = 0.8)
SigmaT <- runif(n = G, min = 0.1, max = 0.8)
## Simulate Tumor Proportion
pi <- matrix(0, 3, My)
pi[1,] <- runif(n = My, min = 0.01, max = 0.97)
for(j in 1:My){
pi[2, j] <- runif(n = 1, min = 0.01, max = 0.98 - pi[1,j])
pi[3, j] <- 1 - sum(pi[,j])
}
## Simulate Data
for(k in 1:G){
data.N1[k,] <- 2^rnorm(M1, MuN1[k], SigmaN1[k]); # normal reference 1
data.N2[k,] <- 2^rnorm(M2, MuN2[k], SigmaN2[k]); # normal reference 1
True.data.T[k,] <- 2^rnorm(My, MuT[k], SigmaT[k]); # True Tumor
True.data.N1[k,] <- 2^rnorm(My, MuN1[k], SigmaN1[k]); # True Normal 1
True.data.N2[k,] <- 2^rnorm(My, MuN2[k], SigmaN2[k]); # True Normal 1
data.Y[k,] <- pi[1,]*True.data.N1[k,] + pi[2,]*True.data.N2[k,] +
pi[3,]*True.data.T[k,] # Mixture Tumor
}
````
where $G1$ is the number of genes that $\mu_{N1}$ is close to $\mu_{N2}$.
```{r read_data_3comp, warning=FALSE, message=FALSE}
data("test.data.3comp")
# res.S1 <- DeMixT_S1(data.Y = test.data.3comp$data.Y, data.N1 = test.data.3comp$data.N1,
# data.N2 = test.data.3comp$data.N2, if.filter = TRUE)
load('Res_3comp/res.S1.RData');
```
```{r sim_3comp, fig.height = 4, fig.width = 6, fig.align='center', warning=FALSE}
res.3comp= as.data.frame(cbind(round(t(matrix(t(test.data.3comp$pi), nrow = 1)),2),
round(t(matrix(t(res.S1$pi), nrow = 1)),2),
rep(c('N1','N2','T'), each = 20)))
res.3comp$V1 <- as.numeric(as.character(res.3comp$V1))
res.3comp$V2 <- as.numeric(as.character(res.3comp$V2))
res.3comp$V3 = as.factor(res.3comp$V3)
names(res.3comp) = c('True.Proportion', 'Estimated.Proportion', 'Component')
## Plot
ggplot(res.3comp, aes(x=True.Proportion, y=Estimated.Proportion, group = Component, color=Component, shape=Component)) +
geom_point() +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "black", lwd = 0.5) +
xlim(0,1) + ylim(0,1) +
scale_shape_manual(values=c(seq(1:3))) +
labs(x = 'True Proportion', y = 'Estimated Proportion')
```
## 5.4 PRAD in TCGA Dataset
```{r read_data_PRAD, warning=FALSE, message=FALSE}
load('res.PRAD.RData');
```
```{r PRAD, fig.height = 4, fig.width = 6, fig.align='center', warning=FALSE}
res.PRAD.df = as.data.frame(cbind(res.PRAD$res.GS.PRAD, res.PRAD$res.GS.SP.PRAD))
res.PRAD.df$V1 <- as.numeric(as.character(res.PRAD.df$V1))
res.PRAD.df$V2 <- as.numeric(as.character(res.PRAD.df$V2))
names(res.PRAD.df) = c('Estimated.Proportion', 'Estimated.Proportion.SP')
## Plot
ggplot(res.PRAD.df, aes(x=Estimated.Proportion, y=Estimated.Proportion.SP)) +
geom_point() +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "black", lwd = 0.5) +
xlim(0,1) + ylim(0,1) +
scale_shape_manual(values=c(seq(1:3))) +
labs(x = 'Estimated Proportion', y = 'Estimated Proportion After Spike-in')
```
Tumor proportion estimation remains consistence with and without spike-in normal, when its distribution is roughly symmetric around 0.5.
# 6. Session Info
```{r}
sessionInfo(package = "DeMixT")
```