Tutorial: Single Cell RNA-Seq Analysis of NK, CD4+ and CD8+ T cells isolated from LCMV infected WT and Ifng - CTCF binding site mutant mice
1. Quality Control
This document takes the cellranger count data for two samples (a WT and an Ifng CBS mutant mice) and performs standard single cell rnaseq analysis.
# Install hdf5 library if its not already installed
#sudo apt-get install libhdf5-dev or brew install hdf5
#install.packages("hdf5r")
# Load below libraries
library(Seurat)
library(ggplot2)
library(plotly)
library(tidyverse)
library(cowplot)
library(scProportionTest)
library(rollama)
library(SCassist)
library(httr)
library(jsonlite)1.1. Data setup
# Set the input, output paths
setwd("/singlecellassistant/")
datapath="exampleData/LCMV"
outputpath="/singlecellassistant/exampleOutput_G/LCMV/"
# Load RData for example dataset
load("/singlecellassistant/exampleOutput_G/LCMV/LCMV.RData")
# Specify samples for analysis.
# Define LCMV Day4, CD4, CD8, NK samples
samples=c("WT","KO")
# Read the data file
WT <- Read10X_h5(file.path("/singlecellassistant/exampleData/LCMV/GSM6625298_scRNA_LCMV_Day4_CD4_CD8_NK_WT_filtered_feature_bc_matrix.h5"), use.names = T)
KO <- Read10X_h5(file.path("/singlecellassistant/exampleData/LCMV/GSM6625299_scRNA_LCMV_Day4_CD4_CD8_NK_KO_filtered_feature_bc_matrix.h5"), use.names = T)
# Add sample names to columnnames
colnames(WT$`Gene Expression`) <- paste(sapply(strsplit(colnames(WT$`Gene Expression`),split="-"),'[[',1L),"WT",sep="-")
colnames(KO$`Gene Expression`) <- paste(sapply(strsplit(colnames(KO$`Gene Expression`),split="-"),'[[',1L),"KO",sep="-")
# Create a Seurat data object from the gex matrix
WT <- CreateSeuratObject(counts = WT[["Gene Expression"]], names.field = 2,names.delim = "\\-")
KO <- CreateSeuratObject(counts = KO[["Gene Expression"]], names.field = 2,names.delim = "\\-")
# Merge seurat objects
allsamples <- merge(WT, KO, project = "LCMV")
allsamples## An object of class Seurat
## 32285 features across 10551 samples within 1 assay
## Active assay: RNA (32285 features, 0 variable features)
## 2 layers present: counts.1, counts.2head(allsamples)## orig.ident nCount_RNA nFeature_RNA
## AAACCCAAGAGTTGTA-WT WT 17880 4120
## AAACCCAAGATTTGCC-WT WT 7916 2395
## AAACCCAAGTTGTAGA-WT WT 35993 5094
## AAACCCACAAGTATCC-WT WT 9059 2865
## AAACCCAGTCCAGGTC-WT WT 10948 2939
## AAACCCATCCTGGGAC-WT WT 16768 3484
## AAACGAACAGAGAGGG-WT WT 10467 3159
## AAACGAACAGCGACAA-WT WT 12996 2853
## AAACGAACAGTGTGGA-WT WT 26674 4654
## AAACGAATCTAGCCTC-WT WT 11434 3100table(Idents(allsamples))##
## WT KO
## 5666 4885# Merge counts data
allsamples = JoinLayers(allsamples)1.2. Quality Analysis
This part of the workflow uses standard plots to identify potential filtering cutoff values.
# Generate Percent mito for each of the cells
grep("^mt-",rownames(allsamples@assays$RNA$counts),value = TRUE)## [1] "mt-Nd1" "mt-Nd2" "mt-Co1" "mt-Co2" "mt-Atp8" "mt-Atp6" "mt-Co3"
## [8] "mt-Nd3" "mt-Nd4l" "mt-Nd4" "mt-Nd5" "mt-Nd6" "mt-Cytb"allsamples$percent.mito <- PercentageFeatureSet(allsamples, pattern = "^mt-")
summary(allsamples$percent.mito)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.819 2.425 5.053 3.310 83.600head(allsamples)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCAAGAGTTGTA-WT WT 17880 4120 3.333333
## AAACCCAAGATTTGCC-WT WT 7916 2395 2.627590
## AAACCCAAGTTGTAGA-WT WT 35993 5094 2.008724
## AAACCCACAAGTATCC-WT WT 9059 2865 1.953858
## AAACCCAGTCCAGGTC-WT WT 10948 2939 2.831567
## AAACCCATCCTGGGAC-WT WT 16768 3484 0.942271
## AAACGAACAGAGAGGG-WT WT 10467 3159 3.668673
## AAACGAACAGCGACAA-WT WT 12996 2853 1.569714
## AAACGAACAGTGTGGA-WT WT 26674 4654 1.570818
## AAACGAATCTAGCCTC-WT WT 11434 3100 1.687948# Generate Percent ribo for each of the cells
allsamples$percent.ribo <- PercentageFeatureSet(allsamples, pattern = "^Rp[sl]")
summary(allsamples$percent.ribo)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.002 18.156 24.731 24.894 32.868 56.275head(allsamples)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCAAGAGTTGTA-WT WT 17880 4120 3.333333
## AAACCCAAGATTTGCC-WT WT 7916 2395 2.627590
## AAACCCAAGTTGTAGA-WT WT 35993 5094 2.008724
## AAACCCACAAGTATCC-WT WT 9059 2865 1.953858
## AAACCCAGTCCAGGTC-WT WT 10948 2939 2.831567
## AAACCCATCCTGGGAC-WT WT 16768 3484 0.942271
## AAACGAACAGAGAGGG-WT WT 10467 3159 3.668673
## AAACGAACAGCGACAA-WT WT 12996 2853 1.569714
## AAACGAACAGTGTGGA-WT WT 26674 4654 1.570818
## AAACGAATCTAGCCTC-WT WT 11434 3100 1.687948
## percent.ribo
## AAACCCAAGAGTTGTA-WT 22.80761
## AAACCCAAGATTTGCC-WT 16.28348
## AAACCCAAGTTGTAGA-WT 21.85147
## AAACCCACAAGTATCC-WT 16.81201
## AAACCCAGTCCAGGTC-WT 25.30142
## AAACCCATCCTGGGAC-WT 35.07276
## AAACGAACAGAGAGGG-WT 22.92921
## AAACGAACAGCGACAA-WT 42.22068
## AAACGAACAGTGTGGA-WT 26.22029
## AAACGAATCTAGCCTC-WT 28.44149# Generate Percent hemoglobin for each of the cells
allsamples$percent.hb <- PercentageFeatureSet(allsamples, pattern = "^Hb[^(p)]")
summary(allsamples$percent.hb)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.000000 0.000000 0.004744 0.005529 2.162162head(allsamples)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCAAGAGTTGTA-WT WT 17880 4120 3.333333
## AAACCCAAGATTTGCC-WT WT 7916 2395 2.627590
## AAACCCAAGTTGTAGA-WT WT 35993 5094 2.008724
## AAACCCACAAGTATCC-WT WT 9059 2865 1.953858
## AAACCCAGTCCAGGTC-WT WT 10948 2939 2.831567
## AAACCCATCCTGGGAC-WT WT 16768 3484 0.942271
## AAACGAACAGAGAGGG-WT WT 10467 3159 3.668673
## AAACGAACAGCGACAA-WT WT 12996 2853 1.569714
## AAACGAACAGTGTGGA-WT WT 26674 4654 1.570818
## AAACGAATCTAGCCTC-WT WT 11434 3100 1.687948
## percent.ribo percent.hb
## AAACCCAAGAGTTGTA-WT 22.80761 0.011185682
## AAACCCAAGATTTGCC-WT 16.28348 0.000000000
## AAACCCAAGTTGTAGA-WT 21.85147 0.002778318
## AAACCCACAAGTATCC-WT 16.81201 0.000000000
## AAACCCAGTCCAGGTC-WT 25.30142 0.009134088
## AAACCCATCCTGGGAC-WT 35.07276 0.005963740
## AAACGAACAGAGAGGG-WT 22.92921 0.009553836
## AAACGAACAGCGACAA-WT 42.22068 0.000000000
## AAACGAACAGTGTGGA-WT 26.22029 0.000000000
## AAACGAATCTAGCCTC-WT 28.44149 0.000000000# Check the total number of cells in this experiment
table(allsamples$orig.ident)##
## KO WT
## 4885 56661.3. Quality Filtering
Here we look at the distribution of the data on MT, Rb and Hb genes and identify the cutoff values to filter and visualize the before and after QC data.
# Generate before filtering quality plots
qcbefore=VlnPlot(allsamples,features = c("nFeature_RNA", "nCount_RNA","percent.mito","percent.ribo"),ncol = 2, pt.size = 0) +
NoLegend()
qcbefore
# Save the plot in current working directory
ggsave(paste(outputpath,"LCMV-before-qc-violineplot.pdf",sep=""), plot = qcbefore, width = 15, height = 20, units = "cm")
# Filter data based on the above plot.
allsamplesgood <- subset(allsamples, subset = nCount_RNA > 2000 & nCount_RNA < 37000 & nFeature_RNA > 500 & nFeature_RNA < 6500 & percent.mito < 8 & percent.ribo < 50 & percent.hb < 2)
# Counts before filtering
table(allsamples$orig.ident)##
## KO WT
## 4885 5666# Counts after filtering
table(allsamplesgood$orig.ident)##
## KO WT
## 4130 47401.4. Generate after QC plots
# Generate after filtering quality plots
qcafter=VlnPlot(allsamplesgood,features = c("nFeature_RNA", "nCount_RNA","percent.mito","percent.ribo"),ncol = 2, pt.size = 0) +
NoLegend()
qcafter
# Save the plot in current working directory
ggsave(paste(outputpath,"LCMV-after-qc-violineplot.pdf",sep=""), plot = qcafter, width = 15, height = 20, units = "cm")2. Normalization
Here we use the standard normalization method.
# Cell level normalization - accounts for sequencing depth
allsamplesgood <- NormalizeData(
object = allsamplesgood,
normalization.method = "LogNormalize",
scale.factor = 10000)3. Variable features
Here we identify variable genes and perform scaling, including regressing to remove any technical variations
##################### Find variable genes
# Returns top 2000 variable genes
# VST is uses LOESS method
# FindVariableFeatures needs normalization
allsamplesgood <- FindVariableFeatures(
object = allsamplesgood,
selection.method = "vst")
top10 <- head(VariableFeatures(allsamplesgood), 10)
top10## [1] "Ccl4" "Hist1h1b" "Xcl1" "Hist1h2ap" "Hist1h2ae" "Spp1"
## [7] "Hist1h3c" "Ccl3" "H2-Aa" "Cxcl10"## Scaling & Batch Correction
##################### SCALE DATA
# PCA needs scaled data
# Zero centers and scales (mean/sd) gene/feature data in each cell (for across sample comparison), so extreme ranges in expression do not affect clustering (done for making cells with similar expression cluster together)
dim(allsamplesgood)## [1] 32285 8870allsamplesgood <- ScaleData(
object = allsamplesgood,
# Scale all genes - not just variable genes
features = rownames(allsamplesgood),
# Remove unwanted sources of variation (technical, batch etc.)
vars.to.regress = c("percent.mito", "nFeature_RNA", "percent.ribo", "nCount_RNA"))
# List top 10 variable genes
top10 <- head(VariableFeatures(allsamplesgood), 10)
top10## [1] "Ccl4" "Hist1h1b" "Xcl1" "Hist1h2ap" "Hist1h2ae" "Spp1"
## [7] "Hist1h3c" "Ccl3" "H2-Aa" "Cxcl10"# Plot variable features, label the top 10, save the plot
vfp1 <- VariableFeaturePlot(allsamplesgood)
vfp1 <- LabelPoints(plot = vfp1, points = top10, repel = TRUE)
vfp1
ggsave(paste(outputpath,"LCMV-variablefeatureplot.pdf",sep=""), plot = vfp1, width = 20, height = 15, units = "cm")4. PCA
Here we perform PCA to use for downstream analysis
# Run PCA to generate 50 pcs
allsamplesgood <- RunPCA(object = allsamplesgood, npcs=50)
# Print genes in the top 5 pcs
print(allsamplesgood[["pca"]], dims = 1:5, nfeatures = 5)## PC_ 1
## Positive: Nusap1, Ccna2, Knl1, Top2a, Birc5
## Negative: Hopx, Bcl2a1b, Ctsb, Tnfrsf1b, Ly6a
## PC_ 2
## Positive: Klrb1c, Fcer1g, Nkg7, Klrd1, Plek
## Negative: Inpp4b, Cd6, Hif1a, Rgs10, Tnfrsf4
## PC_ 3
## Positive: Aif1, C1qa, Slc40a1, C1qb, Hmox1
## Negative: Hopx, Tnfrsf18, S100a10, Cd6, Sdf4
## PC_ 4
## Positive: E2f1, Cdc6, Lig1, Ccne2, Pcna
## Negative: Ube2c, Cenpf, Ccnb1, Cdc20, Aspm
## PC_ 5
## Positive: Gm42418, Ly6c1, Nsg2, Actn1, Sell
## Negative: Tnfrsf4, Tnfrsf18, Foxp3, Maf, Ctla44.1. Choosing appropriate number of PCs
Here we look at the elbow plot and identify the number of PCs to use for downstream analysis
# Identify number of PCs that explains majority of variations
ElbowPlot(allsamplesgood, ndims=50, reduction = "pca")
# Visualize the genes in the first PC
VizDimLoadings(allsamplesgood, dims = 1, ncol = 1) + theme_minimal(base_size = 8)
# Visualize the genes in the second PC
VizDimLoadings(allsamplesgood, dims = 2, ncol = 1) + theme_minimal(base_size = 8)
# Set the identity column to WT vs KO
Idents(allsamplesgood) <- "orig.ident"
# Visualize the cells after pca
DimPlot(object = allsamplesgood, reduction = "pca")
# Plot heatmap with cells=500 plotting cells with extreme cells on both ends of spectrum
DimHeatmap(object = allsamplesgood, dims = 1:5, cells = 500, balanced = TRUE)
5. FindNeighbors
Here we run FindNeighbors function, given the number of pcs (looking at the elbow plot) we are going to use for the downstream analysis.
# Run findneighbors
allsamplesgood <- FindNeighbors(allsamplesgood, dims = 1:20)6. Clustering
Here we do clustering, randomnly choosing some resolution numbers, to start with.
# Perform clustering using the different resolutions, to identify the numbers, to match the 23 clusters we got in the manuscript.
allsamplesgood <- FindClusters(allsamplesgood, resolution = seq(1.5,2,0.1))## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8870
## Number of edges: 310666
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7853
## Number of communities: 19
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8870
## Number of edges: 310666
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7785
## Number of communities: 21
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8870
## Number of edges: 310666
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7710
## Number of communities: 22
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8870
## Number of edges: 310666
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7656
## Number of communities: 23
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8870
## Number of edges: 310666
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7596
## Number of communities: 23
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8870
## Number of edges: 310666
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7533
## Number of communities: 26
## Elapsed time: 0 secondshead(allsamplesgood)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCAAGAGTTGTA-WT WT 17880 4120 3.333333
## AAACCCAAGATTTGCC-WT WT 7916 2395 2.627590
## AAACCCAAGTTGTAGA-WT WT 35993 5094 2.008724
## AAACCCACAAGTATCC-WT WT 9059 2865 1.953858
## AAACCCAGTCCAGGTC-WT WT 10948 2939 2.831567
## AAACCCATCCTGGGAC-WT WT 16768 3484 0.942271
## AAACGAACAGAGAGGG-WT WT 10467 3159 3.668673
## AAACGAACAGCGACAA-WT WT 12996 2853 1.569714
## AAACGAACAGTGTGGA-WT WT 26674 4654 1.570818
## AAACGAATCTAGCCTC-WT WT 11434 3100 1.687948
## percent.ribo percent.hb RNA_snn_res.1.5 RNA_snn_res.1.6
## AAACCCAAGAGTTGTA-WT 22.80761 0.011185682 14 15
## AAACCCAAGATTTGCC-WT 16.28348 0.000000000 0 0
## AAACCCAAGTTGTAGA-WT 21.85147 0.002778318 9 8
## AAACCCACAAGTATCC-WT 16.81201 0.000000000 15 16
## AAACCCAGTCCAGGTC-WT 25.30142 0.009134088 13 5
## AAACCCATCCTGGGAC-WT 35.07276 0.005963740 12 13
## AAACGAACAGAGAGGG-WT 22.92921 0.009553836 2 2
## AAACGAACAGCGACAA-WT 42.22068 0.000000000 6 4
## AAACGAACAGTGTGGA-WT 26.22029 0.000000000 9 8
## AAACGAATCTAGCCTC-WT 28.44149 0.000000000 12 13
## RNA_snn_res.1.7 RNA_snn_res.1.8 RNA_snn_res.1.9
## AAACCCAAGAGTTGTA-WT 16 17 16
## AAACCCAAGATTTGCC-WT 1 1 1
## AAACCCAAGTTGTAGA-WT 9 9 10
## AAACCCACAAGTATCC-WT 19 19 20
## AAACCCAGTCCAGGTC-WT 7 16 13
## AAACCCATCCTGGGAC-WT 12 14 15
## AAACGAACAGAGAGGG-WT 2 2 2
## AAACGAACAGCGACAA-WT 4 3 4
## AAACGAACAGTGTGGA-WT 9 9 10
## AAACGAATCTAGCCTC-WT 6 10 8
## RNA_snn_res.2 seurat_clusters
## AAACCCAAGAGTTGTA-WT 15 15
## AAACCCAAGATTTGCC-WT 1 1
## AAACCCAAGTTGTAGA-WT 16 16
## AAACCCACAAGTATCC-WT 22 22
## AAACCCAGTCCAGGTC-WT 18 18
## AAACCCATCCTGGGAC-WT 12 12
## AAACGAACAGAGAGGG-WT 0 0
## AAACGAACAGCGACAA-WT 2 2
## AAACGAACAGTGTGGA-WT 17 17
## AAACGAATCTAGCCTC-WT 12 12# Count number of clusters at each resolution
sapply(grep("res",colnames(allsamplesgood@meta.data),value = TRUE),
function(x) length(unique(allsamplesgood@meta.data[,x])))## RNA_snn_res.1.5 RNA_snn_res.1.6 RNA_snn_res.1.7 RNA_snn_res.1.8 RNA_snn_res.1.9
## 19 21 22 23 23
## RNA_snn_res.2
## 26# change default identity
Idents(allsamplesgood) <- "RNA_snn_res.1.8"
# list cell number in each cluster for HC vs UV
table(Idents(allsamplesgood),allsamplesgood$orig.ident)##
## KO WT
## 0 620 250
## 1 384 397
## 2 382 335
## 3 68 540
## 4 323 284
## 5 348 163
## 6 321 173
## 7 280 182
## 8 86 343
## 9 138 263
## 10 83 313
## 11 193 197
## 12 183 177
## 13 157 172
## 14 165 125
## 15 96 194
## 16 66 198
## 17 48 189
## 18 50 73
## 19 45 65
## 20 45 45
## 21 38 48
## 22 11 14set.seed(20)
# Run UMAP with identified dimensions
allsamplesgood <- RunUMAP(allsamplesgood, dims = 1:20)
# Plot umap
DimPlot(object = allsamplesgood, pt.size=0.5, reduction = "umap", label = T)
# Color by WT vs KO
DimPlot(object = allsamplesgood, pt.size=0.5, reduction = "umap", group.by = "orig.ident")
# Plot umap
dimplot1<-DimPlot(object = allsamplesgood, pt.size=0.5, reduction = "umap", label = T)
# Color by WT vs KO
dimplot2<-DimPlot(object = allsamplesgood, pt.size=0.5, reduction = "umap", group.by = "orig.ident")
# Save umaps
ggsave(paste(outputpath,"LCMV-umap-clusters.pdf",sep=""), plot = dimplot1, width = 20, height = 15, units = "cm")
ggsave(paste(outputpath,"LCMV-umap-clusters-type.pdf",sep=""), plot = dimplot2, width = 20, height = 15, units = "cm")
# Custom list of genes to visualize
mygenes = c("Ly6c1", "Tbc1d4","Foxp3","Klrg1", "Irf8","Fcer1g","Gm44175","Cd8a","Cd8b1","Mcm3","Dut","Pclaf","Cd4","Klrb1c")
RidgePlot(allsamplesgood, features = mygenes)