Contents

1 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.

2 Install MICSQTL

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("MICSQTL")

3 Load packages

library(MICSQTL)

Load packages for making plots.

library(reshape2)
library(GGally)
library(ggplot2)

4 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:

This example data can be loaded by calling:

data(se)

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

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.

4.1 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, MuSiC, 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.

4.2 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.

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)

Boxplots of estimated cell proportions using cross-source.

4.3 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.

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.

4.3.1 Comparison to 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"
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.

Top3 principal components of transcriptome.

4.4 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.

head(rowData(se))
#> DataFrame with 6 rows and 4 columns
#>               Chr     Start       End      Symbol
#>       <character> <integer> <integer> <character>
#> AAGAB          15  67202823  67254631       AAGAB
#> AARS2           6  44300549  44313323       AARS2
#> AASS            7 122076491 122133726        AASS
#> ABAT           16   8735739   8781427        ABAT
#> ABCA1           9 104784317 104903679       ABCA1
#> ABCA2           9 137007931 137028140       ABCA2

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.

head(metadata(se)$anno_SNP)
#>        CHROM       POS          ID
#> 332373     9 137179658 9:137179658
#> 237392     9 104596634 9:104596634
#> 106390     9  28487163  9:28487163
#> 304108     9 126307371 9:126307371
#> 295846     9 122787821 9:122787821
#> 126055     9  33975396  9:33975396

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.

target_protein <- rowData(se)[rowData(se)$Chr == 9, ][seq_len(3), "Symbol"]
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.

unlist(lapply(metadata(se)$choose_SNP_list, length))
#>   ABCA1   ABCA2 AGTPBP1 
#>      26      22       7

4.5 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.

system.time(se <- csQTL(se))

We can check the results from csQTL analysis for one of target proteins:

res <- metadata(se)$TOAST_output[[2]]
head(res[order(apply(res, 1, min)), ])

5 Licenses of the analysis methods

method citation
TCA 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 Feng, Qing, et al. “Angle-based joint and individual variation explained.” Journal of multivariate analysis 166 (2018): 241-265.
TOAST 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

#> R version 4.4.0 (2024-04-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB              LC_COLLATE=C              
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] MICSQTL_1.2.2               SummarizedExperiment_1.34.0
#>  [3] Biobase_2.64.0              GenomicRanges_1.56.0       
#>  [5] GenomeInfoDb_1.40.0         IRanges_2.38.0             
#>  [7] S4Vectors_0.42.0            BiocGenerics_0.50.0        
#>  [9] MatrixGenerics_1.16.0       matrixStats_1.3.0          
#> [11] BiocStyle_2.32.0           
#> 
#> loaded via a namespace (and not attached):
#>   [1] pbapply_1.7-2           formatR_1.14            rlang_1.1.3            
#>   [4] magrittr_2.0.3          ggridges_0.5.6          e1071_1.7-14           
#>   [7] compiler_4.4.0          gdata_3.0.0             vctrs_0.6.5            
#>  [10] quadprog_1.5-8          stringr_1.5.1           pkgconfig_2.0.3        
#>  [13] crayon_1.5.2            fastmap_1.1.1           backports_1.4.1        
#>  [16] XVector_0.44.0          utf8_1.2.4              rmarkdown_2.26         
#>  [19] pracma_2.4.4            nloptr_2.0.3            UCSC.utils_1.0.0       
#>  [22] purrr_1.0.2             xfun_0.43               zlibbioc_1.50.0        
#>  [25] cachem_1.0.8            jsonlite_1.8.8          EpiDISH_2.20.0         
#>  [28] DelayedArray_0.30.1     BiocParallel_1.38.0     gmodels_2.19.1         
#>  [31] broom_1.0.5             parallel_4.4.0          R6_2.5.1               
#>  [34] bslib_0.7.0             stringi_1.8.4           RColorBrewer_1.1-3     
#>  [37] limma_3.60.0            GGally_2.2.1            car_3.1-2              
#>  [40] jquerylib_0.1.4         Rcpp_1.0.12             bookdown_0.39          
#>  [43] iterators_1.0.14        knitr_1.46              Matrix_1.7-0           
#>  [46] nnls_1.5                splines_4.4.0           tidyselect_1.2.1       
#>  [49] abind_1.4-5             yaml_2.3.8              doParallel_1.0.17      
#>  [52] codetools_0.2-20        lattice_0.22-6          tibble_3.2.1           
#>  [55] plyr_1.8.9              evaluate_0.23           lambda.r_1.2.4         
#>  [58] futile.logger_1.4.3     ggstats_0.6.0           proxy_0.4-27           
#>  [61] ggpubr_0.6.0            pillar_1.9.0            BiocManager_1.30.23    
#>  [64] carData_3.0-5           foreach_1.5.2           generics_0.1.3         
#>  [67] ggplot2_3.5.1           munsell_0.5.1           TOAST_1.18.0           
#>  [70] scales_1.3.0            gtools_3.9.5            class_7.3-22           
#>  [73] glue_1.7.0              tools_4.4.0             data.table_1.15.4      
#>  [76] ggsignif_0.6.4          grid_4.4.0              dirmult_0.1.3-5        
#>  [79] tidyr_1.3.1             matrixcalc_1.0-6        colorspace_2.1-0       
#>  [82] GenomeInfoDbData_1.2.12 rsvd_1.0.5              cli_3.6.2              
#>  [85] futile.options_1.0.1    fansi_1.0.6             S4Arrays_1.4.0         
#>  [88] dplyr_1.1.4             corpcor_1.6.10          gtable_0.3.5           
#>  [91] rstatix_0.7.2           sass_0.4.9              digest_0.6.35          
#>  [94] SparseArray_1.4.3       htmltools_0.5.8.1       lifecycle_1.0.4        
#>  [97] httr_1.4.7              TCA_1.2.1               locfdr_1.1-8           
#> [100] statmod_1.5.0           MASS_7.3-60.2