# MarkerPen: Marker Gene Detection via Penalized Principal Component Analysis

#### Yixuan Qiu, Jiebiao Wang, Jing Lei, and Kathryn Roeder

#### 2021-03-14

## Introduction

The `markerpen` R package uses penalized principal component analysis to detect cell-type-specific marker genes from bulk tissue RNA sequencing data. It is a semi-supervised algorithm that uses bulk transcriptome data to refine a prior marker gene list.

## Installation

`markerpen` can be installed directly from CRAN:

    install.packages("markerpen")

A C++ compiler that supports the C++11 standard is needed to install `markerpen` from source.

For best performance, it is strongly suggested linking your R to the OpenBLAS library for matrix computation, although this step is optional. You can achieve this with the help of the ropenblas package.

## Example

In this section we show a real example of selecting marker genes from bulk tissue RNA sequencing data.

First, read in the data file located in the `examples` folder:

    library(markerpen)
    library(scales)

    # A subset of the ROSMAP data
    # Subsampled from the geneExprRaw.txt file located in
    # https://github.com/ellispatrick/CortexCellDeconv/tree/master/CellTypeDeconvAnalysis/Data
    load(system.file("examples", "gene_expr.RData", package = "markerpen"))

    # Convert gene name to Ensembl name
    ind = match(rownames(dat), markerpen::gene_mapping$name)
    ind = na.omit(ind)
    ensembl = markerpen::gene_mapping$ensembl[ind]

    # Get expression matrix - rows are observations, columns are genes
    mat_exp = t(dat[markerpen::gene_mapping$name[ind], ])
    colnames(mat_exp) = ensembl
    print(mat_exp[1:5, 1:5])

    ##           ENSG00000175899 ENSG00000128274 ENSG00000158122 ENSG00000103591
    ## X50107871        53.48904       12.058472        5.969030        32.64904
    ## X50300408        33.57736       11.149030        4.301594        35.66727
    ## X17929065       123.91786        9.063111        4.769154        33.99096
    ## X50105437        33.62506       10.395273        3.439445        45.69874
    ## X50107583       107.42628       22.870649        7.812406        38.85704
    ##           ENSG00000115977
    ## X50107871        328.9962
    ## X50300408        179.5002
    ## X17929065        326.8588
    ## X50105437        254.7298
    ## X50107583        252.5395

Next, read in the prior marker gene lists, collected from the published literature. We also restrict the search range for each cell type to a subset of the whole genome.

The two R data files are provided in the `examples` folder.

    # Read in prior marker genes
    load(system.file("examples", "published_markers.RData", package = "markerpen"))
    load(system.file("examples", "markers_range.RData", package = "markerpen"))

Below is the main part of the analysis: selecting marker genes for major cell types. After obtaining the marker genes for one cell type, we remove those genes from the search range for the next cell type, in order to make markers for different cell types non-overlapping. Theoretically the order of the cell types in computing has an impact on the final result, but in practice the impact is small if the cell types can be well separated.

    # Markers for astrocytes
    ast_re = refine_markers(mat_exp, markers_range$astrocytes, pub_markers$astrocytes,
                            lambda = 0.45, w = 1.5, maxit = 500, eps = 1e-3, verbose = 0)
    # Remove selected markers from the expression matrix
    mat_rest = mat_exp[, setdiff(colnames(mat_exp), ast_re$markers)]

    # Markers for oligodendrocytes
    oli_re = refine_markers(mat_rest, markers_range$oligodendrocytes, pub_markers$oligodendrocytes,
                            lambda = 0.45, w = 1.5, maxit = 500, eps = 1e-3, verbose = 0)
    mat_rest = mat_rest[, setdiff(colnames(mat_rest), oli_re$markers)]

    # Markers for microglia
    mic_re = refine_markers(mat_rest, markers_range$microglia, pub_markers$microglia,
                            lambda = 0.45, w = 1.5, maxit = 500, eps = 1e-3, verbose = 0)
    mat_rest = mat_rest[, setdiff(colnames(mat_rest), mic_re$markers)]

    # Markers for endothelial
    end_re = refine_markers(mat_rest, markers_range$endothelial, pub_markers$endothelial,
                            lambda = 0.45, w = 1.5, maxit = 500, eps = 1e-3, verbose = 0)
    mat_rest = mat_rest[, setdiff(colnames(mat_rest), end_re$markers)]

    # Markers for neurons
    neu_re = refine_markers(mat_rest, markers_range$neurons, pub_markers$neurons,
                            lambda = 0.45, w = 1.5, maxit = 500, eps = 1e-3, verbose = 0)

    # Refined markers
    markers_re = list(astrocytes       = ast_re$markers,
                      oligodendrocytes = oli_re$markers,
                      microglia        = mic_re$markers,
                      endothelial      = end_re$markers,
                      neurons          = neu_re$markers)

We post-process the selected genes by ordering them and selecting the top 50 markers for each cell type.

    # Post-process selected markers
    # Pick the first 50 ordered markers
    cor_markers = cor(mat_exp[, unlist(markers_re)])
    markers_ord = sort_markers(cor_markers, markers_re)
    markers_ord = lapply(markers_ord, head, n = 50)

Finally, visualize the sample correlation matrices on published markers and refined markers.

    # Function to visualize the sample correlation matrix
    vis_cor = function(mat_exp, markers)
    {
        all_genes = colnames(mat_exp)
        markers = intersect(unlist(markers), all_genes)
        cor_markers = cor(mat_exp[, unlist(markers)])
        p = nrow(cor_markers)

        cols = c("#08306b", "#08519c", "#2171b5", "#6baed6", "#9ecae1", "#c6dbef", "#deebf7",
                 "#ffffff",
                 "#fcf1f1", "#fae1e1", "#facdcd", "#f49c9c", "#f56566", "#f13a3c", "#d00003")

        ncols = length(cols)
        cols = scales::gradient_n_pal(cols, values = (0:ncols) / ncols)((1:100) / 100)
        op = par(mar = c(0, 0, 0, 0))
        image(cor_markers[, p:1], col = cols, breaks = (-50:50) / 50, asp = 1, axes = FALSE)
        par(op)
    }

Published markers:

    vis_cor(mat_exp, pub_markers)

Refined markers:

    vis_cor(mat_exp, markers_re)

Ordered top 50 markers for each cell type:

    vis_cor(mat_exp, markers_ord)