Introduction

Single-cell RNA sequencing is becoming an increasingly common tool to investigate the cellular population and patients’ outcome in cancer research. However, due to the sparse data and the complex tumor microenvironment, it is challenging to identify neoplastic cells that play important roles in tumor growth and disease progression. This challenge is exaggerated in the research of blood cancer patients, from whom the neoplastic cells can be highly similar to the normal cells.

In this package we present partCNV/partCNVH, a statistical framework for rapid and accurate detection of aneuploid cells with local copy number deletion or amplification. Our method uses an EM algorithm with mixtures of Poisson distributions while incorporating cytogenetics information to guide the classification (partCNV). When applicable, we further improve the accuracy by integrating a Hidden Markov Model for feature selection (partCNVH).

Installation

The package can be installed using BiocManager by the following commands

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("partCNV")

Alternatively, the package can be installed using devtools and launched by

library(devtools)
install_github("rx-li/partCNV")

Understand your cytogenetics data

First of all, our method can only be used to infer cell status if the patient has a part of the cells with cytogenetics alterations. For example, if a patient has cytogenetics data as 46,XY,del(20)(q11.1q13.1)[5]/46,XY[15], it means that out of 20 metaphases, 5 (i.e., 25%) of the cells has deletion on chromosome 20 in region q11.1 to q13.1. If the patient has a normal cytogenetics, e.g., 46,XY[20], or all altered cells, e.g., 46,XY,del(20)(q11.1q13.1)[25], there won’t be any need to apply the proposed method.

Second, when you have a complicated cytogenetics feature, use them one by one to identify the desired cell group. For example, 47,XY,+8[5]/46,idem,del(8)(q21.2q24.3)/46,XY[7], we start with the chromosome amplification on chromosome 8 excluding the region (q21.2q24.3) to identify the cells with this alteration using the proposed method. After the cells with chromosome 8 amplification are identified, we apply the proposed method to identify cells with del(8)(q21.2q24.3).

Lastly, our tool is primarily designed for analyzing the scRNA-seq data from patients with hematologic malignancies, such as MDS and AML. For solid tumors such as lung cancer and breast cancer, it is better to use CNV based method such as copyKAT and inferCNV inferCNV.

Load example data

The analysis can start from whole-genome scRNA-seq data or a subset of scRNA-seq data based on the location of interest. In this package, we prepared a whole-genome scRNA-seq data:

library(partCNV)
data(SimData)
dim(SimData)

## [1] 24519   400

SimData[1:5,1:5]

## 5 x 5 sparse Matrix of class "dgCMatrix"
##            C1 C2 C3 C4 C5
## AL627309.1  .  .  .  .  .
## AL627309.5  .  .  .  .  .
## AL627309.4  .  .  .  .  .
## AP006222.2  .  .  .  .  .
## LINC01409   .  .  .  .  .

This simple example dataset contains 24519 genes and 400 cells. If you start with Seurat object (after quality control), run:

library(Seurat)
Seurat_obj <- NormalizeData(Your_SeuratObj, normalization.method = "LogNormalize", scale.factor = 10000)
Counts = Seurat_obj@assays$RNA@counts

If you start with a SingleCellExperiment object, the counts can be normalized using function NormalizeCounts.

Counts <- NormalizeCounts(Your_SingleCellExperimentObj, scale_factor=10000)

Run partCNV

Since this is simulation data, prior knowledge (e.g., cytogenetics data in real studies) shows that this example data has 40% of cells with deletion on chromosome 20 q11.1 to q13.1. Let’s get started with locating this region.

res <- GetCytoLocation(cyto_feature = "chr20(q11.1-q13.1)")

## downloading 1 resources

## retrieving 1 resource

## loading from cache

## require("GenomicRanges")

## Interested region: chr20:28100001-51200000.

## A total of 381 genes are located in this region.

Then let’s subset the data and normalize the data:

GEout <- GetExprCountCyto(cytoloc_output = res, Counts = as.matrix(SimData), normalization = TRUE, qt_cutoff = 0.99)

For this function, the qt_cutoff is to filter out the cells with very low expressions. 0.99 here means that we filter out cells that only express in 1% (1-0.99) cells. Make this qt_cutoff larger if your total gene number within the region is small.

Now we apply partCNV:

pcout <- partCNV(int_counts = GEout$ProcessedCount,
                 cyto_type = "del",
                 cyto_p = 0.40)

Understand the results:

table(pcout)

## pcout
##   0   1 
## 242 158

sum(pcout==1)/length(pcout)

## [1] 0.395

p1 <- sum(pcout==1)/length(pcout)

39.5% of cells are labeled as locally aneuploid (1) and others are diploid (0).

Let’s visualize it:

library(Seurat)
sim_seurat <- CreateSeuratObject(counts = SimData)
sim_seurat <- NormalizeData(sim_seurat, normalization.method = "LogNormalize", scale.factor = 10000)
sim_seurat <- FindVariableFeatures(sim_seurat, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(sim_seurat)
sim_seurat <- ScaleData(sim_seurat, features = all.genes)
sim_seurat <- RunPCA(sim_seurat, features = VariableFeatures(object = sim_seurat))
sim_seurat <- RunUMAP(sim_seurat, dims = 1:10)
sim_seurat <- AddMetaData(
    object = sim_seurat,
    metadata = pcout,
    col.name = "partCNV_label"
)
sim_seurat$partCNV_label <- factor(sim_seurat$partCNV_label, levels = c(1,0), labels = c(
    "aneuploid", "diploid"
))

library(ggplot2)
DimPlot(sim_seurat, reduction = "umap", group = "partCNV_label") + ggtitle(paste0("partCNV (", signif(p1,2)*100, "%)"))

Run partCNVH

Compared with partCNV, partCNVH added an additional step of feature selection. This is especially helpful if your cytogenetics provide a very broad region and part of it does not have chromosomal alterations. The first two steps of using partCNVH are exactly the same as using partCNV.

res <- GetCytoLocation(cyto_feature = "chr20(q11.1-q13.1)")

## loading from cache

## Interested region: chr20:28100001-51200000.

## A total of 381 genes are located in this region.

GEout <- GetExprCountCyto(cytoloc_output = res, Counts = as.matrix(SimData), normalization = TRUE, qt_cutoff = 0.99)

Now we apply partCNVH:

pcHout <- partCNVH(int_counts = GEout$ProcessedCount,
                 cyto_type = "del",
                 cyto_p = 0.40)

Understand the results (in pcHout, EMlabel is the partCNV label and EMHMMlabel is the partCNVH label).

table(pcHout$EMHMMlabel)

## 
##   0   1 
## 235 165

sum(pcHout$EMHMMlabel==1)/length(pcHout$EMHMMlabel)

## [1] 0.4125

p2 <- sum(pcHout$EMHMMlabel==1)/length(pcHout$EMHMMlabel)

41.25% of cells are labeled as locally aneuploid (1) and others are diploid (0).