Ecoli Pigs Single Cell Processing

Set working directory

knitr::opts_knit$set(root.dir = ".")

Load libraries

# load libraries
library(cowplot)  # plot_grid()
library(DoubletFinder)  # paramSweep()
library(dplyr)  # ungroup()
library(ggrepel) # geom_text_repel()
library(ggplot2)  # ggplot()
library(grid)  # grid.arrange()
library(gridExtra)  # grid.arrange()
library(harmony) # RunHarmony()
library(reshape2) # melt()
library(Seurat)  # Read10X_h5()
library(Seurat.utils) # RenameGenesSeurat()
library(stringr) # str_match()

Set variables and thresholds

pathToRef <- "/research/labs/neurology/fryer/projects/references/pig"
tissue <- "Brain"
control <- "Saline"
treatment <- "Ecoli"
control_color <- "gray"
treatment_color <- "green"
sample_colors <- c(control_color,treatment_color)
treatment_colors <- c(control_color, treatment_color)
myContrasts <- c("Ecoli - Control")
tool <- "cellranger"

nCount.min <- 500
nFeature.min <- 250
complexity.cutoff <- 0.85
mt.cutoff <- 15
hb.cutoff <- 5

Save functions

These functions with help simultaneously save plots as a png and pdf.

saveToPDF <- function(...) {
    d = dev.copy(pdf,...)
    dev.off(d)
}
saveToPNG <- function(...) {
    d = dev.copy(png,...)
    dev.off(d)
}

Setup seurat object

Using CellBender filtered output.

prefix <- "../../cellranger/"
suffix <- "/outs/filtered_feature_bc_matrix.h5"

if (tissue == "Brain" && file.exists("../../rObjects/brain_merged_h5.rds")) {
  pigs <- readRDS("../../rObjects/brain_merged_h5.rds")
} else if (tissue == "Brain") {
  # individual sample objects
  br1 <- CreateSeuratObject(Read10X_h5(paste0(prefix,"S1_BR",suffix)))
  br23 <- CreateSeuratObject(Read10X_h5(paste0(prefix,"E23_BR",suffix)))
  
  # merge objects
  pigs <- merge(x = br1, 
                y = c(br23), 
                add.cell.ids = c("1.Saline","23.Ecoli"), 
                project = paste0("Ecoli Pigs ", tissue, " Single Cell"))
  
  # cleanup and save
  remove(br1,br23)
  saveRDS(pigs, "../../rObjects/brain_merged_h5.rds")
  
} 
# add manual pig annotations (KK_gene_name)
genes <- read.csv(paste0(pathToRef, "/manual_pig_annotations.csv"))
new.names <- str_replace_all(rownames(pigs),
                             genes$ENSEMBL_ID,
                             genes$OUR_ID)

## Warning in stri_replace_all_regex(string, pattern,
## fix_replacement(replacement), : longer object length is not a multiple of
## shorter object length

pigs <- RenameGenesSeurat(pigs, new.names)

## [1] "Run this before integration. It only changes obj@assays$RNA@counts, @data and @scale.data."

# preview
pigs

## An object of class Seurat 
## 21303 features across 6941 samples within 1 assay 
## Active assay: RNA (21303 features, 0 variable features)

nCount_RNA = total number of transcripts (UMIs) in a single cell nFeature_RNA = number of unique genes (features)

# create sample column
barcodes <- colnames(pigs)
pattern <- "(.+)_[ACGT]+-(\\d+)"
sample <- str_match(barcodes, pattern)[,2]
table(sample)

## sample
## 1.Saline 23.Ecoli 
##     4925     2016

pigs$sample <- factor(sample, levels = c("1.Saline", "23.Ecoli"))
table(pigs$sample)  # check

## 
## 1.Saline 23.Ecoli 
##     4925     2016

Idents(pigs) <- pigs$sample

Add treatment column to metadata

# create treatment column
treat <- gsub("1.Saline", control, pigs$sample)
treat <- gsub("23.Ecoli", treatment, treat)

pigs$treatment <- factor(treat, levels = c("Saline","Ecoli"))
table(pigs$treatment)

## 
## Saline  Ecoli 
##   4925   2016

# cell.complexity
pigs$cell.complexity <- log10(pigs$nFeature_RNA) / log10(pigs$nCount_RNA)

# percent.mt
mt.genes <- c("ND1","ND2","COX1","COX2","ATP8","ATP6","COX3","ND3",
              "ND4L","ND4","ND5","ND6","CYTB")
pigs$percent.mt <- PercentageFeatureSet(pigs, features = mt.genes)

# percent.ribo
# ribosomal proteins begin with 'RPL' or 'RSL' for this annotation file
gene.names <- rownames(pigs)
ribo.genes <- gene.names[grep("^R[SP]L", gene.names)]
pigs$percent.ribo <- PercentageFeatureSet(pigs, features = ribo.genes)

# percent.hb
# percent.hb - hemoglobin proteins begin with 'HB' or 'HBP' for pig
hb.genes <- gene.names[grep("^HB[^(P)]", gene.names)]
pigs$percent.hb <- PercentageFeatureSet(pigs, features = hb.genes)

Pre-filtering QC

Number of cells

# Visualize the number of cell counts per sample
data <- as.data.frame(table(pigs$sample))
colnames(data) <- c("sample","frequency")

ncells1 <- ggplot(data, aes(x = sample, y = frequency, fill = sample)) + 
  geom_col() +
  theme_classic() +
  geom_text(aes(label = frequency), 
            position=position_dodge(width=0.9), 
            vjust=-0.25) +
  scale_fill_manual(values = sample_colors) + 
  scale_y_continuous(breaks = seq(0,5000, by = 1000), limits = c(0,5000)) +
  ggtitle("Raw: cells per sample") +
  theme(legend.position =  "none") + 
  theme(axis.text.x = element_text(angle = 45, hjust=1))
ncells1

Density plots

# set graphical parameter
par(mfrow = c(3,1))

# Visualize nCount_RNA
den1 <- ggplot(pigs@meta.data,
       aes(color = sample,
           x = nCount_RNA,
           fill = sample)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  scale_x_log10() +
  scale_color_manual(values = sample_colors) +
  scale_fill_manual(values = sample_colors) +
  xlab("nCount_RNA") +
  ylab("Density") +
  geom_vline(xintercept = nCount.min)

# Visualize percent.mt
den2 <- ggplot(pigs@meta.data,
       aes(color = sample,
           x = percent.mt,
           fill = sample)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  scale_x_continuous(n.breaks = 4) +
  geom_vline(xintercept = mt.cutoff) +
  scale_color_manual(values = sample_colors) +
  scale_fill_manual(values = sample_colors) +
  xlab("% Mitochondrial Genes") +
  ylab("Density")

# Visualize cell complexity
# Quality cells are usually above 0.85
den3 <- ggplot(pigs@meta.data,
       aes(color = sample,
           x = cell.complexity,
           fill = sample)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  scale_color_manual(values = sample_colors) +
  scale_fill_manual(values = sample_colors) +
  xlab("Cell Complexity (log10(nFeature/nCount))") +
  ylab("Density") +
  geom_vline(xintercept = complexity.cutoff)

# Arrange graphs in grid
plots1 <- list(den1,den2,den3)
layout1 <- rbind(c(1),c(2),c(3))
grid1 <- grid.arrange(grobs = plots1, layout_matrix = layout1)

Violin plots

# nFeature, nCount, and cell.complexity violins
v1 <- VlnPlot(pigs,
              features = c("nFeature_RNA", "nCount_RNA","cell.complexity"),
              ncol = 3,
              group.by = 'sample',
              cols = sample_colors,
              pt.size = 0)
v1

#  percent violins
v2 <- VlnPlot(pigs,
              features = c("percent.mt","percent.ribo","percent.hb"),
              ncol = 3,
              group.by = 'sample',
              cols = sample_colors,
              pt.size = 0)
v2

Scatter plots

s1 <- ggplot(
  pigs@meta.data,
  aes(x = nCount_RNA, y = nFeature_RNA, color = percent.mt)) + 
  geom_point() + 
  stat_smooth(method=lm) +
    scale_x_log10() +       
  scale_y_log10() + 
  theme_classic() +
  geom_vline(xintercept = nCount.min) + 
  geom_hline(yintercept = nFeature.min) + 
  facet_wrap(~sample) +
  scale_colour_gradient(low = "gray90", high = "black", limits =c(0,100))
  #geom_rect(aes(xmin=300, xmax=300, ymin=1000,
   #             ymax=3000), color="transparent", fill="orange", alpha=0.3)
s1

## `geom_smooth()` using formula 'y ~ x'

s2 <- FeatureScatter(pigs,
               feature1 = "nCount_RNA",
               feature2 = "percent.mt",
               group.by = 'sample',
               cols = sample_colors,
               shuffle = TRUE)
s2

Filtering

Cell-level filtering

We want to be careful filtering because removing things can easily lead to misinterpretation. For example, cells with high percent.mt could actually just be involved in respiratory processes.

We will filter based on 6 conditions:
– nCount_RNA > 500 – nFeature_RNA > 250 – cell.complexity > 0.85 – percent.mt < 1 – percent.hb < 5

And removing MT genes as to not alter down stream differential expression

# filter
pigs.filtered <- subset(pigs,
                        subset = (nCount_RNA > nCount.min) &
                          (nFeature_RNA > nFeature.min) &
                          (cell.complexity > complexity.cutoff) &
                          (percent.mt < mt.cutoff) &
                          (percent.hb < hb.cutoff))

# print cells removed
print(paste0(dim(pigs)[2] - dim(pigs.filtered)[2]," cells removed"))

## [1] "943 cells removed"

Gene-level filtering

Remove lowly expressed genes. We will keep genes that have at least 1 count in 10 cells.

# filter genes
counts <- GetAssayData(object = pigs.filtered, slot = "counts")
nonzero <- counts > 0  # produces logical
keep <- Matrix::rowSums(nonzero) >= 10  # sum the true/false
counts.filtered <- counts[keep,]  # keep certain genes

# overwrite pigs.filtered
pigs.filtered <- CreateSeuratObject(counts.filtered, 
                                    meta.data = pigs.filtered@meta.data)

# print features removed
print(paste0(dim(counts)[1] - dim(counts.filtered)[1], " features removed"))

## [1] "6930 features removed"

Doublet filtering

https://github.com/chris-mcginnis-ucsf/DoubletFinder\

Heterotypic - doublets derived from transcriptionally distinct cells. DoubletFinder works best on this type of doublet.

Homotopic - Transcriptionally similar cell doublets. DoubletFinder does not work as great on this type of doublet.

pANN - proportion of artificial nearest neighbors (pANN)

BCMVN - mean-variance normalized bimodality coefficient of pANN distributions produced during pN -pK parameter sweeps. The BCMVN may be used to identify the pK parameter.

Overview of steps:
A. Prepare each sample
B. pK Identification (no ground-truth) - defines the PC neighborhood size used to compute pANN
C. Homotypic Doublet Proportion Estimate - homotypic doublets may not be a problem depending on the type of analysis you are performing. If you have some doublets of the same type and their counts are normalized, they will generally represent the profile of single cells of the same type.
D. DoubletFinder
E. Visualize where the doublets are located

# split object by sample
pigs.split <- SplitObject(pigs.filtered, split.by = "sample") 

for (i in 1:length(pigs.split)) {
  # normalize and find PCs
  print(i)
  pig_sample <- NormalizeData(pigs.split[[i]])
  sampleID <- levels(droplevels(pig_sample@meta.data$sample))
  pig_sample <- FindVariableFeatures(pig_sample, selection.method = "vst", nfeatures = 2000)
  pig_sample <- ScaleData(pig_sample)
  pig_sample <- RunPCA(pig_sample)
  
  # get significant PCs
  stdv <- pig_sample[["pca"]]@stdev
  sum.stdv <- sum(pig_sample[["pca"]]@stdev)
  percent.stdv <- (stdv / sum.stdv) * 100
  cumulative <- cumsum(percent.stdv)
  co1 <- which(cumulative > 90 & percent.stdv < 5)[1]
  co2 <- sort(which((percent.stdv[1:length(percent.stdv) - 1] - 
                       percent.stdv[2:length(percent.stdv)]) > 0.1), 
              decreasing = T)[1] + 1
  min.pc <- min(co1, co2)
  min.pc
  
  # run umap
  pig_sample <- RunUMAP(pig_sample, dims = 1:min.pc, reduction = "pca")
  
  # cluster
  pig_sample <- FindNeighbors(object = pig_sample, dims = 1:min.pc)                           
  pig_sample <- FindClusters(object = pig_sample, resolution = 0.2)
  
  # Assign identity of clusters
  Idents(object = pig_sample) <- "RNA_snn_res.0.2"
  d1 <- DimPlot(pig_sample,
        reduction = "umap",
        label = TRUE,
        label.size = 6)
  path <- paste0("../../results/doubletFinder/",treatment,"_",tolower(tissue),
               "_doubletFinder_UMAP_res0.02_",sampleID)
  pdf(paste0(path, ".pdf"), width = 5, height = 4)
  print(d1)
  dev.off()
  
  # number of cells in each cluster
  n_cells <- FetchData(pig_sample, vars = c("ident")) %>% dplyr::count(ident) %>%tidyr::spread(ident, n)
  
  ## pK Identification (no ground-truth) 
  sweep.res.list <- paramSweep_v3(pig_sample, PCs = 1:min.pc, sct = FALSE)
  sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)
  bcmvn <- find.pK(sweep.stats)
  
  # Optimal pK for any scRNA-seq data can be manually discerned as maxima in BCmvn distributions
  bcmvn_max <- bcmvn[which.max(bcmvn$BCmetric),]
  pK_value <- bcmvn_max$pK
  pK_value <- as.numeric(levels(pK_value))[pK_value]
  
  # Homotypic Doublet Proportion Estimate 
  annotations <- pig_sample@meta.data$seurat_clusters
  homotypic.prop <- modelHomotypic(annotations) 
  nExp_poi <- round(pK_value*nrow(pig_sample@meta.data))
  nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))
  # Run DoubletFinder with varying classification  
  pig_sample <- doubletFinder_v3(pig_sample, PCs = 1:min.pc, 
                pN = 0.25, pK = pK_value, nExp = nExp_poi.adj, 
                reuse.pANN = FALSE, sct = FALSE)
  
  # set DF class for calling doublets
  DF_class <- pig_sample@meta.data[, grep("DF.classifications",colnames(pig_sample@meta.data)),]
  DF_class[which(DF_class == "Doublet")] <- "Doublet"
  table(DF_class)
  
  # table showing the number of doublets and singlets
  write.table(table(DF_class), paste0("../../results/doubletFinder/",treatment,"_",tolower(tissue),
               "_doubletFinder_table_",sampleID), sep = "\t", 
              row.names = FALSE, quote = FALSE)
  pig_sample@meta.data[,"CellTypes_DF"] <- DF_class
  
  # plot
  d2 <- DimPlot(pig_sample, group.by="CellTypes_DF", reduction="umap",
          order=c("Coll.Duct.TC","Doublet"), 
          cols=c("#66C2A5","black"))
  path <- paste0("../../results/doubletFinder/",treatment,"_",tolower(tissue),
               "_doubletFinder_UMAP_",sampleID)
  pdf(paste0(path, ".pdf"), width = 5,height = 4)
  print(d2)
  dev.off()
  
  # plot
  f1 <- FeaturePlot(pig_sample, 
            reduction = "umap", 
            features = c("nFeature_RNA", "nCount_RNA", 
                         "cell.complexity", "percent.mt"),
            pt.size = 0.4, 
            order = TRUE,
            label = TRUE)
  path <- paste0("../../results/doubletFinder/",treatment,"_",tolower(tissue),
               "_doubletFinder_FeaturePlot_",sampleID)
  pdf(paste0(path, ".pdf"), width = 7, height = 7)
  print(f1)
  dev.off()
  
  #only keep singlets
  pig_sample_singlets <- subset(pig_sample, subset = CellTypes_DF == "Singlet")
  
  # inspect
  d3 <- DimPlot(pig_sample_singlets, group.by="CellTypes_DF", reduction="umap",
          order=c("Coll.Duct.TC","Doublet"),
          cols=c("#66C2A5","black"))
  path <- paste0("../../results/doubletFinder/",treatment,"_",tolower(tissue),
               "_doubletFinder_UMAP_singlets_",sampleID)
  pdf(paste0(path, ".pdf"), width = 5, height = 4)
  print(d3)
  dev.off()
  
  # number of cells in each cluster per and post removing doublets
  n_cells_singlets <- FetchData(pig_sample_singlets, vars = c("ident")) %>% dplyr::count(ident) %>% tidyr::spread(ident, n)
  n_cells_singlets
  ncells_per_cluster <- rbind(n_cells, n_cells_singlets)
  row.names(ncells_per_cluster) <- c("Doublets and singlets", "Singlets only")
  ncells_per_cluster
  difference <- diff(as.matrix(ncells_per_cluster))
  difference <- as.data.frame(difference)
  row.names(difference) <- c("difference")
  cbind(difference, ncells_per_cluster)
  write.table(ncells_per_cluster, paste0(
    "../../results/doubletFinder/",treatment,"_",tolower(tissue),
    "_doubletFinder_table_ncells_per_cluster",sampleID, ".txt"), sep = "\t", 
    row.names = FALSE, quote = FALSE)
  
  # plot the number of cells in each cluster per and post doubletFinder
  ncell_matrix <- as.matrix(ncells_per_cluster)
  ncells_melt <- melt(ncell_matrix)
  colnames(ncells_melt) <- c("doublet type","cluster","number of cells")
  ncell_max <- ncells_melt[which.max(ncells_melt$`number of cells`),]
  ncell_max_value <- ncell_max$`number of cells`
  cellmax <- ncell_max_value + 800 # so that the figure doesn't cut off the text
  b1 <- ggplot(ncells_melt, aes(x = factor(cluster), y = `number of cells`,
                          fill = `doublet type`)) + 
    geom_bar(stat="identity", colour="black", width=1, position = position_dodge(width=0.8)) +
    geom_text(aes(label = `number of cells`), 
              position=position_dodge(width=0.9), vjust=-0.25, angle = 45, hjust=-.01) + 
    theme_classic() + scale_fill_manual(values = c("gray", "#66C2A5")) +
    ggtitle("Number of cells per cluster") +  xlab("cluster") +
    theme(axis.text.x = element_text(angle = 45, hjust=1)) +
    scale_y_continuous(limits = c(0,cellmax))
  path <- paste0("../../results/doubletFinder/",treatment,"_",tolower(tissue),
               "_doubletFinder_barplot_ncells_per_cluster",sampleID)
  pdf(paste0(path, ".pdf"), width = 7,height = 5)
  print(b1)
  dev.off()
  f2 <- FeaturePlot(pig_sample_singlets, 
            reduction = "umap", 
            features = c("nFeature_RNA", "nCount_RNA", 
                         "cell.complexity", "percent.mt"),
            pt.size = 0.4, 
            order = TRUE,
            label = TRUE)
  path <- paste0("../../results/doubletFinder/",treatment,"_",tolower(tissue),
               "_doubletFinder_FeaturePlot_singlets",sampleID)
  pdf(paste0(path, ".pdf"), width = 7,height = 7)
  print(f2)
  dev.off()
  
  # put the pigs together again
  pigs.split[[i]] <- pig_sample_singlets
}

## [1] 1
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 4319
## Number of edges: 147414
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9081
## Number of communities: 10
## Elapsed time: 0 seconds

## [1] "Creating artificial doublets for pN = 5%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 10%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 15%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 20%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 25%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 30%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.005..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."

## NULL
## [1] "Creating 1440 artificial doublets..."
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Computing pANN..."
## [1] "Classifying doublets.."

## [1] 2
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1679
## Number of edges: 55343
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8865
## Number of communities: 6
## Elapsed time: 0 seconds

## [1] "Creating artificial doublets for pN = 5%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 10%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 15%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 20%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 25%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."
## [1] "Creating artificial doublets for pN = 30%"
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Defining neighborhoods..."
## [1] "Computing pANN across all pK..."
## [1] "pK = 0.01..."
## [1] "pK = 0.02..."
## [1] "pK = 0.03..."
## [1] "pK = 0.04..."
## [1] "pK = 0.05..."
## [1] "pK = 0.06..."
## [1] "pK = 0.07..."
## [1] "pK = 0.08..."
## [1] "pK = 0.09..."
## [1] "pK = 0.1..."
## [1] "pK = 0.11..."
## [1] "pK = 0.12..."
## [1] "pK = 0.13..."
## [1] "pK = 0.14..."
## [1] "pK = 0.15..."
## [1] "pK = 0.16..."
## [1] "pK = 0.17..."
## [1] "pK = 0.18..."
## [1] "pK = 0.19..."
## [1] "pK = 0.2..."
## [1] "pK = 0.21..."
## [1] "pK = 0.22..."
## [1] "pK = 0.23..."
## [1] "pK = 0.24..."
## [1] "pK = 0.25..."
## [1] "pK = 0.26..."
## [1] "pK = 0.27..."
## [1] "pK = 0.28..."
## [1] "pK = 0.29..."
## [1] "pK = 0.3..."

## NULL
## [1] "Creating 560 artificial doublets..."
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Computing pANN..."
## [1] "Classifying doublets.."

# converge pigs.split
pigs.singlets <- merge(x = pigs.split[[1]],
                       y = c(pigs.split[[2]]),
                       project = paste0("Ecoli Pigs ", tissue, " Single Cell"))

# print how many cells removed
print(paste0(dim(pigs.filtered)[2] - dim(pigs.singlets)[2]," cells removed"))

## [1] "347 cells removed"

# how many removed if we had an upper nCount and nFeature
pigs.upper <- subset(pigs.singlets,
                     subset = (nCount_RNA < 10000) & (nFeature_RNA < 5000))
print(paste0(dim(pigs.filtered)[2] - dim(pigs.upper)[2],
             " cells would have been removed if upper bound applied"))

## [1] "508 cells would have been removed if upper bound applied"

# overwrite pigs.filtered
pigs.filtered <- pigs.singlets

# reset levels
pigs.filtered$treatment <- factor(pigs.filtered$treatment,
                                  levels = c("Saline","Ecoli"))
pigs.filtered$sample <- factor(pigs.filtered$sample,
                               levels = c("1.Saline", "23.Ecoli"))

# cleanup
remove(pigs.singlets, pigs.upper, pigs.split, pig_sample_singlets, pig_sample)
remove(n_cells,n_cells_singlets,ncell_matrix,ncell_max,ncells_per_cluster,ncells_melt)
remove(sweep.res.list, sweep.stats,bcmvn,bcmvn_max,difference)
remove(d1,d2,d3,f1,f2)
remove(counts,counts.filtered, nonzero)

Mitochondrial gene filtering

# remove mt.genes
counts <- GetAssayData(object = pigs.filtered, slot = "counts")
keep <- !rownames(counts) %in% mt.genes # false when mt.gene
counts.filtered <- counts[keep,]

# overwrite pigs.filtered
pigs.filtered <- CreateSeuratObject(counts.filtered, 
                                    meta.data = pigs.filtered@meta.data)

# print features removed
print(paste0(dim(counts)[1] - dim(counts.filtered)[1], " features removed"))

## [1] "13 features removed"

Histogram

# User params
goi <- "KK-MALAT1"
threshold <- 5

# Subset data
log2.threshold <- log2(threshold + 0.01)
counts.df <- FetchData(pigs.filtered, vars = "goi")
colnames(counts.df) <- "counts"
log2.counts.df <- log2(counts.df + 0.01)

# Histogram
title <- paste0("Ecoli Brain Nuclei: ", goi, "\nnCount_RNA > ", threshold)
hist1 <- ggplot(counts.df, aes(x = counts)) + 
  geom_histogram(bins = 100, fill = "gray", color = "black") + 
  labs(title = title, x=NULL, y=NULL) +
  xlab(paste0(goi, " nCount_RNA")) + ylab("# of Samples") + theme_bw() +
  geom_vline(xintercept = threshold, col = "blue") +
  annotate("rect", 
           xmin = -Inf,
           xmax = threshold, 
           ymin = 0, 
           ymax=Inf, 
           alpha=0.2, 
           fill="chocolate4") +
  annotate("rect", 
           xmin = threshold,
           xmax = Inf, 
           ymin = 0, 
           ymax=Inf, 
           alpha=0.2, 
           fill="deepskyblue")

# Histogram log transformed
hist2 <- ggplot(log2.counts.df, aes(x = counts)) + 
  geom_histogram(bins = 100, fill = "gray", color = "black") + 
  labs(title = title, x=NULL, y=NULL) +
  xlab(paste0("Log2(",goi, " nCount_RNA)")) + ylab("# of Samples") + theme_bw() +
  geom_vline(xintercept = log2.threshold, col = "blue") +
  annotate("rect", 
           xmin = -Inf,
           xmax = log2.threshold, 
           ymin = 0, 
           ymax=Inf, 
           alpha=0.2, 
           fill="chocolate4") +
  annotate("rect", 
           xmin = log2.threshold,
           xmax = Inf, 
           ymin = 0, 
           ymax=Inf, 
           alpha=0.2, 
           fill="deepskyblue")

# plot
plots1 <- list(hist1,hist2)
layout1 <- rbind(c(1),c(2))
grid1 <- grid.arrange(grobs = plots1, layout_matrix = layout1)

# number removed
table(counts.df$counts > threshold)

Number of cells

# Visualize the number of cell counts per sample
data <- as.data.frame(table(pigs.filtered$sample))
colnames(data) <- c("sample","frequency")

ncells2 <- ggplot(data, aes(x = sample, y = frequency, fill = sample)) + 
  geom_col() +
  theme_classic() +
  geom_text(aes(label = frequency), 
            position=position_dodge(width=0.9), 
            vjust=-0.25) +
  scale_fill_manual(values = sample_colors) + 
  scale_y_continuous(breaks = seq(0,5000, by = 1000), limits = c(0,5000)) +
  ggtitle("Filtered: cells per sample") +
  theme(legend.position =  "none") + 
  theme(axis.text.x = element_text(angle = 45, hjust=1))

# Arrange graphs in grid
plots2 <- list(ncells1,ncells2)
layout2 <- cbind(c(1),c(2))
grid2 <- grid.arrange(grobs = plots2, layout_matrix = layout2)

Density plots

# set graphical parameter
par(mfrow = c(3,1))

# Visualize the number of counts per cell
den4 <- ggplot(pigs.filtered@meta.data,
       aes(color = sample,
           x = nCount_RNA,
           fill = sample)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  scale_x_log10() +
  scale_color_manual(values = sample_colors) +
  scale_fill_manual(values = sample_colors) +
  xlab("nCount_RNA") +
  ylab("Density") +
  geom_vline(xintercept = nCount.min)

# Visualize percent.mt
den5 <- ggplot(pigs.filtered@meta.data,
       aes(color = sample,
           x = percent.mt,
           fill = sample)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  scale_x_log10() +
  scale_color_manual(values = sample_colors) +
  scale_fill_manual(values = sample_colors) +
  xlab("% Mitochondrial Genes") +
  ylab("Density") +
  geom_vline(xintercept = mt.cutoff)

# Visualize cell complexity
# Quality cells are usually above 0.80
den6 <- ggplot(pigs.filtered@meta.data,
       aes(color = sample,
           x = cell.complexity,
           fill = sample)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  scale_x_log10() +
  scale_color_manual(values = sample_colors) +
  scale_fill_manual(values = sample_colors) +
  xlab("Cell Complexity (log10(nFeature/nCount))") +
  ylab("Density") +
  geom_vline(xintercept = complexity.cutoff)

# Arrange graphs in grid
plots3 <- list(den1,den2,den3,den4,den5,den6)
layout3 <- rbind(c(1,4),c(2,5),c(3,6))
grid3 <- grid.arrange(grobs = plots3, layout_matrix = layout3)

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Removed 9 rows containing non-finite values (stat_density).

Violin plots

# nFeature, nCount, and cell.complexity violins
v3 <- VlnPlot(pigs.filtered,
              features = c("nFeature_RNA", "nCount_RNA","cell.complexity"),
              ncol = 3,
              group.by = 'sample',
              cols = sample_colors,
              pt.size = 0)
v3

#  percent violins
v4 <- VlnPlot(pigs.filtered,
              features = c("percent.mt","percent.ribo","percent.hb"),
              ncol = 3,
              group.by = 'sample',
              cols = sample_colors,
              pt.size = 0)
v4

Scatter plots

s3 <- ggplot(
  pigs.filtered@meta.data,
  aes(x = nCount_RNA, y = nFeature_RNA, color = percent.mt)) + 
  geom_point() + 
  stat_smooth(method=lm) +
    scale_x_log10() +       
  scale_y_log10() + 
  theme_classic() +
  geom_vline(xintercept = nCount.min) + 
  geom_hline(yintercept = nFeature.min) + 
  facet_wrap(~sample) +
  scale_colour_gradient(low = "gray90", high = "black", limits =c(0,100))
  #geom_rect(aes(xmin=300, xmax=300, ymin=1000,
   #             ymax=3000), color="transparent", fill="orange", alpha=0.3)
s3

## `geom_smooth()` using formula 'y ~ x'

s4 <- FeatureScatter(pigs.filtered,
               feature1 = "nCount_RNA",
               feature2 = "percent.mt",
               group.by = 'sample',
               cols = sample_colors,
               shuffle = TRUE)
s4

Box plot

# Visualize the distribution of genes detected per cell via boxplot
b1 <- ggplot(pigs.filtered@meta.data,
       aes(x = sample, 
           y = log10(nFeature_RNA), 
           fill=sample)) + 
  geom_boxplot() + 
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  theme(plot.title = element_text(hjust = 0.5, face="bold")) +
  ggtitle("Unique Genes / Cell / Sample") +
  scale_color_manual(values = sample_colors) +
  scale_fill_manual(values = sample_colors) +
  xlab("Sample")
b1

Top transcripts

df <- data.frame(row.names = rownames(pigs.filtered))
df$rsum <- rowSums(x = pigs.filtered, slot = "counts")
df$gene_name <- rownames(df)
df <- df[order(df$rsum,decreasing = TRUE),]
head(df, 10)

##                      rsum          gene_name
## CST3               845555               CST3
## TMSB4X             348715             TMSB4X
## B2M                206036                B2M
## ENSSSCG00000001229 122442 ENSSSCG00000001229
## C1QB               119084               C1QB
## TMSB10             114200             TMSB10
## APOE               108719               APOE
## ENSSSCG00000036618  94103 ENSSSCG00000036618
## ENSSSCG00000004970  92692 ENSSSCG00000004970
## TPT1                89121               TPT1

Next steps

For something to be informative, it needs to exhibit variation, but not all variation is informative. The goal of our clustering analysis is to keep the major sources of variation in our dataset that should define our cell types, while restricting the variation due to uninteresting sources of variation (sequencing depth, cell cycle differences, mitochondrial expression, batch effects, etc.). Then, to determine the cell types present, we will perform a clustering analysis using the most variable genes to define the major sources of variation in the dataset.

1. Explain sources of unwanted variation: cell cycle, percent.mt (cell stress)
2. Normalizing and regressing out sources of unwanted variation: sctransform
3. Integration (OPTIONAL STEP): recommended when comparing groups, uses shared sources of high variation to identify shared sub-populations across conditions, we visualize after integration to ensure it’s good before clustering
   • Stuart and Butler et al. (2018)
4. Clustering cells: cluster cells based on similarity of their gene expression profiles, cells are assigned to a cluster based on their PCA scores based on expression of integrated highly variable genes
5. Cluster quality evaluation: see what clusters aren’t influenced by sources of uninteresting variation, check major PCs are driving clusters, explore cell identities by looking a known markers

Cell cycle

The most common biological data correction is to remove the effects of the cell cycle on the transcriptome. This data correction can be performed by a simple linear regression against a cell cycle score.

Check cell cycle phase BEFORE doing sctransform. Counts to need to be comparable between cells and each cell has a different number for nCount_RNA.

# summary of counts per cell
summary(pigs.filtered@meta.data$nCount_RNA)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     501    2690    4002    4305    5458   19514

Use the NormalizeData() function with the argument LogNormalze to account for sequencing depth. nCount_RNA for each gene is divided by the total nCount_RNA for that cell. This is done for all cells. This number is then multiplied by the scale.factor so we don’t have to deal with a tiny number. This number is then natural-log transformed using log1p. log1p is the natural logarithm (base e) of 1 + count. The 1 will prevent taking the log of 0.

# EXAMPLE
# TUBB1 gene expression in cell 1
MALAT1_cell1_count <- FetchData(pigs.filtered, vars = "KK-MALAT1", cells = 1)
MALAT1_cell1_count <- MALAT1_cell1_count$`KK-MALAT1`
# get total counts in cell 1
cell1_count <- pigs.filtered$nCount_RNA[1]
# divide, scale, and log1p transform
log1p((MALAT1_cell1_count/cell1_count)*10000)

pigs.phase <- NormalizeData(pigs.filtered,
                            scale.factor = 10000, # default
                            normalization.method = "LogNormalize" # default
                            )

# check
#FetchData(pigs.phase, vars = "KK-MALAT1", cells = 1)

Give each cell a score based on expression of G1, G2/M, and S phase markers. A list of markers is provided for humans. Since the pig genome is poorly annotated we will use the human list. We will use the CellCycleScoring() function in seurat.

Below is a resource for acquiring cell markers for other organisms https://hbctraining.github.io/scRNA-seq_online/lessons/cell_cycle_scoring.html

G1 ~10 hrs S ~5-6 hrs G2 ~3-4 hrs M ~2 hrs

G1 (10 hrs) > G2/M (5-6 hrs) = S (5-6 hrs)

If the score is negative for both S.Score and G2M.Score the phase is G1. Otherwise the the greatest positive value between S.Score and G2M.Score determines the phase.

g2m_genes <- c("NCAPD2","ANLN","KK-TACC3","HMMR","GTSE1","NDC80","KK-AURKA",
               "TPX2","BIRC5","G2E3","CBX5","RANGAP1","CTCF","CDCA3","TTK",
               "SMC4","ECT2","CENPA","CDC20","NEK2","CENPF","TMPO","HJURP",
               "CKS2","DLGAP5","PIMREG","TOP2A","PSRC1","CDCA8","CKAP2",
               "KK-NUSAP1","KIF23","KIF11","KIF20B","CENPE","GAS2L3","KIF2C",
               "NUF2","KK-ANP32E","LBR","MKI67","CCNB2","CDC25C","HMGB2",
               "CKAP2L","BUB1","CDK1","CKS1B","UBE2C","CKAP5","AURKB","CDCA2",
               "TUBB4B","JPT1")
s_genes <- c("KK-UBR7","RFC2","RAD51","MCM2","TIPIN","MCM6","UNG","POLD3",
             "WDR76","CLSPN","CDC45","CDC6","MSH2","MCM5","POLA1","MCM4",
             "RAD51AP1","GMNN","RPA2","CASP8AP2","HELLS","E2F8","GINS2","PCNA",
             "NASP","BRIP1","DSCC1","DTL","CDCA7","CENPU","ATAD2","CHAF1B",
             "USP1","KK-SLBP","RRM1","FEN1","KK-RRM2","EXO1","CCNE2","TYMS",
             "BLM","KK-PRIM1","UHRF1")

# score cells for cell cycle
pigs.phase <- CellCycleScoring(pigs.phase,
                               g2m.features = g2m_genes,
                               s.features = s_genes,
                               set.ident = TRUE)

## Warning: The following features are not present in the object: KK-UBR7, CHAF1B,
## KK-SLBP, KK-RRM2, KK-PRIM1, not searching for symbol synonyms

## Warning: The following features are not present in the object: KK-TACC3, KK-
## AURKA, KK-NUSAP1, KK-ANP32E, MKI67, CKS1B, not searching for symbol synonyms

Bar graph

cellcyclecount_barplot <- 
  as_tibble(pigs.phase[[]]) %>%
  ggplot(aes(Phase, fill = Phase)) + geom_bar()
cellcyclecount_barplot

Pie chart

# pie point 
cellcyclecount_piepoint <- 
  as_tibble(pigs.phase[[]]) %>%
  ggplot(aes(x=S.Score, y=G2M.Score, color=Phase)) + 
  geom_point()
cellcyclecount_piepoint

Top variable genes

# Identify the most variable genes
pigs.phase <- FindVariableFeatures(pigs.phase,
                                   selection.method = "vst",  # default vst
                                   nfeatures = 2000,  # default 2000
                                   verbose = FALSE)

# view top variable genes
top40 <- head(VariableFeatures(pigs.phase), 40)
top40

##  [1] "GNLY"               "CXCL9"              "ENSSSCG00000026043"
##  [4] "CXCL11"             "CD163"              "CCL5"              
##  [7] "HPS5"               "CXCL10"             "S100A12"           
## [10] "MRC1"               "CCL4"               "SDS"               
## [13] "TOP2A"              "MT3"                "HMGB2"             
## [16] "ENSSSCG00000001456" "CAMK2A"             "SPP1"              
## [19] "ENSSSCG00000013788" "IL1R2"              "S100A4"            
## [22] "ACSL4"              "CENPF"              "PLA2G2D"           
## [25] "ENSSSCG00000005638" "CCL2"               "TIMP1"             
## [28] "STMN1"              "XCL1"               "VIM"               
## [31] "CPSF6"              "FOS"                "KIF11"             
## [34] "MMP8"               "APOA1"              "CLU"               
## [37] "ENSSSCG00000000261" "ISG15"              "HLA-DRA"           
## [40] "STEAP4"

# plot variable features with labels
VarFeatPlot <- VariableFeaturePlot(pigs.phase, cols = c("gray47","red"))
VarFeatPlotLabel <- LabelPoints(plot = VarFeatPlot, 
                    points = top40, repel = TRUE, fontface="italic", 
                    xnudge = 0, ynudge = 0, max.overlaps = 12)
VarFeatPlotLabel

Mean-variance

# The variability information can be accessed using the HVFInfo method. 
# The names of the variable features can be accessed with VariableFeatures().
variance.data <- as_tibble(HVFInfo(pigs.phase),rownames = "Gene")
variance.data <- variance.data %>% mutate(hypervariable=Gene %in% VariableFeatures(pigs.phase))

# We can plot out a graph of the variance vs mean and highlight the selected genes 
# this way, we can see whether we think we’re likely to capture what we need.
subset_data <- subset(variance.data, hypervariable == TRUE)
varGeneslog <- variance.data %>% 
  ggplot(aes(log(mean),log(variance),color=hypervariable)) + 
  geom_point() + 
  scale_color_manual(values=c("black","red")) +  geom_text_repel(
    data = subset_data, max.overlaps = 20,
    aes(
      x = log(mean),
      y = log(variance),
      label = Gene, 
      fontface="italic",),segment.alpha = 1,size = 4) +
  theme(legend.position="bottom")
varGeneslog

PCA

See if the cell cycle is a major source of variation using PCA. Choose the most variable gene features (we have already done), then e the data. We scale the data because highly expressed genes exhibit the highest amount of variation and we don’t want our ‘highly variable genes’ only to reflect high expression.

vst: First, fits a line to the relationship of log(variance) and log(mean) using local polynomial regression (loess). Then, feature values are standardized using the observed mean and expected variance (given by the fitted line). Feature variance is then calculated on the standardized values after clipping to a maximum (see clip.max parameter).

The ScaleData() function in Seurat will adjust gene expressions so that the mean expression in each cell is 0. It will also scale each gene to give a variance of 1 for each cell.

# Scale the counts
pigs.phase <- ScaleData(pigs.phase)

## Centering and scaling data matrix

pigs.phase@assays

## $RNA
## Assay data with 14360 features for 5651 cells
## Top 10 variable features:
##  GNLY, CXCL9, ENSSSCG00000026043, CXCL11, CD163, CCL5, HPS5, CXCL10,
## S100A12, MRC1

pigs.phase.pca <- RunPCA(pigs.phase, features = c(s_genes, g2m_genes))

## Warning in PrepDR(object = object, features = features, verbose = verbose): The
## following 38 features requested have not been scaled (running reduction without
## them): KK-UBR7, RFC2, RAD51, TIPIN, WDR76, CDC45, MSH2, MCM5, POLA1, RAD51AP1,
## GMNN, RPA2, GINS2, NASP, BRIP1, DSCC1, DTL, CHAF1B, KK-SLBP, RRM1, FEN1, KK-
## RRM2, CCNE2, TYMS, KK-PRIM1, KK-TACC3, KK-AURKA, CBX5, RANGAP1, CTCF, TMPO, KK-
## NUSAP1, KK-ANP32E, LBR, MKI67, CDC25C, CKS1B, JPT1

## Warning in irlba(A = t(x = object), nv = npcs, ...): You're computing too large
## a percentage of total singular values, use a standard svd instead.

## Warning: Requested number is larger than the number of available items (59).
## Setting to 59.

## Warning: Requested number is larger than the number of available items (59).
## Setting to 59.

## Warning: Requested number is larger than the number of available items (59).
## Setting to 59.

## Warning: Requested number is larger than the number of available items (59).
## Setting to 59.

## Warning: Requested number is larger than the number of available items (59).
## Setting to 59.

## PC_ 1 
## Positive:  PSRC1, CDCA7, CASP8AP2, UNG, BLM, POLD3, G2E3, MCM6, HELLS, EXO1 
##     GAS2L3, USP1, MCM4, CKAP5, CDC6, MCM2, PCNA, ANLN, PIMREG, KIF20B 
##     HJURP, CENPU, TUBB4B, ECT2, AURKB, CDCA2, E2F8, ATAD2, CLSPN, CENPA 
## Negative:  KIF11, HMMR, TPX2, TOP2A, NDC80, CENPE, CKAP2, CENPF, DLGAP5, UBE2C 
##     KIF23, UHRF1, BIRC5, CDCA8, CCNB2, NUF2, SMC4, CDCA3, CKS2, CDK1 
##     NCAPD2, BUB1, CDC20, CKAP2L, HMGB2, KIF2C, TTK, GTSE1, NEK2 
## PC_ 2 
## Positive:  MCM4, MCM2, CENPU, UNG, PCNA, CLSPN, ATAD2, UHRF1, HELLS, MCM6 
##     CDC6, POLD3, CASP8AP2, BLM, EXO1, USP1, CDCA7, E2F8, HMGB2, KIF11 
##     SMC4, PSRC1, CKAP2L, CKAP5, TOP2A, TPX2, ANLN, G2E3, CKAP2, NDC80 
## Negative:  CENPA, PIMREG, UBE2C, NUF2, CCNB2, HJURP, BUB1, CDCA3, GTSE1, KIF2C 
##     CENPF, ECT2, DLGAP5, KIF23, CDK1, NEK2, BIRC5, CDC20, AURKB, TTK 
##     GAS2L3, KIF20B, NCAPD2, CDCA8, TUBB4B, CKS2, CDCA2, CENPE, HMMR 
## PC_ 3 
## Positive:  G2E3, CASP8AP2, MCM6, USP1, KIF20B, CKAP5, HELLS, BLM, ANLN, ECT2 
##     NCAPD2, POLD3, NUF2, CKS2, CKAP2L, NEK2, PCNA, HJURP, PSRC1, CDCA2 
##     PIMREG, KIF2C, BUB1, TUBB4B, KIF23, CENPA, CDK1, GAS2L3, CDCA3, GTSE1 
## Negative:  HMGB2, MCM4, CLSPN, UNG, AURKB, SMC4, MCM2, UHRF1, CENPU, CDC6 
##     TOP2A, HMMR, CKAP2, KIF11, BIRC5, EXO1, CDCA8, TPX2, NDC80, UBE2C 
##     E2F8, CENPF, CDC20, DLGAP5, CENPE, CCNB2, TTK, ATAD2, CDCA7 
## PC_ 4 
## Positive:  USP1, BLM, MCM6, ANLN, CKAP5, MCM4, ECT2, CKAP2L, PSRC1, ATAD2 
##     CKAP2, UHRF1, MCM2, E2F8, CDC6, CENPE, CDCA2, TPX2, KIF11, NDC80 
##     KIF23, CDCA8, EXO1, CDCA7, TOP2A, NEK2, HMMR, CENPU, GAS2L3, GTSE1 
## Negative:  POLD3, HELLS, G2E3, KIF20B, SMC4, CASP8AP2, CLSPN, HMGB2, AURKB, CKS2 
##     CDK1, NCAPD2, TUBB4B, PCNA, BIRC5, UNG, NUF2, UBE2C, CCNB2, HJURP 
##     CENPF, PIMREG, BUB1, TTK, KIF2C, CDCA3, CENPA, CDC20, DLGAP5 
## PC_ 5 
## Positive:  CKAP5, TUBB4B, CASP8AP2, USP1, CKS2, UNG, POLD3, PCNA, CDK1, AURKB 
##     HMGB2, BIRC5, SMC4, BLM, KIF20B, CENPA, TTK, MCM4, CDCA3, CDCA2 
##     CDCA8, PIMREG, KIF23, UBE2C, CDC20, TPX2, KIF11, CCNB2, HMMR, CLSPN 
## Negative:  HELLS, ECT2, NCAPD2, ANLN, G2E3, MCM6, ATAD2, UHRF1, PSRC1, CKAP2L 
##     CENPF, GTSE1, E2F8, CENPE, CDCA7, NUF2, NEK2, CENPU, DLGAP5, EXO1 
##     CDC6, NDC80, MCM2, CKAP2, TOP2A, HJURP, GAS2L3, BUB1, KIF2C

DimPlot(pigs.phase.pca)

If the plots for each phase look very similar to each other, do not regress out variation due to cell cycle. You can plot PC1 vs PC2 before and after regression to see how effective it was. G1 (10 hrs) > G2/M (5-6 hrs) = S (5-6 hrs)

# Perform PCA
pigs.phase <- RunPCA(pigs.phase)

# Plot the PCA colored by cell cycle phase
cycle.pca <- DimPlot(pigs.phase,
        reduction = "pca",
        group.by= "Phase",
        split.by = "Phase")
cycle.pca

SCTransform

Now, we can use the SCTransform method as a more accurate method of normalizing, estimating the variance of the raw filtered data, and identifying the most variable genes. Variation in sequencing depth (total nCount_RNA per cell) is normalized using a regularized negative binomial model. ??Variance is also adjusted based on pooling information across genes with similar abundances??

Sctransform automatically accounts for cellular sequencing depth by regressing out sequencing depth (nUMIs). However, if there are other sources of uninteresting variation identified in the data during the exploration steps we can also include these. We observed little to no effect due to cell cycle phase and so we chose not to regress this out of our data. We observed some effect of mitochondrial expression and so we choose to regress this out from the data.

Since we have four samples in our dataset (from two conditions), we want to keep them as separate objects and transform them as that is what is required for integration. We will first split the cells in seurat.phase object by sample.

# Split seurat object by timepoint to perform SCT on all samples
pigs.split <- SplitObject(pigs.phase, split.by = "sample")

Now we will use a ‘for loop’ to run the SCTransform() on each sample, and regress out mitochondrial expression by specifying in the vars.to.regress argument of the SCTransform() function.

Before we run this for loop, we know that the output can generate large R objects/variables in terms of memory. If we have a large dataset, then we might need to adjust the limit for allowable object sizes within R (Default is 500 * 1024 ^ 2 = 500 Mb) using the following code:

options(future.globals.maxSize = 4000 * 1024^5)

for (i in 1:length(pigs.split)) {
  print(paste0("Sample ", i))
  pigs.split[[i]] <- SCTransform(pigs.split[[i]]
                               #  vars.to.regress = c("percent.mt") 
                                 )
}

## [1] "Sample 1"
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   4%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |======================                                                |  32%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |============================                                          |  39%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |==========================================                            |  61%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |================================================                      |  68%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |====================================================================  |  96%
  |                                                                            
  |======================================================================| 100%
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   4%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |======================                                                |  32%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |============================                                          |  39%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |==========================================                            |  61%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |================================================                      |  68%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |====================================================================  |  96%
  |                                                                            
  |======================================================================| 100%
## [1] "Sample 2"
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |======================================================================| 100%
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |======================================================================| 100%
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |======================================================================| 100%

NOTE: By default, after normalizing, adjusting the variance, and regressing out uninteresting sources of variation, SCTransform will rank the genes by residual variance and output the 3000 most variant genes. If the dataset has larger cell numbers, then it may be beneficial to adjust this parameter higher using the variable.features.n argument.

Note, the last line of output specifies “Set default assay to SCT”. We can view the different assays that we have stored in our seurat object.

A thread about whether or not regress out batch: https://github.com/satijalab/seurat/issues/3270 It is suggested to not regress out batch, and instead use a data integration method

# Check
pigs.split

## $`1.Saline`
## An object of class Seurat 
## 28275 features across 3983 samples within 2 assays 
## Active assay: SCT (13915 features, 3000 variable features)
##  1 other assay present: RNA
##  1 dimensional reduction calculated: pca
## 
## $`23.Ecoli`
## An object of class Seurat 
## 27486 features across 1668 samples within 2 assays 
## Active assay: SCT (13126 features, 3000 variable features)
##  1 other assay present: RNA
##  1 dimensional reduction calculated: pca

Integration

Run harmony

Condition-specific clustering of cells indicates that we need to integrate the cells across conditions to ensure that cells of the same cell type cluster together.

To integrate, use the shared highly variable genes from each condition identified using SCTransform. Then, integrate conditions to overlay cells that are similar or have a “common set of biological features” between groups.

Now, using our SCTransform object as input, let’s perform the integration across conditions.

First, we need to specify that we want to use all of the 3000 most variable genes identified by SCTransform for the integration. By default, this function selects the top 2000 genes.

# Choose the features to use when integrating multiple datasets. 
# will use nfeatures as 3000 as defined by running SCTransform above
var.features <- SelectIntegrationFeatures(object.list = pigs.split, 
                                          nfeatures = 3000)

# merge the pigs
pigs.sct.merged <- merge(x = pigs.split[[1]],
                         y = c(pigs.split[[2]]),
                         project = paste0("Ecoli Pigs ", tissue, " Single Cell"))

# define the variable features 
VariableFeatures(pigs.sct.merged) <- var.features

# run PCA on the merged object
pigs.sct.merged <- RunPCA(object = pigs.sct.merged, assay = "SCT")

# run harmony to harmonize over samples 
pigs.integrated <- RunHarmony(object = pigs.sct.merged, 
                              group.by.vars = "sample", 
                              assay.use = "SCT",
                              reduction = "pca", 
                              plot_convergence = TRUE)

Check output

# first put the pigs back in place
Idents(pigs.integrated) <- pigs.integrated$sample
pigs.integrated$treatment <- factor(pigs.integrated$treatment,
                                    levels = c("Saline","Ecoli"))
pigs.integrated$sample <- factor(pigs.integrated$sample,
                                 levels = c("1.Saline","23.Ecoli"))
# check the embedding
harmony_embeddings <- Embeddings(pigs.integrated, 'harmony')
harmony_embeddings[1:5, 1:5]

##                              harmony_1 harmony_2 harmony_3  harmony_4
## 1.Saline_AAACCCAAGGGCAACT-1   3.528350  2.227911 -0.539476 -3.3992357
## 1.Saline_AAACCCACACGCTGTG-1   7.239787  2.590364 -2.274517 -0.7715327
## 1.Saline_AAACCCAGTCCTGAAT-1   5.811397 -4.977975  1.365228  2.5140437
## 1.Saline_AAACCCAGTCGTTGGC-1 -46.150818 12.528053 -4.853579  6.2060722
## 1.Saline_AAACCCAGTGATGGCA-1   5.090145  1.511635 -2.428197 -0.5712383
##                                harmony_5
## 1.Saline_AAACCCAAGGGCAACT-1   1.82639385
## 1.Saline_AAACCCACACGCTGTG-1   0.07902291
## 1.Saline_AAACCCAGTCCTGAAT-1  -9.99066603
## 1.Saline_AAACCCAGTCGTTGGC-1 -12.52232189
## 1.Saline_AAACCCAGTGATGGCA-1  -1.13015936

# check the PCA plot
p1 <- DimPlot(object = pigs.integrated, 
              reduction = "harmony",
              group.by = "treatment", 
              cols = treatment_colors) + NoLegend()
p2 <- VlnPlot(object = pigs.integrated, 
              features = "harmony_1", 
              group.by = "treatment", 
              pt.size = 0, 
              cols = treatment_colors) + NoLegend()
plot_grid(p1,p2)

Top variable features

Top 20 variable features

top20 <- pigs.integrated@assays$SCT@var.features[1:20]
top20

##  [1] "CCL5"               "CXCL9"              "GNLY"              
##  [4] "CXCL10"             "ENSSSCG00000026043" "CD163"             
##  [7] "CXCL11"             "HPS5"               "ACSL4"             
## [10] "GZMH"               "SDS"                "ENSSSCG00000016184"
## [13] "MRC1"               "STEAP4"             "CCL4"              
## [16] "SPP1"               "IL1R2"              "CCL8"              
## [19] "ENSSSCG00000006588" "S100A12"

PCA plot

After integration, to visualize the integrated data we can use dimensionality reduction techniques, such as PCA and Uniform Manifold Approximation and Projection (UMAP). While PCA will determine all PCs, we can only plot two at a time. In contrast, UMAP will take the information from any number of top PCs to arrange the cells in this multidimensional space. It will take those distances in multidimensional space, and try to plot them in two dimensions. In this way, the distances between cells represent similarity in expression.

To generate these visualizations with the harmony output, use reduction = “harmony”

# Plot PCA
pca1 <- DimPlot(pigs.integrated,
        reduction = "harmony",
        split.by = "treatment",
        group.by = "treatment",
        cols = treatment_colors)
pca1

pca2 <- DimPlot(pigs.integrated,
        reduction = "harmony",
        split.by = "sample",
        group.by = "sample",
        cols = sample_colors)
pca2

pca3 <- DimPlot(pigs.integrated,
        reduction = "harmony",
        group.by = "sample",
        shuffle = TRUE)
pca3

Find significant PCs

To overcome the extensive technical noise in the expression of any single gene for scRNA-seq data, Seurat assigns cells to clusters based on their PCA scores derived from the expression of the integrated most variable genes, with each PC essentially representing a “metagene” that combines information across a correlated gene set. Determining how many PCs to include in the clustering step is therefore important to ensure that we are capturing the majority of the variation, or cell types, present in our dataset.

# Printing out the most variable genes driving PCs
print(x = pigs.integrated[["pca"]], 
      dims = 1:10, 
      nfeatures = 10)

## PC_ 1 
## Positive:  CALCR, APOE, ENSSSCG00000031640, RNF128, C1QB, ENSSSCG00000012238, SMAP2, C1QA, P2RY12, ENSSSCG00000033041 
## Negative:  ENSSSCG00000037142, RAC2, ETS1, CD3E, ENSSSCG00000031838, ENSSSCG00000034448, CCL5, GZMH, CD2, TPT1 
## PC_ 2 
## Positive:  ETS1, GZMH, CALCR, CD3E, CCL5, RAC2, S100A6, CTSW, CST7, RNF128 
## Negative:  HLA-DRA, SLA-DQB1, ENSSSCG00000001455, ENSSSCG00000001456, ENSSSCG00000033262, CD163, MRC1, ENSSSCG00000013788, ENSSSCG00000006350, ENSSSCG00000013100 
## PC_ 3 
## Positive:  MT3, CLU, S100A1, ENSSSCG00000024009, GPM6A, CAMK2A, ENSSSCG00000032480, ENSSSCG00000014540, MAP1B, MAP1A 
## Negative:  TMSB10, C1QB, CALCR, CD3E, GZMH, CCL5, RPS8, CD2, ENSSSCG00000031838, RPL10 
## PC_ 4 
## Positive:  ENSSSCG00000026043, ENSSSCG00000016184, ACSL4, SDS, MXD1, ENSSSCG00000022925, ENSSSCG00000015664, CSF3R, SELL, STEAP4 
## Negative:  ENSSSCG00000033262, HLA-DRA, SLA-DQB1, ENSSSCG00000001455, ENSSSCG00000001456, CD163, MRC1, CCL14, ENSSSCG00000006350, ENSSSCG00000004970 
## PC_ 5 
## Positive:  S100A8, ENSSSCG00000006588, TPT1, ENSSSCG00000033310, ENSSSCG00000004970, ENSSSCG00000014540, RPS12, RPS23, RPS20, RPL23 
## Negative:  ENSSSCG00000024973, CCL5, HLA-DRA, ENSSSCG00000001455, SLA-DQB1, ENSSSCG00000033089, GHSR, GZMH, GNLY, CD2 
## PC_ 6 
## Positive:  HMGB2, STMN1, HMMR, ENSSSCG00000009327, H2AFZ, TOP2A, DUT, TPX2, SMC2, KIF11 
## Negative:  ENSSSCG00000024973, ENSSSCG00000033089, CXCL10, GHSR, HLA-DRA, SAMD9, ENSSSCG00000001455, CD2, CCL5, ENSSSCG00000001229 
## PC_ 7 
## Positive:  HMGB2, STMN1, H2AFZ, DUT, TPX2, HMMR, TUBA1B, TOP2A, ENSSSCG00000035544, KIF11 
## Negative:  MGP, FN1, SLC16A9, SPARC, ENSSSCG00000005638, ENSSSCG00000000261, JAM2, ENSSSCG00000003839, CALD1, COL3A1 
## PC_ 8 
## Positive:  CD163, ENSSSCG00000031640, SELENOP, MRC1, ENSSSCG00000013788, DAB2, ENSSSCG00000006350, CCL14, ENSSSCG00000013100, STAB1 
## Negative:  CXCL10, ENSSSCG00000024973, GHSR, ISG15, ENSSSCG00000033089, ENSSSCG00000009240, IFIT1, ISG12(A), IFIT2, MX1 
## PC_ 9 
## Positive:  GNLY, KLF2, CCL5, GZMA.1, FLNA, NKG7, ENSSSCG00000037364, ADGRG1, ENSSSCG00000023584, S1PR5 
## Negative:  IL7R, CD3E, SIT1, CXCR3, ENSSSCG00000016475, SPOCK2, CD4, KCNA3, CD40LG, ENSSSCG00000029160 
## PC_ 10 
## Positive:  SELENOP, ENSSSCG00000031640, IL7R, LTC4S, CCL14, RNF128, TREM2, TMSB10, B2M, GHSR 
## Negative:  ENSSSCG00000006588, S100A8, ATP1B1, SPP1, MS4A7, MBNL3, GNLY, CCL5, NRP1, GZMH

Quantitative approach to an elbow plot - The point where the principal components only contribute 5% of standard deviation and the principal components cumulatively contribute 90% of the standard deviation. - The point where the percent change in variation between the consecutive PCs is less than 0.1%.

First metric

# Determine percent of variation associated with each PC
stdv <- pigs.integrated[["pca"]]@stdev
sum.stdv <- sum(pigs.integrated[["pca"]]@stdev)
percent.stdv <- (stdv / sum.stdv) * 100

# Calculate cumulative percents for each PC
cumulative <- cumsum(percent.stdv)

# Determine which PC exhibits cumulative percent greater than 90% and
# and % variation associated with the PC as less than 5
co1 <- which(cumulative > 90 & percent.stdv < 5)[1]
co1

## [1] 42

Second metric

# Determine the difference between variation of PC and subsequent PC
co2 <- sort(which(
  (percent.stdv[1:length(percent.stdv) - 1] - 
     percent.stdv[2:length(percent.stdv)]) > 0.1), 
  decreasing = T)[1] + 1

# last point where change of % of variation is more than 0.1%.
co2

## [1] 18

Choose the minimum of these two metrics as the PCs covering the majority of the variation in the data.

# Minimum of the two calculation
min.pc <- min(co1, co2)
min.pc

## [1] 18

Elbow plot

Use min.pc we just calculated to generate the clusters. We can plot the elbow plot again and overlay the information determined using our metrics:

# Create a dataframe with values
plot_df <- data.frame(pct = percent.stdv, 
           cumu = cumulative, 
           rank = 1:length(percent.stdv))

# Elbow plot to visualize 
  ggplot(plot_df, aes(cumulative, percent.stdv, label = rank, color = rank > min.pc)) + 
  geom_text() + 
  geom_vline(xintercept = 90, color = "grey") + 
  geom_hline(yintercept = min(percent.stdv[percent.stdv > 5]), color = "grey") +
  theme_bw()

UMAP

# Run UMAP
pigs.integrated <- RunUMAP(pigs.integrated,
                           dims = 1:min.pc,
                           reduction = "harmony",
                           n.components = 3) # set to 3 to use with VR

# plot UMAP and color based on treatment
DimPlot(pigs.integrated,
        group.by = "treatment",
        split.by = "treatment",
        shuffle = TRUE,
        cols = treatment_colors)

Cluster the cells

Seurat uses a graph-based clustering approach, which embeds cells in a graph structure, using a K-nearest neighbor (KNN) graph (by default), with edges drawn between cells with similar gene expression patterns. Then, it attempts to partition this graph into highly interconnected ‘quasi-cliques’ or ‘communities’ [Seurat - Guided Clustering Tutorial].

We will use the FindClusters() function to perform the graph-based clustering. The resolution is an important argument that sets the “granularity” of the downstream clustering and will need to be optimized for every individual experiment. For datasets of 3,000 - 5,000 cells, the resolution set between 0.4-1.4 generally yields good clustering. Increased resolution values lead to a greater number of clusters, which is often required for larger datasets.

The FindClusters() function allows us to enter a series of resolutions and will calculate the “granularity” of the clustering. This is very helpful for testing which resolution works for moving forward without having to run the function for each resolution.

# Determine the K-nearest neighbor graph
pigs.unannotated <- FindNeighbors(object = pigs.integrated, 
                                 assay = "SCT", # set as default after SCTransform
                                 reduction = "harmony",
                                 dims = 1:min.pc)

# Determine the clusters for various resolutions
pigs.unannotated <- FindClusters(object = pigs.unannotated,
                                 algorithm = 1, # 1= Louvain
                                 resolution = seq(0.1,0.8,by=0.1))

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5651
## Number of edges: 189697
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9359
## Number of communities: 6
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5651
## Number of edges: 189697
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9011
## Number of communities: 9
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5651
## Number of edges: 189697
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8805
## Number of communities: 10
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5651
## Number of edges: 189697
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8609
## Number of communities: 11
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5651
## Number of edges: 189697
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8444
## Number of communities: 13
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5651
## Number of edges: 189697
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8323
## Number of communities: 15
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5651
## Number of edges: 189697
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8215
## Number of communities: 16
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5651
## Number of edges: 189697
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8120
## Number of communities: 16
## Elapsed time: 0 seconds

colors <- c("dodgerblue","yellow","firebrick1","gold","gray40","lightgray",
            "cyan","chocolate4","pink","orange","darkgreen","purple","lightgray","tan")
DimPlot(pigs.unannotated,
        group.by = "SCT_snn_res.0.1",
        label = TRUE,
        cols = colors)

Explore resolutions

# 0.1
umap0.1 <- DimPlot(pigs.unannotated,
        group.by = "SCT_snn_res.0.1",
        label = TRUE)
umap0.1

# 0.2
umap0.2 <- DimPlot(pigs.unannotated,
        group.by = "SCT_snn_res.0.2",
        label = TRUE)
umap0.2

# 0.3
umap0.3 <- DimPlot(pigs.unannotated,
        group.by = "SCT_snn_res.0.3",
        label = TRUE)
umap0.3

# 0.4
umap0.4 <- DimPlot(pigs.unannotated,
        group.by = "SCT_snn_res.0.4",
        label = TRUE)
umap0.4

Clustering QC

Treatment, sample, phase

# treatment
u1 <- DimPlot(pigs.unannotated, 
        label = FALSE, 
        group.by = "SCT_snn_res.0.1",
        split.by = "treatment") +
  NoLegend()
u1

# sample
u2 <- DimPlot(pigs.unannotated, 
        label = FALSE,
        group.by = "SCT_snn_res.0.1",
        split.by = "sample") +
  NoLegend()
u2

# phase
u3 <- DimPlot(pigs.unannotated, 
        label = FALSE,
        group.by = "SCT_snn_res.0.1",
        split.by = "Phase") +
  NoLegend()
u3

Revisit QC metrics

# nCount
f1 <- FeaturePlot(pigs.unannotated, 
            features = "nCount_RNA",
            pt.size = 0.4, 
            order = TRUE,
            min.cutoff = 'q10',
            label = TRUE)
f1

# nFeature
f2 <- FeaturePlot(pigs.unannotated, 
            features = "nFeature_RNA",
            pt.size = 0.4, 
            order = TRUE,
            min.cutoff = 'q10',
            label = TRUE)
f2

# percent.mt
f3 <- FeaturePlot(pigs.unannotated, 
            features = "percent.mt",
            pt.size = 0.4, 
            order = TRUE,
            min.cutoff = 'q10',
            label = TRUE)
f3

# cell.complexity
f4 <- FeaturePlot(pigs.unannotated, 
            features = "cell.complexity",
            pt.size = 0.4, 
            order = TRUE,
            min.cutoff = 'q10',
            label = TRUE)
f4

f5 <- FeaturePlot(pigs.unannotated, 
            features = "S.Score",
            pt.size = 0.4, 
            order = TRUE,
            min.cutoff = 'q10',
            label = TRUE)
f5

f6 <- FeaturePlot(pigs.unannotated, 
            features = "G2M.Score",
            pt.size = 0.4, 
            order = TRUE,
            min.cutoff = 'q10',
            label = TRUE)
f6

Bulk DEGs

# cell.complexity
f5 <- FeaturePlot(pigs.unannotated, 
            features = c("ATF3","CCL2","CXCL10","ICAM1"),
            pt.size = 0.4, 
            order = TRUE,
            min.cutoff = 'q1',
            label = FALSE)
f5

#f6 <- FeaturePlot(pigs.unannotated, 
#            features = "KK-MALAT1",
#            pt.size = 0.4, 
#            order = TRUE,
#            min.cutoff = 'q1',
 #           label = TRUE)
#f6

Percent cells per cluster

pigs.unannotated@meta.data$seurat_clusters <- 
  pigs.unannotated@meta.data$SCT_snn_res.0.1

# treatment
b1 <- pigs.unannotated@meta.data %>%
  group_by(seurat_clusters, treatment) %>%
  dplyr::count() %>%
  group_by(seurat_clusters) %>%
  dplyr::mutate(percent = 100*n/sum(n)) %>%
  ungroup() %>%
  ggplot(aes(x=seurat_clusters,y=percent, fill=treatment)) +
  geom_col() +
  scale_fill_manual(values = treatment_colors) +
  ggtitle("Percentage of treatment per cluster")
b1

# sample
b2 <- pigs.unannotated@meta.data %>%
  group_by(seurat_clusters, sample) %>%
  dplyr::count() %>%
  group_by(seurat_clusters) %>%
  dplyr::mutate(percent = 100*n/sum(n)) %>%
  ungroup() %>%
  ggplot(aes(x=seurat_clusters,y=percent, fill=sample)) +
  geom_col() +
  #scale_fill_manual(values = sample_colors) +
  ggtitle("Percentage of sample per cluster")
b2

Number cells per cluster

treatment_ncells <- FetchData(pigs.unannotated, 
                     vars = c("ident", "treatment")) %>%
  dplyr::count(ident, treatment) %>%
  tidyr::spread(ident, n)
write.table(treatment_ncells, 
            paste0("../../results/nCells/",
                   treatment, "_",tolower(tissue),
                   "_cells_per_cluster_treatment.txt"),
            quote = FALSE, sep = "\t")

sample_ncells <- FetchData(pigs.unannotated, 
                     vars = c("ident", "sample")) %>%
  dplyr::count(ident,sample) %>%
  tidyr::spread(ident, n)
write.table(sample_ncells, 
            paste0("../../results/nCells/",
                   treatment, "_",tolower(tissue),
                   "_cells_per_cluster_sample.txt"),
            quote = FALSE, sep = "\t")