In WilsonImmunologyLab/LinQView: Integrate information from both transcriptome and cell surface protein signals to better identify cell heterogeneity

options(tinytex.verbose = TRUE)
knitr::opts_chunk$set(
  cache = TRUE,
  cache.lazy = FALSE,
  tidy = TRUE
)

Load packages

library(LinQView)
library(Seurat)
library(cowplot)
library(ggplot2)

CITE-seq dataset

In this tutorial, we use a public CITE-seq dataset to illustrate Joint analysis using LinQ-seq. Data could be download from NCBI: RNA (ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE100nnn/GSE100866/suppl/GSE100866_CBMC_8K_13AB_10X-RNA_umi.csv.gz) ADT (ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE100nnn/GSE100866/suppl/GSE100866_CBMC_8K_13AB_10X-ADT_umi.csv.gz)

Step 1 Load data from 10X folder

# Load in the RNA UMI matrix
cbmc.rna <- as.sparse(x = read.csv(file = "../../../Data/citeseq/GSE100866_CBMC_8K_13AB_10X-RNA_umi.csv.gz", sep = ",", header = TRUE, row.names = 1))
cbmc.rna <- CollapseSpeciesExpressionMatrix(object = cbmc.rna)
# Load in the ADT UMI matrix
cbmc.adt <- as.sparse(x = read.csv(file = "../../../Data/citeseq/GSE100866_CBMC_8K_13AB_10X-ADT_umi.csv.gz", sep = ",", header = TRUE, row.names = 1))
cbmc.adt <- cbmc.adt[setdiff(x = rownames(x = cbmc.adt), y = c("CCR5", "CCR7", "CD10")), ]

Step 2 Create object

t1 <- Sys.time()
cbmc <- createObject(rna = cbmc.rna, adt = cbmc.adt)
t2 <- Sys.time()
t2 - t1

Step 3 Pre-process

1) Filter out unwanted cells (optional)

for this dataset, we don't need to filter out unwanted cells

cbmc <- subset(cbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < softThreshold(cbmc))

2) Remove unwanted genes (optional)

for this dataset, we don't need to filter out unwanted genes

# remove Ig genes
#cbmc <- removeGene(object = cbmc,pattern = '^IG[HKL]')

3) Normalization

data Normalization for both ADT (CLR) and RNA (log)

t1 <- Sys.time()
cbmc <- dataNormalization(object = cbmc)
t2 <- Sys.time()
t2 - t1

4) Indentify HVGs for RNA data

Call seurat function to identify highly variable genes (HVG) for RNA data

t1 <- Sys.time()
cbmc <- FindVariableFeatures(object = cbmc)   # directly use Seurat Function
t2 <- Sys.time()
t2 - t1

5) Data scaling

Scale data for both ADT and RNA

t1 <- Sys.time()
cbmc <- dataScaling(object = cbmc) 
t2 <- Sys.time()
t2 - t1

Step 4 Linear dimension reduction (PCA)

directly call Seurat function for linear dimension reduction (PCA)

t1 <- Sys.time()
cbmc <- RunPCA(cbmc, features = VariableFeatures(object = cbmc), verbose = FALSE)   # directly use Seurat Function
t2 <- Sys.time()
t2 - t1

Step 5 Determine number of PCs

call Seurat function JackStraw to determine number of PCs

#cbmc <- JackStraw(cbmc, num.replicate = 100)
#cbmc <- ScoreJackStraw(cbmc, dims = 1:20)
#JackStrawPlot(cbmc, dims = 1:20)
#ElbowPlot(cbmc)

Step 6 Distance calculation and joint distance calculation

calculate cell-cell distances for RNA, ADT and joint. number of PC was set to 20 by default.

t1 <- Sys.time()
cbmc <- jointDistance(object = cbmc, keep.rna = TRUE, keep.adt = TRUE, dims = 25)
t2 <- Sys.time()
t2 - t1

Step 7 Non-linear dimension reduction (UMAP and t-SNE)

run UMAP as Non-linear dimension reduction for RNA, ADT and joint analysis.

t1 <- Sys.time()
cbmc <- tsneFromDistane(object = cbmc, assay = "All")
t2 <- Sys.time()
t2 - t1

Step 8 Clustering

t1 <- Sys.time()
cbmc <- clusteringFromDistance(object = cbmc, assay = "All", resolution = c(0.9,0.9,0.9))
t2 <- Sys.time()
t2 - t1
# contribution of two modalities
distHeatMap(object = cbmc)

Step 9 Visualization ADT vs RNA vs Joint

1) Cell clusters

plots <- generateGridDimPlot(cbmc, legend = FALSE, darkTheme = FALSE)

listPlot(object = plots, align = "h")

###### user also can only plot some of those plots by index, figure ident or figure map info
#listPlot(object = plots, fig.ident = "RNA")
#listPlot(object = plots, fig.ident = "RNA", fig.map = "RNA")
#user can use plotInfo() function to get index, figure ident and figure map information, then plot figures by index
plotInfo(plots)
#listPlot(object = plots, fig.id = 1)

2) Heat maps

# Heatmap for joint clusters
heatMapPlot(object = cbmc, group.by = "jointClusterID", height.rel = 1, adt.label = TRUE)
# Heatmap for RNA clusters
heatMapPlot(object = cbmc, group.by = "rnaClusterID", height.rel = 1, adt.label = TRUE)
# Heatmap for ADT clusters
heatMapPlot(object = cbmc, group.by = "adtClusterID", height.rel = 1, adt.label = TRUE)

3) RNA and ADT expression

VlnPlot(cbmc, features = c("rna_CD8A", "adt_CD8", "rna_NCAM1", "adt_CD56", "rna_CD3G", "adt_CD3"), group.by = 'jointClusterID', pt.size = 0, ncol = 2)

WilsonImmunologyLab/LinQView documentation built on Jan. 3, 2022, 10 p.m.