---
title: "MICSQTL: Multi-omic deconvolution, Integration and Cell-type-specific Quantitative Trait Loci"
author:
- name: Yue Pan
  affiliation: Department of Biostatistics, St. Jude Children's Research Hospital
  email: Yue.Pan@STJUDE.ORG
- name: Qian Li
  affiliation: Department of Biostatistics, St. Jude Children's Research Hospital
package: MICSQTL
output: BiocStyle::html_document
vignette: >
  %\VignetteIndexEntry{MICSQTL: Multi-omic deconvolution, Integration and Cell-type-specific Quantitative Trait Loci}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
abstract: >
  MICSQTL (Multi-omic deconvolution, Integration and Cell-type-specific Quantitative Trait Loci) is a tool that estimates cell type proportions in bulk proteomes by leveraging shared information in matched transcriptomes. Based on the deconvoluted cell fractions, MICSQTL further integrates and visualizes multi-source profiles at bulk or cell type levels, as well as performs Quantitative Trait Loci (QTL) mapping at cellular resolution. This is a novel pipeline for multi-omic integrative data analysis that fulfills the need for comprehensive analysis of diverse data sources. This pipeline enables valuable insights into cellular composition, facilitates cell-type-specific protein QTL mapping, and streamlines multi-modal data integration and dimension reduction.
---
    
```{r set, echo = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)
```

# Introduction

Our pipeline, `MICSQTL`, integrates RNA and protein expressions to detect potential cell marker proteins and estimate cell abundance in mixed proteomes without a reference signature matrix. `MICSQTL` enables cell-type-specific quantitative trait loci (QTL) mapping for proteins or transcripts using bulk expression data and estimated cellular composition per molecule type, eliminating the necessity for single-cell sequencing. We use matched transcriptome-proteome from human brain frontal cortex tissue samples to demonstrate the input and output of our tool.

![MICSQTL workflow.](fig1.PNG)


# Install MICSQTL

```{r setup, message = FALSE, warning = FALSE, eval = FALSE}
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("MICSQTL")
```

# Load packages

```{r lib, message = FALSE, warning = FALSE}
library(MICSQTL)
```

Load packages for making plots.

```{r lib2, message = FALSE, warning = FALSE, eval = FALSE}
library(reshape2)
library(GGally)
library(ggplot2)
```

# Quick start

To conduct the analysis, the initiation involves the use of a 
`SummarizedExperiment` object that contains bulk protein expression data, which 
can be rescaled using either log or MinMax transformations, in the assays slot. 
The row metadata (rowData slot) should contain information about the protein 
features. Additionally, incorporation of bulk gene expression data, which can 
also be rescaled using either log or MinMax transformations and needs to be 
consistent with the bulk protein expression data, as well as a reference file 
(depending on the chosen method), as integral elements within the metadata slot, 
are imperative. For more accurate cell-type fraction estimations, it's 
recommended to include only marker genes. 

For easier illustration purposes, we provide an example `SummarizedExperiment` 
object within this package, which contains the following elements:
    
-   protein_data (assay): An example proteomics data (on log scale) 
with 2,242 rows (protein) and 127 columns (sample).

-   anno_protein (rowData): A data frame with 2,242 rows and 4 columns 
(Chr, Start, End, Symbol) as annotations of each protein from `protein_data`.

-   ref_protein (in metadata): A signature matrix with 2,242 rows (protein) and 
4 columns (cell types), which serves as a reference of known cellular signatures 
(on log scale).

-   gene_data (in metadata): A data frame with 2,867 rows (genes) and 127 
columns (sample) (on log scale).

-   ref_gene (in metadata): A signature matrix with 4,872 rows (genes) and 5 
columns (cell types), which serves as a reference of known cellular signatures 
(on log scale).

-   prop_gene (in metadata): A pre-defined deconvoluted transcriptome proportion 
matrix.

-   SNP_data (in metadata): A sparse matrix with 2,000 rows (SNP), which stores 
the information of genetic variants at each location from one chromosome and 127 
columns (sample, should match the sample in `protein_data`). Each matrix entry 
corresponds to the genotype group indicator (0, 1 or 2) for a sample at a 
genetic location.

-   anno_SNP (in metadata): A data frame with 2,000 rows and 3 columns 
(CHROM, POS, ID), which stores Annotations of each SNP from `SNP_data`.

-   meta (in metadata):A data frame with 127 rows (sample) and 2 columns 
(disease status and gender) as metadata.

-   cell_counts (in metadata): A matrix containing cell counts across multiple 
subjects, where subjects are represented as rows and cell types as columns. 
Each entry (i, j) in the matrix indicates the count of cells belonging to the ith 
subject and jth cell type.

This example data can be loaded by calling:

```{r obj}
data(se)
```

Below is an example code for building the `SummarizedExperiment` object from raw 
data frames or matrices. 

```{r eg, eval = FALSE}
se <- SummarizedExperiment(
    assays = list(protein = your_protein_data),
    rowData = your_anno_protein
)
metadata(se) <- list(
    gene_data = your_gene_data
)
```

Additional metadata can be incorporated using a command such as 
`metadata(se)$new_data <- new_data` if further information is necessary for 
visualization or csQTL (cell-type specific quantitative trait loci) analysis. 
For detailed instructions, please refer to the following sections and the 
function documentation.




## Cell-type proportion deconvolution

This step estimates the proportions of cell types for each molecule type.

In this current version, only `nnls` is supported as a single-source deconvolution method. Users can utilize other methods such as [CIBERSORT](https://cibersortx.stanford.edu/cshome.php), [MuSiC](https://www.nature.com/articles/s41467-018-08023-x), etc., to obtain the proportion estimates. These estimates will be useful as initial values in the subsequent deconvolution based on cross-source. It is important to note that the samples used in the cell-type proportion estimates must match the samples in the bulk protein expression data.

## Cross-source cell-type proportion deconvolution

The reference matrix for pure cell proteomics may be incomplete due to limitations inherent in single-cell proteomics technologies. To address this issue, we propose a novel cross-source cell-type fraction deconvolution method "Joint Non-negative Matrix Factorization" (JNMF) that capitalizes on matched bulk transcriptome-proteome data. In the following example, we illustrate the process of estimating protein proportions by integrating information from deconvoluted transcriptomes.

There are multiple options available for initializing cellular fractions and purified proteomics on a sample-wise basis, each with different input requirements. The following example illustrates using the CIBERSORT method to estimate initial proportions (`pinit`), coupled with an external reference matrix (`ref_pnl`) containing gene expression profiles or molecular signatures of different cell types. These references are typically obtained from small-scale single-cell or flow cytometry experiments. Please note that `ref_pnl` should have the same rescaling transformation as the bulk transcriptomes/proteomics and 
have be non-negative.

It is recommended to use the `ajive_decomp` function (further discussed in the following section) with the `refactor_loading = TRUE` to enhance joint deconvolution. This setting enables cross-source feature selection aimed at identifying potential protein cell markers.

If the estimated cell-type proportions contain any invalid or non-finite values, consider adjusting the step size.

```{r cross, message = FALSE, warning = FALSE, results = FALSE, eval = FALSE}
se <- ajive_decomp(se, use_marker = FALSE, refactor_loading = TRUE)
se <- deconv(se, source = "cross", method = "JNMF", 
             Step = c(10^(-9), 10^(-7)),
             use_refactor = 1000,
             pinit = se@metadata$prop_gene,
             ref_pnl = se@metadata$ref_gene)
```

```{r plot2, echo = FALSE, warning = FALSE, eval = FALSE}
ggplot(
    cbind(data.frame(melt(metadata(se)$prop), metadata(se)$meta)),
    aes(x = Var2, y = value, fill = Var2)
) +
    geom_point(
        position = position_jitterdodge(
            jitter.width = 0.1,
            dodge.width = 0.7
        ),
        aes(fill = Var2, color = Var2),
        pch = 21, alpha = 0.5
    ) +
    geom_boxplot(lwd = 0.7, outlier.shape = NA) +
    theme_classic() +
    facet_wrap(~disease) +
    xlab("Cell type") +
    ylab("Estimated proportion") +
    theme(legend.position = "none")
```

![Boxplots of estimated cell proportions using cross-source.](plot2-1.png)


## Integrative visualization

AJIVE (Angle based Joint and Individual Variation Explained) is useful when there are multiple data matrices measured on the same set of samples. It decomposes each data matrix as three parts: (1) Joint variation across data types (2) Individual structured variation for each data type and (3) Residual noise.

It is similar as principal component analysis (PCA), but principal component analysis only takes a single data set and decomposes it into modes of variation that maximize variation. AJIVE finds joint modes of variation from multiple data sources.

Common normalized scores are one of the desirable output to explore the joint behavior that is shared by different data sources. Below we show the visualization of common normalized scores. It is clear that the disease status of these samples are well separated by the first common normalized scores.

```{r ajive, eval = FALSE}
se <- ajive_decomp(se, plot = TRUE,
                   group_var = "disease",
                   scatter = TRUE, scatter_x = "cns_1", scatter_y = "cns_2")
metadata(se)$cns_plot
```

![Sample group separation based on common normalized scores from AJIVE.](ajive-1.png)

### Comparison to PCA

```{r pca}
pca_res <- prcomp(t(assay(se)), rank. = 3, scale. = FALSE)
pca_res_protein <- data.frame(pca_res[["x"]])
pca_res_protein <- cbind(pca_res_protein, metadata(se)$meta$disease)
colnames(pca_res_protein)[4] <- "disease"
```

```{r pcaplot, eval = FALSE}
ggpairs(pca_res_protein,
        columns = seq_len(3), aes(color = disease, alpha = 0.5),
        upper = list(continuous = "points")
) + theme_classic()


pca_res <- prcomp(t(metadata(se)$gene_data), rank. = 3, scale. = FALSE)
pca_res_gene <- data.frame(pca_res[["x"]])
pca_res_gene <- cbind(pca_res_gene, metadata(se)$meta$disease)
colnames(pca_res_gene)[4] <- "disease"
ggpairs(pca_res_gene,
        columns = seq_len(3), aes(color = disease, alpha = 0.5),
        upper = list(continuous = "points")
) + theme_classic()
```

![Top3 principal components of proteome.](pcaplot-1.png)

![Top3 principal components of transcriptome.](pcaplot-2.png)

## Feature filtering

The feature filtering can be applied at both proteins/genes and SNPs. This step is optional but highly recommended to filter out some features that are not very informative or do not make much sense biologically. Note that this function is required to run even no filtering is expected to be done (just set `filter_method = "null"`) to obtain a consistent object format for downstream analysis.

To apply feature filtering, annotation files for protein/gene and SNPs are required. The annotation file for proteins/genes should be stored in `rowData()`, where each row corresponds to a protein/gene with it's symbol as row names. The first column should be a character vector indicating which chromosome each protein or gene is on. In addition, it should contain at least a "Start" column with numeric values indicating the start position on that chromosome, a "End" column with numeric values indicating the end position on that chromosome and a "Symbol" column as a unique name for each protein or gene.

```{r filter1}
head(rowData(se))
```

The information from genetic variants should be stored in a P (the number of SNP) by N (the number of samples, should match the sample in `counts` slot) matrix contained as an element (`SNP_data`) in `metadata` slot. Each matrix entry corresponds to the genotype group indicator (0 for 0/0, 1 for 0/1 and 2 for 1/1) for a sample at a genetic location. The annotations of these SNP should be stored as an element (`anno_SNP`) in `metadata` slot. It should include at least the following columns: (1) "CHROM" (which chromosome the SNP is on); (2) "POS" (position of that SNP) and (3) "ID" (a unique identifier for each SNP, usually a combination of chromosome and its position).

The example SNP data provided here were restricted to chromosome 9 only. In practice, the SNPs may from multiple or even all chromosomes.

```{r filter2}
head(metadata(se)$anno_SNP)
```

For filtering at protein or gene level, only those symbols contained in `target_SNP` argument will be kept and if not provided, all SNPs will be used for further filtering.

For filtering at SNP level, there are three options: (1) filter out the SNPs that have minor allele frequency below the threshold defined by `filter_allele` argument (`filter_method = "allele"`); (2) filter out the SNPs that the fraction of samples in the smallest genotype group below the threshold defined by `filter_geno` argument (`filter_method = "allele"`) and (3) restrict to cis-regulatory variants (`filter_method = "distance"`): the SNPs up to 1 Mb proximal to the start of the gene. Both filtering methods can be applied simultaneously by setting `filter_method = c("allele", "distance")`.

To simplify the analysis, we only test 3 targeted proteins from chromosome 9 as an example.

```{r filter3}
target_protein <- rowData(se)[rowData(se)$Chr == 9, ][seq_len(3), "Symbol"]
```

```{r filter 4}
se <- feature_filter(se,
    target_protein = target_protein,
    filter_method = c("allele", "distance"),
    filter_allele = 0.15,
    filter_geno = 0.05,
    ref_position = "TSS"
)
```

The results after filtering will be stored as an element (`choose_SNP_list`) in `metadata` slot. It is a list with the length of the number of proteins for downstream analysis. Each element stores the index of SNPs to be tested for corresponding protein. The proteins with no SNPs correspond to it will be removed from the returned list.

In this example, the number of SNPs corresponding to each protein after filtering ranges from 7 to 26.

```{r filter5}
unlist(lapply(metadata(se)$choose_SNP_list, length))
```

## csQTL analysis

In this step, the `TOAST` method is implemented for cell-type-specific differential expression analysis based on samples' genotype.

The result will be stored as an element (`TOAST_output`) in `metadata` slot. It is a list with the same length as tested proteins or genes where each element consists of a table including protein or gene symbol, SNP ID and p-values from each cell type. A significant p-value indicates that the protein or gene expression is different among the sample from different genotype groups.

```{r csQTL1, eval = FALSE}
system.time(se <- csQTL(se))
```

We can check the results from csQTL analysis for one of target proteins:
    
```{r csQTL2, eval = FALSE}
res <- metadata(se)$TOAST_output[[2]]
head(res[order(apply(res, 1, min)), ])
```


# Licenses of the analysis methods

| method                                                                 | citation                                                                                                                                                              |
|------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [TCA](https://cran.r-project.org/web/packages/TCA/index.html)          | Rahmani, Elior, et al. "Cell-type-specific resolution epigenetics without the need for cell sorting or single-cell biology." Nature communications 10.1 (2019): 3417. |
| [AJIVE](https://github.com/idc9/r_jive)                                | Feng, Qing, et al. "Angle-based joint and individual variation explained." Journal of multivariate analysis 166 (2018): 241-265.                                      |
| [TOAST](http://bioconductor.org/packages/release/bioc/html/TOAST.html) | Li, Ziyi, and Hao Wu. "TOAST: improving reference-free cell composition estimation by cross-cell type differential analysis." Genome biology 20.1 (2019): 1-17.       |

# Session info {.unnumbered}

```{r sessionInfo, echo=FALSE}
sessionInfo()
```
