# CiteFuse: getting started

Yingxin Lin1 and Hani Jieun Kim1,2,3

1School of Mathematics and Statistics, The University of Sydney, Australia
2Charles Perkins Centre, The University of Sydney, Australia
3Computational Systems Biology Group, Children’s Medical Research Institute, Faculty of Medicine and Health, The University of Sydney, Australia

#### 13 April 2021

# 1 Introduction

`CiteFuse` is a computational framework that implements a suite of methods and tools for CITE-seq data from pre-processing through to integrative analytics. This includes doublet detection, network-based modality integration, cell type clustering, differential RNA and protein expression analysis, ADT evaluation, ligand-receptor interaction analysis, and interactive web-based visualisation of the analyses. This vignette demonstrates the usage of `CiteFuse` on a subset data of CITE-seq data from human PBMCs as an example (Mimitou et al., 2019).

First, install `CiteFuse` using `BiocManager`.

    if (!requireNamespace("BiocManager", quietly = TRUE)) {
     install.packages("BiocManager")
    }
    BiocManager::install("CiteFuse")

    library(CiteFuse)
    library(scater)
    library(SingleCellExperiment)
    library(DT)

    data("CITEseq_example", package = "CiteFuse")
    names(CITEseq_example)
    #> [1] "RNA" "ADT" "HTO"
    lapply(CITEseq_example, dim)
    #> $RNA
    #> [1] 19521   500
    #>
    #> $ADT
    #> [1]  49 500
    #>
    #> $HTO
    #> [1]   4 500

Here, we start from a list of three matrices of unique molecular identifier (UMI), antibody derived tags (ADT) and hashtag oligonucleotide (HTO) count, which have common cell names. There are 500 cells in our subsetted dataset. And characteristically of CITE-seq data, the matrices are matched, meaning that for any given cell we know the expression level of their RNA transcripts (genome-wide) and its corresponding cell surface protein expression. The `preprocessing` function will utilise the three matrices and its common cell names to create a `SingleCellExperiment` object, which stores RNA data in an `assay` and `ADT` and `HTO` data within in the `altExp` slot.

    sce_citeseq <- preprocessing(CITEseq_example)
    sce_citeseq
    #> class: SingleCellExperiment
    #> dim: 19521 500
    #> metadata(0):
    #> assays(1): counts
    #> rownames(19521): hg19_AL627309.1 hg19_AL669831.5 ... hg19_MT-ND6
    #>   hg19_MT-CYB
    #> rowData names(0):
    #> colnames(500): AAGCCGCGTTGTCTTT GATCGCGGTTATCGGT ... TTGGCAACACTAGTAC
    #>   GCTGCGAGTTGTGGCC
    #> colData names(0):
    #> reducedDimNames(0):
    #> altExpNames(2): ADT HTO

# 2 Detecting both cross- and within-sample doublets using `CiteFuse`

## 2.1 HTO Normalisation and Visualisation

The function `normaliseExprs` is used to scale the alternative expression. Here, we used it to perform log-transformation of the `HTO` count, by setting `transform = "log"`.

    sce_citeseq <- normaliseExprs(sce = sce_citeseq,
                                  altExp_name = "HTO",
                                  transform = "log")

Then we can perform dimension reduction on the `HTO` count by using `runTSNE` or `runUMAP`, then use `visualiseDim` function to visualise the reduced dimension plot. Our CITE-seq dataset contain data from four samples that were pooled before sequencing. The samples were multiplexed through cell hashing (Stoekius et al., 2018). The four clusters observed on reduced dimension plots equate to the four different samples.

    sce_citeseq <- scater::runTSNE(sce_citeseq,
                                   altexp = "HTO",
                                   name = "TSNE_HTO",
                                   pca = TRUE)

    visualiseDim(sce_citeseq,
                 dimNames = "TSNE_HTO") + labs(title = "tSNE (HTO)")

    sce_citeseq <- scater::runUMAP(sce_citeseq,
                                   altexp = "HTO",
                                   name = "UMAP_HTO")

    visualiseDim(sce_citeseq,
                 dimNames = "UMAP_HTO") + labs(title = "UMAP (HTO)")

## 2.2 Doublet identification step 1: cross-sample doublet detection

An important step in single cell data analysis is the removal of doublets. Doublets form as a result of co-encapsulation of cells within a droplet, leading to a hybrid transcriptome from two or more cells. In CiteFuse, we implement a step-wise doublet detection approach to remove doublets. We first identify the cross-sample doublets via the `crossSampleDoublets` function.

    sce_citeseq <- crossSampleDoublets(sce_citeseq)
    #> number of iterations= 20
    #> number of iterations= 24
    #> number of iterations= 46
    #> number of iterations= 50

The results of the cross sample doublets are then saved in `colData` as `doubletClassify_between_label` and `doubletClassify_between_class`.

    table(sce_citeseq$doubletClassify_between_label)
    #>
    #>                 1                 2                 3                 4
    #>               115               121                92               129
    #> doublet/multiplet
    #>                43
    table(sce_citeseq$doubletClassify_between_class)
    #>
    #>           Singlet doublet/multiplet
    #>               457                43

We can then highlight the cross-sample doublets in our tSNE plot of HTO count.

    visualiseDim(sce_citeseq,
                 dimNames = "TSNE_HTO",
                 colour_by = "doubletClassify_between_label")

Furthermore, `plotHTO` function allows us to plot the pairwise scatter HTO count. Any cells that show co-expression of orthologocal HTOs (red) are considered as doublets.

    plotHTO(sce_citeseq, 1:4)

## 2.3 Doublet identification step 1: within-sample doublet detection

We then identify the within-sample doublets via the `withinSampleDoublets` function.

    sce_citeseq <- withinSampleDoublets(sce_citeseq,
                                        minPts = 10)

The results of the cross sample doublets are then saved in the `colData` as `doubletClassify_within_label` and `doubletClassify_within_class`.

    table(sce_citeseq$doubletClassify_within_label)
    #>
    #>  Doublets(Within)_1  Doublets(Within)_2  Doublets(Within)_3  Doublets(Within)_4
    #>                   3                   7                   4                   6
    #> NotDoublets(Within)
    #>                 480
    table(sce_citeseq$doubletClassify_within_class)
    #>
    #> Doublet Singlet
    #>      20     480

Again, we can visualise the within-sample doublets in our tSNE plot.

    visualiseDim(sce_citeseq,
                 dimNames = "TSNE_HTO",
                 colour_by = "doubletClassify_within_label")

Finally, we can filter out the doublet cells (both within and between batches) for the downstream analysis.

    sce_citeseq <- sce_citeseq[, sce_citeseq$doubletClassify_within_class == "Singlet" & sce_citeseq$doubletClassify_between_class == "Singlet"]
    sce_citeseq
    #> class: SingleCellExperiment
    #> dim: 19521 437
    #> metadata(3): doubletClassify_between_threshold
    #>   doubletClassify_between_resultsMat doubletClassify_within_resultsMat
    #> assays(1): counts
    #> rownames(19521): hg19_AL627309.1 hg19_AL669831.5 ... hg19_MT-ND6
    #>   hg19_MT-CYB
    #> rowData names(0):
    #> colnames(437): GATCGCGGTTATCGGT GGCTGGTAGAGGTTAT ... TTGGCAACACTAGTAC
    #>   GCTGCGAGTTGTGGCC
    #> colData names(5): doubletClassify_between_label
    #>   doubletClassify_between_class nUMI doubletClassify_within_label
    #>   doubletClassify_within_class
    #> reducedDimNames(2): TSNE_HTO UMAP_HTO
    #> altExpNames(2): ADT HTO

# 3 Clustering

## 3.1 Performing SNF

The first step of analysis is to integrate the RNA and ADT matrix. We use a popular integration algorithm called similarity network fusion (SNF) to integrate the multiomic data.

    sce_citeseq <- scater::logNormCounts(sce_citeseq)
    sce_citeseq <- normaliseExprs(sce_citeseq, altExp_name = "ADT", transform = "log")
    system.time(sce_citeseq <- CiteFuse(sce_citeseq))
    #> Calculating affinity matrix
    #> Performing SNF
    #>    user  system elapsed
    #>    4.76    0.19    4.95

We now proceed with the fused matrix, which is stored as `SNF_W` in our `sce_citeseq` object.

## 3.2 Performing spectral clustering

CiteFuse implements two different clustering algorithms on the fused matrix, spectral clustering and Louvain clustering. First, we perform spectral clustering with sufficient numbers of `K` and use the eigen values to determine the optimal number of clusters.

    SNF_W_clust <- spectralClustering(metadata(sce_citeseq)[["SNF_W"]], K = 20)
    #> Computing Spectral Clustering
    plot(SNF_W_clust$eigen_values)

    which.max(abs(diff(SNF_W_clust$eigen_values)))
    #> [1] 5

Using the optimal cluster number defined from the previous step, we can now use the `spectralClutering` function to cluster the single cells by specifying the number of clusters in `K`. The function takes a cell-to-cell similarity matrix as an input. We have already created the fused similarity matrix from `CiteFuse`. Since the `CiteFuse` function creates and stores the similarity matries from ADT and RNA expression, as well the fused matrix, we can use these two to compare the clustering outcomes by data modality.

    SNF_W_clust <- spectralClustering(metadata(sce_citeseq)[["SNF_W"]], K = 5)
    #> Computing Spectral Clustering
    sce_citeseq$SNF_W_clust <- as.factor(SNF_W_clust$labels)

    SNF_W1_clust <- spectralClustering(metadata(sce_citeseq)[["ADT_W"]], K = 5)
    #> Computing Spectral Clustering
    sce_citeseq$ADT_clust <- as.factor(SNF_W1_clust$labels)

    SNF_W2_clust <- spectralClustering(metadata(sce_citeseq)[["RNA_W"]], K = 5)
    #> Computing Spectral Clustering
    sce_citeseq$RNA_clust <- as.factor(SNF_W2_clust$labels)

## 3.3 Visualisation

The outcome of the clustering can be easily visualised on a reduced dimensions plot by highlighting the points by cluster label.

    sce_citeseq <- reducedDimSNF(sce_citeseq,
                                 method = "tSNE",
                                 dimNames = "tSNE_joint")

    g1 <- visualiseDim(sce_citeseq, dimNames = "tSNE_joint", colour_by = "SNF_W_clust") +
      labs(title = "tSNE (SNF clustering)")
    g2 <- visualiseDim(sce_citeseq, dimNames = "tSNE_joint",  colour_by = "ADT_clust") +
      labs(title = "tSNE (ADT clustering)")
    g3 <- visualiseDim(sce_citeseq, dimNames = "tSNE_joint",  colour_by = "RNA_clust") +
      labs(title = "tSNE (RNA clustering)")

    library(gridExtra)
    grid.arrange(g3, g2, g1, ncol = 2)

The expression of genes and proteins can be visualised by changing the `colour_by` parameter to assess the clusters. As an example, we highlight the plot by the RNA and ADT expression level of CD8.

    g1 <- visualiseDim(sce_citeseq, dimNames = "tSNE_joint",
                       colour_by = "hg19_CD8A",
                       data_from = "assay",
                       assay_name = "logcounts") +
      labs(title = "tSNE: hg19_CD8A (RNA expression)")
    g2 <- visualiseDim(sce_citeseq,dimNames = "tSNE_joint",
                       colour_by = "CD8",
                       data_from = "altExp",
                       altExp_assay_name = "logcounts") +
      labs(title = "tSNE: CD8 (ADT expression)")

    grid.arrange(g1, g2, ncol = 2)

## 3.4 Louvain clustering

As well as spectral clustering, CiteFuse can implement Louvain clustering if users wish to use another clustering method. We use the `igraph` package, and any community detection algorithms available in their package can be selected by changing the `method` parameter.

    SNF_W_louvain <- igraphClustering(sce_citeseq, method = "louvain")
    table(SNF_W_louvain)
    #> SNF_W_louvain
    #>   1   2   3   4   5   6   7
    #>  32  65  88  51 137  36  28

    sce_citeseq$SNF_W_louvain <- as.factor(SNF_W_louvain)

    visualiseDim(sce_citeseq, dimNames = "tSNE_joint", colour_by = "SNF_W_louvain") +
      labs(title = "tSNE (SNF louvain clustering)")

    visualiseKNN(sce_citeseq, colour_by = "SNF_W_louvain")

# 4 Differential Expression Analysis

## 4.1 Exploration of feature expression

CiteFuse has a wide range of visualisation tools to facilitate exploratory analysis of CITE-seq data. The `visualiseExprs` function is an easy-to-use function to generate boxplots, violinplots, jitter plots, density plots, and pairwise scatter/density plots of genes and proteins expressed in the data. The plots can be grouped by using the cluster labels stored in the `sce_citeseq` object.

    visualiseExprs(sce_citeseq,
                   plot = "boxplot",
                   group_by = "SNF_W_louvain",
                   feature_subset = c("hg19_CD2", "hg19_CD4", "hg19_CD8A", "hg19_CD19"))

    visualiseExprs(sce_citeseq,
                   plot = "violin",
                   group_by = "SNF_W_louvain",
                   feature_subset = c("hg19_CD2", "hg19_CD4", "hg19_CD8A", "hg19_CD19"))

    visualiseExprs(sce_citeseq,
                   plot = "jitter",
                   group_by = "SNF_W_louvain",
                   feature_subset = c("hg19_CD2", "hg19_CD4", "hg19_CD8A", "hg19_CD19"))

    visualiseExprs(sce_citeseq,
                   plot = "density",
                   group_by = "SNF_W_louvain",
                   feature_subset = c("hg19_CD2", "hg19_CD4", "hg19_CD8A", "hg19_CD19"))

    visualiseExprs(sce_citeseq,
                   altExp_name = "ADT",
                   group_by = "SNF_W_louvain",
                   plot = "violin", n = 5)

    visualiseExprs(sce_citeseq, altExp_name = "ADT",
                   plot = "jitter",
                   group_by = "SNF_W_louvain",
                   feature_subset = c("CD2", "CD8", "CD4", "CD19"))

    visualiseExprs(sce_citeseq, altExp_name = "ADT",
                   plot = "density",
                   group_by = "SNF_W_louvain",
                   feature_subset = c("CD2", "CD8", "CD4", "CD19"))

    visualiseExprs(sce_citeseq, altExp_name = "ADT",
                   plot = "pairwise",
                   feature_subset = c("CD4", "CD8"))
    #> number of iterations= 14
    #> number of iterations= 26

    visualiseExprs(sce_citeseq, altExp_name = "ADT",
                   plot = "pairwise",
                   feature_subset = c("CD45RA", "CD4", "CD8"),
                   threshold = rep(4, 3))

## 4.2 Perform DE Analysis with Wilcoxon Rank Sum test

CiteFuse also calculates differentially expressed (DE) genes through the `DEgenes` function. The cluster grouping to use must be specified in the `group` parameter. If `altExp_name` is not specified, RNA expression will be used as the default expression matrix.

Results form the DE analysis is stored in `sce_citeseq` as `DE_res_RNA_filter` and `DE_res_ADT_filter` for RNA and ADT expression, respectively.

### 4.2.1 For RNA expression

    # DE will be performed for RNA if altExp_name = "none"
    sce_citeseq <- DEgenes(sce_citeseq,
                           altExp_name = "none",
                           group = sce_citeseq$SNF_W_louvain,
                           return_all = TRUE,
                           exprs_pct = 0.5)

    sce_citeseq <- selectDEgenes(sce_citeseq,
                                 altExp_name = "none")
    datatable(format(do.call(rbind, metadata(sce_citeseq)[["DE_res_RNA_filter"]]),
                     digits = 2))

### 4.2.2 For ADT count

    sce_citeseq <- DEgenes(sce_citeseq,
                           altExp_name = "ADT",
                           group = sce_citeseq$SNF_W_louvain,
                           return_all = TRUE,
                           exprs_pct = 0.5)

    sce_citeseq <- selectDEgenes(sce_citeseq,
                                 altExp_name = "ADT")
    datatable(format(do.call(rbind, metadata(sce_citeseq)[["DE_res_ADT_filter"]]),
                     digits = 2))

## 4.3 Visualising DE Results

The DE genes can be visualised with the `DEbubblePlot` and `DEcomparisonPlot`. In each case, the gene names must first be extracted from the DE result objects.

### 4.3.1 circlepackPlot

The `circlepackPlot` takes a list of all DE genes from RNA and ADT DE analysis and will plot only the top most significant DE genes to plot.

    rna_DEgenes <- metadata(sce_citeseq)[["DE_res_RNA_filter"]]
    adt_DEgenes <- metadata(sce_citeseq)[["DE_res_ADT_filter"]]

    rna_DEgenes <- lapply(rna_DEgenes, function(x){
      x$name <- gsub("hg19_", "", x$name)
      x})
    DEbubblePlot(list(RNA = rna_DEgenes, ADT = adt_DEgenes))

### 4.3.2 DEcomparisonPlot

For the `DEcomparisonPlot`, as well as a list containing the DE genes for RNA and ADT, a `feature_list` specifying the genes and proteins of interest is required.

    rna_list <- c("hg19_CD4",
                  "hg19_CD8A",
                  "hg19_HLA-DRB1",
                  "hg19_ITGAX",
                  "hg19_NCAM1",
                  "hg19_CD27",
                  "hg19_CD19")

    adt_list <- c("CD4", "CD8", "MHCII (HLA-DR)", "CD11c", "CD56", "CD27", "CD19")

    rna_DEgenes_all <- metadata(sce_citeseq)[["DE_res_RNA"]]
    adt_DEgenes_all <- metadata(sce_citeseq)[["DE_res_ADT"]]

    feature_list <- list(RNA = rna_list, ADT = adt_list)
    de_list <- list(RNA = rna_DEgenes_all, ADT = adt_DEgenes_all)

    DEcomparisonPlot(de_list = de_list,
                     feature_list = feature_list)

# 5 ADT Importance Evaluation

An important evaluation in CITE-seq data analysis is to assess the quality of each ADT and to evaluate the contribution of ADTs towards clustering outcome. CiteFuse calculates the relative importance of ADT towards clustering outcome by using a random forest model. The higher the score of an ADT, the greater its importance towards the final clustering outcome.

    set.seed(2020)
    sce_citeseq <- importanceADT(sce_citeseq,
                                 group = sce_citeseq$SNF_W_louvain,
                                 subsample = TRUE)

    visImportance(sce_citeseq, plot = "boxplot")

    visImportance(sce_citeseq, plot = "heatmap")

    sort(metadata(sce_citeseq)[["importanceADT"]], decreasing = TRUE)[1:20]
    #>              CD27               CD8               CD4              CD28
    #>         40.833831         37.153051         32.351592         12.295992
    #>               CD5      PECAM (CD31)               CD7             CD11b
    #>         11.590230         11.563673         11.318339         10.820221
    #> IL7Ralpha (CD127)               CD2    MHCII (HLA-DR)              CD44
    #>          9.688986          9.096942          8.642220          8.354667
    #>      CD366 (tim3)         HLA-A,B,C             CD11c               CD3
    #>          5.936975          5.543170          5.085463          4.542166
    #>            CD45RA  CD26 (Adenosine)      PD-1 (CD279)              CD69
    #>          4.033837          3.861237          3.812794          3.733581

The importance scores can be visualised in a boxplot and heatmap. Our evaluation of ADT importance show that unsurprisingly CD4 and CD8 are the top two discriminating proteins in PBMCs.

Let us try clustering with only ADTs with a score greater than 5.

    subset_adt <- names(which(metadata(sce_citeseq)[["importanceADT"]] > 5))
    subset_adt
    #>  [1] "CD11b"             "CD11c"             "CD2"
    #>  [4] "CD27"              "CD28"              "CD366 (tim3)"
    #>  [7] "CD4"               "CD44"              "CD5"
    #> [10] "CD7"               "CD8"               "HLA-A,B,C"
    #> [13] "IL7Ralpha (CD127)" "MHCII (HLA-DR)"    "PECAM (CD31)"

    system.time(sce_citeseq <- CiteFuse(sce_citeseq,
                                        ADT_subset = subset_adt,
                                        metadata_names = c("W_SNF_adtSubset1",
                                                           "W_ADT_adtSubset1",
                                                           "W_RNA")))
    #> Calculating affinity matrix
    #> Performing SNF
    #>    user  system elapsed
    #>    3.88    0.08    3.95

    SNF_W_clust_adtSubset1 <- spectralClustering(metadata(sce_citeseq)[["W_SNF_adtSubset1"]], K = 5)
    #> Computing Spectral Clustering
    sce_citeseq$SNF_W_clust_adtSubset1 <- as.factor(SNF_W_clust_adtSubset1$labels)

    library(mclust)
    adjustedRandIndex(sce_citeseq$SNF_W_clust_adtSubset1, sce_citeseq$SNF_W_clust)
    #> [1] 0.8646855

When we compare between the two clustering outcomes, we find that the adjusted rand index is approximately 0.93, where a value of 1 denotes complete concordance.

# 6 Gene - ADT network

The `geneADTnetwork` function plots an interaction network between genes identified from the DE analysis. The nodes denote proteins and RNA whilst the edges denote positive and negative correlation in expression.

    RNA_feature_subset <- unique(as.character(unlist(lapply(rna_DEgenes_all, "[[", "name"))))
    ADT_feature_subset <- unique(as.character(unlist(lapply(adt_DEgenes_all, "[[", "name"))))

    geneADTnetwork(sce_citeseq,
                   RNA_feature_subset = RNA_feature_subset,
                   ADT_feature_subset = ADT_feature_subset,
                   cor_method = "pearson",
                   network_layout = igraph::layout_with_fr)

    #> IGRAPH 0d46a57 UN-B 72 134 --
    #> + attr: name (v/c), label (v/c), class (v/c), type (v/l), shape (v/c),
    #> | color (v/c), size (v/n), label.cex (v/n), label.color (v/c), value
    #> | (e/n), color (e/c), weights (e/n)
    #> + edges from 0d46a57 (vertex names):
    #>  [1] RNA_hg19_CCL5   --ADT_CD4
    #>  [2] RNA_hg19_NKG7   --ADT_CD4
    #>  [3] RNA_hg19_CST7   --ADT_CD4
    #>  [4] RNA_hg19_GNLY   --ADT_CD4
    #>  [5] RNA_hg19_CD8A   --ADT_CD4
    #>  [6] RNA_hg19_GZMB   --ADT_CD4
    #> + ... omitted several edges

# 7 RNA Ligand - ADT Receptor Analysis

With the advent of CITE-seq, we can now predict ligand-receptor interactions by using cell surface protein expression. CiteFuse implements a `ligandReceptorTest` to find ligand receptor interactions between sender and receiver cells. Importantly, the ADT count is used to predict receptor expression within receiver cells. Note that the setting `altExp_name = "RNA"` would enable users to predict ligand-receptor interaction from RNA expression only.

    data("lr_pair_subset", package = "CiteFuse")
    head(lr_pair_subset)
    #>      [,1]          [,2]
    #> [1,] "hg19_IL17RA" "CD45"
    #> [2,] "hg19_FAS"    "CD11b"
    #> [3,] "hg19_GZMK"   "CD62L"
    #> [4,] "hg19_CD40LG" "CD11b"
    #> [5,] "hg19_FLT3LG" "CD62L"
    #> [6,] "hg19_GZMA"   "CD19"

    sce_citeseq <- normaliseExprs(sce = sce_citeseq,
                                  altExp_name = "ADT",
                                  transform = "zi_minMax")

    sce_citeseq <- normaliseExprs(sce = sce_citeseq,
                                  altExp_name = "none",
                                  exprs_value = "logcounts",
                                  transform = "minMax")

    sce_citeseq <- ligandReceptorTest(sce = sce_citeseq,
                                      ligandReceptor_list = lr_pair_subset,
                                      cluster = sce_citeseq$SNF_W_louvain,
                                      RNA_exprs_value = "minMax",
                                      use_alt_exp = TRUE,
                                      altExp_name = "ADT",
                                      altExp_exprs_value = "zi_minMax",
                                      num_permute = 1000)
    #> 100 ......200 ......300 ......400 ......500 ......600 ......700 ......800 ......900 ......1000 ......

    visLigandReceptor(sce_citeseq,
                      type = "pval_heatmap",
                      receptor_type = "ADT")

    visLigandReceptor(sce_citeseq,
                      type = "pval_dotplot",
                      receptor_type = "ADT")

    visLigandReceptor(sce_citeseq,
                      type = "group_network",
                      receptor_type = "ADT")

    visLigandReceptor(sce_citeseq,
                      type = "group_heatmap",
                      receptor_type = "ADT")

    visLigandReceptor(sce_citeseq,
                      type = "lr_network",
                      receptor_type = "ADT")

# 8 Between-sample analysis

Lastly, we will jointly analyse the current PBMC CITE-seq data, taken from healthy human donors, and another subset of CITE-seq data from patients with cutaneous T-cell lymphoma (CTCL), again from Mimitou et al. (2019). The data `sce_ctcl_subset` provided in our `CiteFuse` package already contains the clustering information.

    data("sce_ctcl_subset", package = "CiteFuse")

To visualise and compare gene or protein expression data, we can use `visualiseExprsList` function.

    visualiseExprsList(sce_list = list(control = sce_citeseq,
                                       ctcl = sce_ctcl_subset),
                       plot = "boxplot",
                       altExp_name = "none",
                       exprs_value = "logcounts",
                       feature_subset = c("hg19_S100A10", "hg19_CD8A"),
                       group_by = c("SNF_W_louvain", "SNF_W_louvain"))

    visualiseExprsList(sce_list = list(control = sce_citeseq,
                                       ctcl = sce_ctcl_subset),
                       plot = "boxplot",
                       altExp_name = "ADT",
                       feature_subset = c("CD19", "CD8"),
                       group_by = c("SNF_W_louvain", "SNF_W_louvain"))

We can then perform differential expression analysis of the RNA expression level across the two clusters that have high CD19 expression in ADT.

    de_res <- DEgenesCross(sce_list = list(control = sce_citeseq,
                                           ctcl = sce_ctcl_subset),
                           colData_name = c("SNF_W_louvain", "SNF_W_louvain"),
                           group_to_test = c("2", "6"))
    de_res_filter <- selectDEgenes(de_res = de_res)
    de_res_filter
    #> $control
    #>             stats.W         pval     p.adjust meanExprs.1 meanExprs.2
    #> hg19_GNLY     737.0 3.794736e-07 3.794736e-07  0.07220401   1.6760852
    #> hg19_GZMH     881.5 8.746317e-06 8.746317e-06  0.00000000   0.7784718
    #> hg19_CST7     802.5 1.654137e-05 1.654137e-05  0.22512314   1.2305384
    #> hg19_ADGRE5   847.0 1.899854e-05 1.899854e-05  0.05559960   0.5042437
    #> hg19_ZEB2     967.5 6.934945e-05 6.934945e-05  0.00000000   0.4358523
    #> hg19_ITGB2    793.0 8.992606e-05 8.992606e-05  0.44185664   1.1053787
    #> hg19_NKG7     821.0 9.749228e-05 9.749228e-05  0.51066894   1.8319464
    #> hg19_CD99     801.0 1.345244e-04 1.345244e-04  0.58551572   1.1849517
    #> hg19_CCL4    1010.5 1.861793e-04 1.861793e-04  0.00000000   0.4388688
    #> hg19_RNF166   925.0 2.809741e-04 2.809741e-04  0.15303165   0.5466327
    #>              meanPct.1 meanPct.2  meanDiff   pctDiff        name   group
    #> hg19_GNLY   0.02325581 0.4923077 1.6038812 0.4690519   hg19_GNLY control
    #> hg19_GZMH   0.00000000 0.3692308 0.7784718 0.3692308   hg19_GZMH control
    #> hg19_CST7   0.13953488 0.5230769 1.0054153 0.3835420   hg19_CST7 control
    #> hg19_ADGRE5 0.06976744 0.4461538 0.4486441 0.3763864 hg19_ADGRE5 control
    #> hg19_ZEB2   0.00000000 0.3076923 0.4358523 0.3076923   hg19_ZEB2 control
    #> hg19_ITGB2  0.44186047 0.7076923 0.6635220 0.2658318  hg19_ITGB2 control
    #> hg19_NKG7   0.30232558 0.6000000 1.3212775 0.2976744   hg19_NKG7 control
    #> hg19_CD99   0.46511628 0.7846154 0.5994359 0.3194991   hg19_CD99 control
    #> hg19_CCL4   0.00000000 0.2769231 0.4388688 0.2769231   hg19_CCL4 control
    #> hg19_RNF166 0.09302326 0.4461538 0.3936010 0.3531306 hg19_RNF166 control
    #>
    #> $ctcl
    #>                stats.W         pval     p.adjust meanExprs.1 meanExprs.2
    #> hg19_RPS26         0.0 1.607680e-18 1.607680e-18   1.1562711   5.0724705
    #> hg19_HBB         812.5 1.590275e-08 1.590275e-08   0.0000000   0.4310260
    #> hg19_LEPROTL1    585.0 1.879911e-07 1.879911e-07   0.5706986   1.3007454
    #> hg19_LINC02446   753.0 3.615473e-06 3.615473e-06   0.2753145   1.0040831
    #> hg19_FOS         776.5 1.208420e-05 1.208420e-05   0.3562232   1.1253374
    #> hg19_TMEM123     778.0 5.002093e-05 5.002093e-05   0.4784153   1.0120341
    #> hg19_PPDPF       774.0 8.616326e-05 8.616326e-05   0.9165180   1.5015294
    #> hg19_LTB         787.0 1.106547e-04 1.106547e-04   1.1302569   1.9787302
    #> hg19_UXT         795.0 1.150403e-04 1.150403e-04   0.6335830   1.1605302
    #> hg19_ATP5IF1     852.0 2.127103e-04 2.127103e-04   0.3641278   0.7876142
    #>                meanPct.1 meanPct.2  meanDiff   pctDiff           name group
    #> hg19_RPS26     0.7538462 1.0000000 3.9161994 0.2461538     hg19_RPS26  ctcl
    #> hg19_HBB       0.0000000 0.4186047 0.4310260 0.4186047       hg19_HBB  ctcl
    #> hg19_LEPROTL1  0.4769231 0.9069767 0.7300469 0.4300537  hg19_LEPROTL1  ctcl
    #> hg19_LINC02446 0.2153846 0.6279070 0.7287686 0.4125224 hg19_LINC02446  ctcl
    #> hg19_FOS       0.2461538 0.6511628 0.7691142 0.4050089       hg19_FOS  ctcl
    #> hg19_TMEM123   0.4307692 0.7674419 0.5336189 0.3366726   hg19_TMEM123  ctcl
    #> hg19_PPDPF     0.7230769 0.9069767 0.5850115 0.1838998     hg19_PPDPF  ctcl
    #> hg19_LTB       0.6153846 0.9069767 0.8484732 0.2915921       hg19_LTB  ctcl
    #> hg19_UXT       0.5384615 0.8372093 0.5269472 0.2987478       hg19_UXT  ctcl
    #> hg19_ATP5IF1   0.3538462 0.6511628 0.4234864 0.2973166   hg19_ATP5IF1  ctcl

# 9 Read data from 10X Genomics

Readers unfamiliar with the workflow of converting a count matrix into a `SingleCellExperiment` object may use the `readFrom10X` function to convert count matrix from a 10X experiment into an object that can be used for all functions in CiteFuse.

    tmpdir <- tempdir()
    download.file("http://cf.10xgenomics.com/samples/cell-exp/3.1.0/connect_5k_pbmc_NGSC3_ch1/connect_5k_pbmc_NGSC3_ch1_filtered_feature_bc_matrix.tar.gz", file.path(tmpdir, "/5k_pbmc_NGSC3_ch1_filtered_feature_bc_matrix.tar.gz"))
    untar(file.path(tmpdir, "5k_pbmc_NGSC3_ch1_filtered_feature_bc_matrix.tar.gz"),
          exdir = tmpdir)
    sce_citeseq_10X <- readFrom10X(file.path(tmpdir, "filtered_feature_bc_matrix/"))
    sce_citeseq_10X

# 10 SessionInfo

    sessionInfo()
    #> R version 4.0.5 (2021-03-31)
    #> Platform: x86_64-w64-mingw32/x64 (64-bit)
    #> Running under: Windows Server 2012 R2 x64 (build 9600)
    #>
    #> Matrix products: default
    #>
    #> locale:
    #> [1] LC_COLLATE=C
    #> [2] LC_CTYPE=English_United States.1252
    #> [3] LC_MONETARY=English_United States.1252
    #> [4] LC_NUMERIC=C
    #> [5] LC_TIME=English_United States.1252
    #>
    #> attached base packages:
    #> [1] parallel  stats4    stats     graphics  grDevices utils     datasets
    #> [8] methods   base
    #>
    #> other attached packages:
    #>  [1] mclust_5.4.7                gridExtra_2.3
    #>  [3] DT_0.17                     scater_1.18.6
    #>  [5] ggplot2_3.3.3               SingleCellExperiment_1.12.0
    #>  [7] SummarizedExperiment_1.20.0 Biobase_2.50.0
    #>  [9] GenomicRanges_1.42.0        GenomeInfoDb_1.26.7
    #> [11] IRanges_2.24.1              S4Vectors_0.28.1
    #> [13] BiocGenerics_0.36.0         MatrixGenerics_1.2.1
    #> [15] matrixStats_0.58.0          CiteFuse_1.2.1
    #> [17] BiocStyle_2.18.1
    #>
    #> loaded via a namespace (and not attached):
    #>   [1] ggbeeswarm_0.6.0          Rtsne_0.15
    #>   [3] colorspace_2.0-0          ellipsis_0.3.1
    #>   [5] ggridges_0.5.3            scuttle_1.0.4
    #>   [7] bluster_1.0.0             XVector_0.30.0
    #>   [9] BiocNeighbors_1.8.2       farver_2.1.0
    #>  [11] graphlayouts_0.7.1        ggrepel_0.9.1
    #>  [13] RSpectra_0.16-0           fansi_0.4.2
    #>  [15] splines_4.0.5             sparseMatrixStats_1.2.1
    #>  [17] knitr_1.31                polyclip_1.10-0
    #>  [19] jsonlite_1.7.2            kernlab_0.9-29
    #>  [21] pheatmap_1.0.12           uwot_0.1.10
    #>  [23] ggforce_0.3.3             BiocManager_1.30.12
    #>  [25] compiler_4.0.5            dqrng_0.2.1
    #>  [27] assertthat_0.2.1          Matrix_1.3-2
    #>  [29] limma_3.46.0              tweenr_1.0.2
    #>  [31] BiocSingular_1.6.0        htmltools_0.5.1.1
    #>  [33] tools_4.0.5               rsvd_1.0.3
    #>  [35] igraph_1.2.6              gtable_0.3.0
    #>  [37] glue_1.4.2                GenomeInfoDbData_1.2.4
    #>  [39] reshape2_1.4.4            dplyr_1.0.5
    #>  [41] Rcpp_1.0.6                jquerylib_0.1.3
    #>  [43] vctrs_0.3.7               rhdf5filters_1.2.0
    #>  [45] crosstalk_1.1.1           DelayedMatrixStats_1.12.3
    #>  [47] ggraph_2.0.5              xfun_0.22
    #>  [49] stringr_1.4.0             beachmat_2.6.4
    #>  [51] lifecycle_1.0.0           irlba_2.3.3
    #>  [53] statmod_1.4.35            edgeR_3.32.1
    #>  [55] zlibbioc_1.36.0           MASS_7.3-53.1
    #>  [57] scales_1.1.1              tidygraph_1.2.0
    #>  [59] rhdf5_2.34.0              RColorBrewer_1.1-2
    #>  [61] yaml_2.2.1                sass_0.3.1
    #>  [63] segmented_1.3-3           stringi_1.5.3
    #>  [65] highr_0.8                 randomForest_4.6-14
    #>  [67] scran_1.18.6              BiocParallel_1.24.1
    #>  [69] rlang_0.4.10              pkgconfig_2.0.3
    #>  [71] bitops_1.0-6              evaluate_0.14
    #>  [73] lattice_0.20-41           purrr_0.3.4
    #>  [75] Rhdf5lib_1.12.1           labeling_0.4.2
    #>  [77] htmlwidgets_1.5.3         cowplot_1.1.1
    #>  [79] tidyselect_1.1.0          plyr_1.8.6
    #>  [81] magrittr_2.0.1            bookdown_0.21
    #>  [83] R6_2.5.0                  magick_2.7.1
    #>  [85] generics_0.1.0            DelayedArray_0.16.3
    #>  [87] DBI_1.1.1                 withr_2.4.1
    #>  [89] pillar_1.6.0              mixtools_1.2.0
    #>  [91] survival_3.2-10           RCurl_1.98-1.3
    #>  [93] tibble_3.1.0              crayon_1.4.1
    #>  [95] utf8_1.2.1                rmarkdown_2.7
    #>  [97] viridis_0.5.1             locfit_1.5-9.4
    #>  [99] grid_4.0.5                FNN_1.1.3
    #> [101] propr_4.2.6               digest_0.6.27
    #> [103] tidyr_1.1.3               dbscan_1.1-6
    #> [105] munsell_0.5.0             beeswarm_0.3.1
    #> [107] viridisLite_0.4.0         vipor_0.4.5
    #> [109] bslib_0.2.4