4. Single cell RNAseq

This document contain codes to process and analyse 10X Genomics gene expression data. Every project has different requirements, so tailor the pipeline according to the needs of the project. These are the key steps of the pipeline:

Align sequencing reads to the genome and quantify gene expression data
Import data into R and create Seurat object
Perform quality control
Dimensional reduction and cell clustering
Extended analysis
1. Differential expression analysis
2. Data integration and label transfer
3. Pseudotime analyses
4. Cross-data correlation
5. RNA velocity

Constructing custom CellRanger transcriptomes

Resource link:

Building a custom transcriptome

Some experiments contain transgenic markers that expresses mRNA at detectable amounts. This section describes the code to generate a custom CellRanger transcriptome that contain these transgenic mRNA(s).

You need:

FASTA file containing the sequence of the transgene. You may use sequence of the entire plasmid or just the sequence of the desired transgene (preferred). Make sure that the 3'UTR of the gene is included, as most of the reads will cluster within the 3'UTR. The FASTA headers ( > header) will serve as chromosome names, so try to keep it short and simple.
GTF transcriptome file. Use the same GTF that comes with CellRanger's pre-built reference

Create custom GTF file

We first need to create a custom GTF file that contain the coordinates of the transgene. The code below is done in R:

# load dependencies
library(tidyverse)
library(GenomicRanges)

# load CellRanger reference GTF
gtf <- import("/PATH/TO/genes.gtf")

## extract the first 3 lines of the GTF and modify the relebant metadata
### I left the values below as it is as an example
endo.gtf <- gtf %>% 
    as.data.frame() %>% 
    head(3) %>% 
    mutate(seqnames = "Ngn2_DSRed",          # set this to the same name as the FASTA sequence header
           start = 1,                        # if the FASTA file contain transgene sequence, this will be 1
           end = 3740,                       # coordinate of the end of transgene
           width = 3740,                     # length of transgene
           strand = "+",                     # keep to "+" strand
           source = "ENDO",                  # set to any character value
           gene_id = "ENDO.1",               # set to an easily recognizable character value
           gene_name = "Ngn2_DSRed",         # Name of gene. This will be displayed as feature after CellRanger run
           havana_gene = NA, 
           transcript_id = c(NA, "ENDO.1.1", "ENDO.1.1"), # Use same id as `gene_id`, with additional suffix
           transcript_name = c(NA, "Ngn2_DSRed.1","Ngn2_DSRed.1"), # Use same name as `gene_name`, with additional suffix
           havana_transcript = NA,
           exon_id = c(NA,NA,"Ngn2_DSred.1.1")) %>% # Use same name as `gene_name`, with additional suffix
    GenomicRanges::makeGRangesFromDataFrame(keep.extra.columns = T)
c(gtf, endo.gtf) %>% 
    rtracklayer::export("../outputs/GTF/mm10_Ngn2DSRed.gtf")

Create custom genome sequence file

Next, concatenate the genome fasta file with the transgene fasta file. This is done in Bash:

#
cat /PATH/TO/genome.fa /PATH/TO/transgene.fa > /PATH/TO/custom_genome.fa

Construct custom CellRanger ref

Lastly, construct the custom reference using cellranger mkref. Output the reference into the project directory.

cd /PATH/TO/PROJECT/outputs

# create folder to contain custom reference
mkdir -p cellrangerref && cd $_

cellranger mkref --genome=NEW_GENOME_NAME \
  --fasta=/path/to/custom_genome.fa \
  --genes=/path/to/custom_transcriptome.gtf \
  --nthreads=35

Aligning 10X FASTQ files using `cellranger count`

Resource link:

10X Single Cell CellRanger Software

cd /PATH/TO/PROJECT/DIRECTORY/outputs
mkdir -p cellrangerout && cd $_

cellranger count --id OUTNAME --sample SAMPLENAME \
    --fastqs /PATH/TO/FASTQ/DIRECTORY \
    --transcriptome /PATH/TO/TRANSCRIPTOME \
    --localcores=40 --localmem=120 --disable-ui

Aligning 10X 3' Multiplexing data using `cellranger multi`

Resource link:

10X Cell Multiplexing with cellranger multi

You need:

A CSV file containing the configuration of the Multiplex experiment. An example can be found at /media/cdn-bc/RAID/Projects/FH026_Youran_iNs/originals/Pilot_Multi_Config.csv

cd /PATH/TO/PROJECT/outputs

# run 
cellranger multi --id=Experiment_name \
        --csv=/path/to/csv \
        --localcores=40

Processing gene expression data using Seurat

Resource link:

Seurat processing workflow

These are the dependencies required for this part of the analysis

## load dependencies
library(tidyverse)
library(Seurat)
library(scDblFinder)
library(SingleCellExperiment)
library(SeuratWrappers)
library(SeuratDisk)
library(patchwork)
library(clustree)

Import count matrices

For every experiment, we will prepare a samples metadata file (CSV) that should contain at least the following columns:

Path : Path to the folder containing the gene count matrix
ID : Experiment name or ID that will be passed to orig.ident

The metadata file can contain any other information pertaining to the experiment which will be passed into the metadata dataframe of the Seurat object.

We start by importing the metadata file and create a list of Seurat objects for each experiment:

## Import samples meta file
samples.meta <- read_csv("/PATH/TO/SAMPLES.csv")

## check if all paths are present
all(file.exists(samples.meta$Path))

samples.seuratlist <- apply(samples.meta, 1, function(line) {
  
  in.file <- str_subset(list.files(line[["Path"]]), "matrix.h5$|.mtx.gz|matrix.mtx")
  
  
  # Determine whether to use Read10X_h5 or Read10X
  if (any(grepl("\\.h5$", in.file))) {
    # load H5
    sample.data <- Read10X_h5(file.path(line[["Path"]], "filtered_feature_bc_matrix.h5"))
  } else if (any(grepl("\\.mtx.gz$", in.file))) {
    # Load directory containing the matrix.mtx, features.tsv, and barcodes.tsv files provided by 10X.
    sample.data <- Read10X(line[["Path"]])
  } else if (any(grepl("\\matrix.mtx$", in.file))) {
    # Load directory containing the matrix.mtx, features.tsv, and barcodes.tsv files provided by 10X.
    sample.data <- ReadParseBio(line[["Path"]])
  }
  
  sample.data
  
  # Create Seurat object
  sample.seurat <- CreateSeuratObject(counts = sample.data,
                                      project = line[["ID"]])
  sample.seurat$experiment <- "RNA"

  # Adding additional metadata
  sample.seurat@meta.data <- mutate(sample.seurat@meta.data, !!!line)

  sample.seurat
})

names(samples.seuratlist) <- samples.meta$ID

Quantify QC metrics

Next, add more information on the quality of each cell by quantifying the following metrics:

Doublet annotation
Percent mitochondria
Percent ribosomal
Cell cycle phase

These values will be stored in the meta.data slot of each Seurat object

## annotate doublets using scDblFinder
samples.seuratlist <- lapply(samples.seuratlist, function(sample){
    sce <- scDblFinder(GetAssayData(sample, assay = "RNA", layer = "counts"))
    sample$scDblFinder.class <- sce$scDblFinder.class
    sample
})

## calculate percent mitochondria
samples.seuratlist <- lapply(samples.seuratlist, function(sample){
    sample$percent_mito <- PercentageFeatureSet(sample, pattern = "^MT-",assay = "RNA") #change regex accordingly
    Idents(sample) <- "experiment"
    sample
})

## plot mitochondria % distribution per experiment
do.call(wrap_plots,lapply(samples.seuratlist, function(sample){
    VlnPlot(sample, "percent_mito", group.by = "experiment", layer = "counts") +
        ggtitle(sample$orig.ident[1]) +
        scale_y_continuous(limits = c(0,65))
}))


## calculate percent ribosomal
samples.seuratlist <- lapply(samples.seuratlist, function(sample){
    sample$percent_ribo <- PercentageFeatureSet(sample, pattern = "^RP[SL]",assay = "RNA") 
    Idents(sample) <- "experiment"
    sample
})


## plot ribosome % distribution per experiment
do.call(wrap_plots,lapply(samples.seuratlist, function(sample){
    VlnPlot(sample, "percent_ribo", group.by = "experiment", layer = "counts") +
        ggtitle(sample$orig.ident[1]) +
        scale_y_continuous(limits = c(0,65))
}))


## quantify cell cycle phase
### run cell cycle scoring
cc.genes <- cc.genes.updated.2019

##################
# To convert cc.genes to other organism of choice, run the following:
cc.genes <- lapply(cc.genes, function(genes) {
   gprofiler2::gorth(genes, 
        target_organism = "mmusculus"   #change target_organism accordingly
   )$ortholog_name
})
##################


samples.seuratlist <- lapply(samples.seuratlist, function(sample){
	sample <- NormalizeData(sample)
    CellCycleScoring(sample, 
                     s.features = cc.genes$s.genes, 
                     g2m.features = cc.genes$g2m.genes, 
                     set.ident = FALSE)
})

do.call(wrap_plots,lapply(samples.seuratlist, function(sample){
    sample@meta.data %>% 
        ggplot(aes(x=Phase, fill = Phase)) + 
        geom_bar() +
        theme_minimal() +
        scale_fill_brewer(palette = "Set2") +
        ggtitle(sample$orig.ident[1]) 
}))

Normalise, scale and regress data

## Run SCTransform and regress out CC-score
samples.seuratlist <- lapply(samples.seuratlist, function(sample){
    SCTransform(sample, method = "glmGamPoi", 
    		vars.to.regress = c("S.Score","G2M.Score"), verbose = FALSE)
})

Filtering out poor quality cells

We performed a strigent 2-pass quality control:

1st-pass:
- Remove cells with >3 MAD lower number of genes than median
- Remove cells with >3 MAD higher percent mitochondria than median
2nd-pass:
- Remove clusters with concomitantly low features (<1000) and low mito % (<2.5)
- Remove clusters with concomitantly low features (<2000) and low rna count (<4000)

### check number of cells before qc
ncell_before.v <- lapply(samples.seuratlist, ncol)

## 1st qc pass using median absolute deviation from the median
qc.df <- do.call(rbind,lapply(samples.seuratlist, function(sample){
    feature.cutoff <- median(sample@meta.data$nFeature_SCT) - 3*mad(sample@meta.data$nFeature_SCT)
    feature.cutoff <- ifelse(feature.cutoff < 100, 100, feature.cutoff)
    mito.cutoff <- median(sample@meta.data$percent_mito) + 3*mad(sample@meta.data$percent_mito)
    mito.cutoff <- ifelse(mito.cutoff > 10, 10, mito.cutoff)
    return(data.frame(feature.cutoff = feature.cutoff,
                      mito.cutoff = mito.cutoff))
}))

samples.seuratlist.filt <- lapply(samples.seuratlist, function(sample) {
    sample.id <- as.character(sample$orig.ident[1])
    subset(sample, 
           subset = nFeature_RNA > qc.df[sample.id,]$feature.cutoff & percent_mito < qc.df[sample.id,]$mito.cutoff)
})


## check number of cells
ncell_after.v <- lapply(samples.seuratlist.filt, ncol)
unlist(ncell_before.v) - unlist(ncell_after.v)

# check metric distribution
do.call(wrap_plots,lapply(samples.seuratlist.filt, function(sample){
    VlnPlot(sample, "percent_mito", group.by = "experiment") +
        ggtitle(sample$orig.ident[1]) 
}))

do.call(wrap_plots,lapply(samples.seuratlist.filt, function(sample){
    VlnPlot(sample, "nFeature_RNA", group.by = "experiment") +
        ggtitle(sample$orig.ident[1]) 
}))

# overwrite original object
samples.seuratlist <- samples.seuratlist.filt
rm(samples.seuratlist.filt) 

## briefly normalize, scale, get variable features, and pca
samples.seuratlist <- lapply(samples.seuratlist, function(sample){
    SCTransform(sample, method = "glmGamPoi", verbose = FALSE)

})

## 2nd qc pass: remove low quality clusters
### Briefly cluster cells using higher resolution
samples.seuratlist <- lapply(samples.seuratlist, function(sample){
	sample <- RunPCA(sample)
    sample <- FindNeighbors(sample, dims = 1:50)
    sample <- FindClusters(sample, resolution = 1.8)
    sample
})

## check distribution of mito% and nFeature
do.call(wrap_plots,lapply(samples.seuratlist, function(sample){
    VlnPlot(sample, "percent_mito") +
        ggtitle(sample$orig.ident[1]) 
})) + plot_layout(ncol=1)

do.call(wrap_plots,lapply(samples.seuratlist, function(sample){
    VlnPlot(sample, "nFeature_RNA") +
        ggtitle(sample$orig.ident[1]) 
})) + plot_layout(ncol=1)


## first round of removing cluster
### done by removing clusters with concomitantly low features (<1000)
### and high mito % (>2.5)
samples.seuratlist <- lapply(samples.seuratlist, function(sample){
    clusters <- sample@meta.data %>% 
        group_by(seurat_clusters) %>% 
        summarise(maxfeat = max(nFeature_RNA), 
                  maxmito = max(percent_mito)) %>% 
        filter(maxfeat >1000 | maxmito < 2.5) %>% 
        pull(seurat_clusters)
    subset(sample, idents = clusters)
})

## check distribution of mito% and nFeature
do.call(wrap_plots,lapply(samples.seuratlist, function(sample){
    VlnPlot(sample, "percent_mito") +
        ggtitle(sample$orig.ident[1]) 
})) + plot_layout(ncol=1)

do.call(wrap_plots,lapply(samples.seuratlist, function(sample){
    VlnPlot(sample, "nFeature_RNA") +
        ggtitle(sample$orig.ident[1]) 
})) + plot_layout(ncol=1)


## recluster and reperform qc
samples.seuratlist <- lapply(samples.seuratlist, function(sample){
    sample <- FindNeighbors(sample, dims = 1:50)
    sample <- FindClusters(sample, resolution = 1.8)
    sample
})

## second round of removing cluster
### done by removing clusters with concomitantly low features (<2000)
### and low rna count (<4000)
samples.seuratlist <- lapply(samples.seuratlist, function(sample){
    clusters <- sample@meta.data %>% 
        group_by(seurat_clusters) %>% 
        summarise(maxfeat = max(nFeature_RNA), 
                  maxcount = max(nCount_RNA)) %>% 
        filter(maxfeat >2000 | maxcount > 4000) %>% 
        pull(seurat_clusters)
    subset(sample, idents = clusters)
})

## check distribution of nCount and nFeature
do.call(wrap_plots,lapply(samples.seuratlist, function(sample){
    VlnPlot(sample, "nCount_RNA") +
        ggtitle(sample$orig.ident[1]) 
})) + plot_layout(ncol=1)

do.call(wrap_plots,lapply(samples.seuratlist, function(sample){
    VlnPlot(sample, "nFeature_RNA") +
        ggtitle(sample$orig.ident[1]) 
})) + plot_layout(ncol=1)

## Final SCTransform
samples.seuratlist <- lapply(samples.seuratlist, function(sample){
    SCTransform(sample, method = "glmGamPoi", 
    		vars.to.regress = c("S.Score","G2M.Score"), verbose = FALSE)
})

Dimensional reduction

The following code reduces the complexity of the data and projects the cell onto a 2-dimensional cartesian coordinate.

## Perform UMAP reduction
samples.seuratlist <- lapply(samples.seuratlist, function(sample){
	sample <- RunPCA(sample)
	sample <- RunTSNE(sample, dims = 1:50)
    sample <- RunUMAP(sample, dims = 1:50)
    sample 
})

Cell clustering

Now that we have high quality cells, we can cluster cells with similar transcriptome. The number of obtained clusters can be controlled using the resolution parameter, and the optimal number of clusters vary between experiments.

The code below performs cell clustering over a range of resolutions. The optimal resolution can be determined visually using the clustree package or by calculating the weighted sum-of-squares at each resolution.

## Determine the best clustering resolution
### Set range of resolutions
res.range <- seq(0.4, 3, 0.2)

samples.seuratlist <- lapply(samples.seuratlist, function(obj){
    DefaultAssay(obj) <- "SCT"
    obj <- FindNeighbors(obj, dims = 1:50)
    for(res in res.range){
        obj <- FindClusters(obj, resolution = res)
    }
    obj
})

# plot clustree for each sample
lapply(samples.seuratlist, function(dat){
    clustree::clustree(dat, prefix = "SCT_snn_res.")   # change the prefix accordingly
})  

### Calculate weighted sum of squares for each experiment
wss.all <- do.call(bind_rows,lapply(samples.seuratlist, function(sample){
  wss.out <- do.call(bind_rows, lapply(res.range, function(k){

    # get centroid per cluster
    centroids <- aggregate(Embeddings(sample, reduction = "pca"),
          sample[[paste0("SCT_snn_res.", k)]], 
          mean) %>% 
    column_to_rownames(paste0("SCT_snn_res.", k))
    diffsq <- (Embeddings(sample, reduction = "pca") - 
        centroids[sample@meta.data[[paste0("SCT_snn_res.", k)]],])^2
    sdsq <- colSds(Embeddings(sample, reduction = "pca"))^2
    wdiff <- diffsq/sdsq
    
    wss <- sum(sqrt(rowSums(wdiff)))
    data.frame(clust=k, wss = wss)
 }))
  
  wss.out %>% 
    mutate(id = sample$orig.ident[1])
}))

### Plot out WSS vs resolution
wss.all %>% 
  group_by(id) %>% 
  mutate(scaled_wss = scale(wss)) %>% 
  ggplot(aes(x=clust,y=scaled_wss,group=id, colour=id)) +
  geom_line()

### Cluster for real
res <- 1 #change accordingly
samples.seuratlist <- lapply(samples.seuratlist, function(sample){
	sample <- FindNeighbors(sample, verbose = F, dims = 1:50)
    sample <- FindClusters(sample, resolution = res, verbose = F)
})

lapply(samples.seuratlist, DimPlot)

Navigation Menu

4. Single cell RNAseq

Constructing custom CellRanger transcriptomes

Create custom GTF file

Create custom genome sequence file

Construct custom CellRanger ref

Aligning 10X FASTQ files using `cellranger count`

Aligning 10X 3' Multiplexing data using `cellranger multi`

Processing gene expression data using Seurat

Import count matrices

Quantify QC metrics

Normalise, scale and regress data

Filtering out poor quality cells

Dimensional reduction

Cell clustering

Differential gene expression analysis

Data integration and label transfer

CCA [Seurat]

Scanorama

Harmony

Trajectory analyses

Monocle3

Slingshot

Quantifying Gene Specficity Index (SI)

RNA velocity analyses

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Clone this wiki locally

4. Single cell RNAseq

Constructing custom CellRanger transcriptomes

Create custom GTF file

Create custom genome sequence file

Construct custom CellRanger ref

Aligning 10X FASTQ files using cellranger count

Aligning 10X 3' Multiplexing data using cellranger multi

Processing gene expression data using Seurat

Import count matrices

Quantify QC metrics

Normalise, scale and regress data

Filtering out poor quality cells

Dimensional reduction

Cell clustering

Differential gene expression analysis

Data integration and label transfer

CCA [Seurat]

Scanorama

Harmony

Trajectory analyses

Monocle3

Slingshot

Quantifying Gene Specficity Index (SI)

RNA velocity analyses

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Clone this wiki locally

Aligning 10X FASTQ files using `cellranger count`

Aligning 10X 3' Multiplexing data using `cellranger multi`