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(cowplot)
library(Seurat)

Step 1 Load data from 10X folder

Data can be downloaded from 10X website https://support.10xgenomics.com/single-cell-gene-expression/datasets/3.0.0/pbmc_1k_protein_v3

# Load in the RNA UMI matrix
cbmc.data <- readDataFrom10X(dir = "../../../Data/1K/filtered_feature_bc_matrix/")

# remove all Ig Proteins
cbmc.data$AntibodyCapture <- cbmc.data$AntibodyCapture[1:14,] 

Step 2 Create object

t1 <- Sys.time()
cbmc <- createObject(data = cbmc.data)
t2 <- Sys.time()
t2 - t1

Step 3 Pre-process

1) Filter out unwanted cells (optional)

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:15)
#ElbowPlot(cbmc)

Step 6 Distance calculation and joint distance calculation

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

t1 <- Sys.time()
cbmc <- jointDistance(object = cbmc, keep.rna = TRUE, keep.adt = TRUE)
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(2,2,2))
t2 <- Sys.time()
t2 - t1

# contribution of two modalities
distHeatMap(object = cbmc, label = TRUE)

Step 9 Visualization ADT vs RNA vs Joint

1) Cell clusters

#gridDimPlot(cbmc, wide.rel = 1.5, legend = FALSE, reduction.prefix = "tsne_", height.rel = 0.5)

plots <- generateGridDimPlot(cbmc, legend = FALSE, darkTheme = FALSE,cluster.lable.size = 8 )
listPlot(object = plots, labels = "")

###### 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 = 3, adt.label = TRUE, rna.label = TRUE)
# Heatmap for RNA clusters
heatMapPlot(object = cbmc, group.by = "rnaClusterID", height.rel = 3, adt.label = TRUE, rna.label = TRUE)
# Heatmap for ADT clusters
heatMapPlot(object = cbmc, group.by = "adtClusterID", height.rel = 3, adt.label = TRUE, rna.label = TRUE)

p1 <- VlnPlot(cbmc, features = "adt_CD4", pt.size = 0, group.by = 'jointClusterID') + NoLegend() 
p2 <- VlnPlot(cbmc, features = "adt_CD3", pt.size = 0, group.by = 'jointClusterID') + NoLegend() 
p3 <- VlnPlot(cbmc, features = "adt_CD16", pt.size = 0, group.by = 'jointClusterID') + NoLegend() 
p4 <- VlnPlot(cbmc, features = "adt_CD56", pt.size = 0, group.by = 'jointClusterID') + NoLegend()
p5 <- VlnPlot(cbmc, features = "adt_CD45RA", pt.size = 0, group.by = 'jointClusterID') + NoLegend() 
p6 <- VlnPlot(cbmc, features = "adt_CD45RO", pt.size = 0, group.by = 'jointClusterID') + NoLegend() 
plot_grid(p1,p2,p3,p4,p5,p6,nrow = 1)

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