Skip to content

Latest commit

 

History

History
564 lines (419 loc) · 21.6 KB

README.md

File metadata and controls

564 lines (419 loc) · 21.6 KB
Raw

R-CMD-check


limiric_logo

Contents

Description | Installation | Basic usage | Extended usage | Output explained | Input file format | Downstream

Description


limiric is a quality control package for identifying damaged cells in scRNA-seq data. Many scRNA-seq quality control tools are designed to identify droplets that contain more than one cell (doublets) or no cells (empty droplets), but few are specialised in identifing those that contain damaged cells. Damaged cell detection is more often achieved by setting thresholds for metrics such as average mitochondrial gene expression and UMI counts per barcode. These thresholds, even when dynamically calculated, vary across samples, cell types, tissues, treatment conditions, and species, resulting in many true cells being excluded from downstream analysis. But searching for true damaged droplets in a sample-specific manner is tedious and not always intuitive, leading to user-defined filtering that lacks reproducibility.


limiric automates the sample-specific detection of damaged cells in one fast acting and highly reproducible function. It operates on the basic principle that damaged and healthy cells can be differentiated by the complexity of their mitochondrial and ribosomal gene expression, and that this complexity can be greatly simplified by clustering in lower dimensional space. This is achieved by applying the Louvain algorithm to a shared nearest neighbor graph constructed for each droplet using the expression of the Ensembl-defined subset of mitochondrial and ribosomal genes.


In addition to predicting damaged cells, limiric can perform other pre-processing tasks including removing red blood cells, correcting for ambient RNA (SoupX), and isolating immune cells. There is also the option to combine results with DropletQC, a community-accepted package that identifies empty and damaged droplets by nuclear fraction scores. Combining the diagnostic power of these nonredundant quality control metrics increases confidence in true damaged cell detection. This being said, limiric performs effectively without DropletQC, so this feature- while advised- is not essential.


The package was developed around the community standard Seurat suite can be incorporated seemlessly into a user's pre-existing Seurat workflow. However, the main output of limiric is a standard, two-column csv that can be used in any platform alongside any exisiting or future scRNA-seq analysis tool.



Installation

Prerequisites

limiric makes use of the following packages


Of these, the following need to be made available in your R environment before limiric can be installed. If not already, you can install them together with the following code

packages <- c("cowplot", "devtools", "dplyr", "ggplot2", "Matrix", "png", "Seurat", "SoupX")

for (pkg in packages) {
    if (!requireNamespace(pkg, quietly = TRUE)) {
      install.packages(pkg)
    }
  }

Be sure to load these into your environment before continuing with limiric installation

packages <- c("cowplot", "devtools", "dplyr", "ggplot2", "Matrix", "png", "Seurat", "SoupX")

for (pkg in packages) {
  if (!require(pkg, character.only = TRUE)) {
    library(pkg)
  }
}

Please check the associated documentation if problems occur in the installation of any of the prerequisite packages



Package installation

After all the prerequisites are installed, you can install the latest development version of limiric from CRAN using

install.packages("limiric")

Alternatively, the package can be installed from GitHub using the devtools package:

devtools::install_github("AlicenJoyHenning/limiric", build_vignettes = TRUE)

Verify installation

To ensure limiric has correctly installed, run the following to see if you can view the package vignette and the function help page.

library(limiric)
help(package = "limiric")
?limiric()

Additionally, you can perform a test run using a small example dataset stored in the testrun directory of this repository, also found here. If you need assistance setting this up, please see Test run for more details.




Quickstart

Basic Usage

  1. Detect damaged cells in a sample using the filtered alignment files (barcodes.tsv.gz, features.tsv.gz, counts.mtx.gz). All you need to input is your sample name, the directory where your alignment files are stored, and a directory where the limiric output can be created.
SRR1234567 <- limiric(
    project_name  = "SRR1234567",
    filtered_path = "/home/user/alignment/SRR1234567/filtered/",
    output_path   = "/home/user/alignment/limiric/"
)

NB Please keep the following in mind:

  • Your files must be zipped, see Input file format for more
  • The above usage assumes the sample is of human origin (See parameter: organism)
  • When compiling the data, all genes will be retained (See parameter: min_cells)
  • Red blood cells will automatically be removed from the sample before damaged cell detection (See parameter: filter_rbc)
  • Your output Seurat object will be filtered, containing only undamaged cells (See parameter: filter_output)


  1. Alternatively, you can use a Seurat object as input
SRR1234567 <- limiric(
    project_name  = "SRR1234567",
    seurat_input  = seurat_object,
    output_path   = "/home/user/alignment/limiric/"
)


  1. If you have more than one sample, you may want to define a sample list instead of running each individually. You can do this through the sample_list parameter.
sample_list <- list(

    list(project_name  = "SRR1234567",
         filtered_path = "/home/user/alignment/SRR1234567/filtered/",
         output_path   = "/home/user/alignment/limiric/"),
    
    list(project_name  = "SRR1234568",
         filtered_path = "/home/user/alignment/SRR1234568/filtered/",
         output_path   = "/home/user/alignment/limiric/"),
    
    list(project_name  = "SRR1234569",
         filtered_path = "/home/user/alignment/SRR1234569/filtered/",
         output_path   = "/home/user/alignment/limiric/")
)

GSE1234567 <- limiric(sample_list = sample_list)

NB Please note that you must define your sample list with each parameter specified, i.e. limiric won't know what to do if your list looks like this :

sample_list <- list(

   list("SRR1234567",
        "/home/user/alignment/SRR1234567/filtered/",
        "/home/user/alignment/limiric/"),
   list("SRR1234568",
        "/home/user/alignment/SRR1234568/filtered/",
        "/home/user/alignment/limiric/"),
   list("SRR1234569",
        "/home/user/alignment/SRR1234569/filtered/",
        "/home/user/alignment/limiric/")
 )

Test run

Download the matrix.mtx, barcodes.tsv. and features.tsv files from this location (files) and store them in a local directory. For demonstration purposes, we will assume this directory is /home/user/scRNA-seq/testrun/. From here, you can run the basic limiric function in your R environment where the terminal output should look as indicated below. For further verification, ensure the output Testrun_CellQC.png is identicial to that shown here.

testrun <- limiric(
    project_name  = "Testrun",
    filter_rbc    = FALSE,
    filtered_path = "/home/user/scRNA-seq/testrun/",
    output_path   = "/home/user/scRNA-seq/testrun/"
)

# Beginning limiric analysis for TestRun ...
# ✔ Seurat object created
# ✔ limiric damaged cell predictions
# ✔ limiric analysis complete.

Note: This shouldn't take more than 10 seconds to run but this may vary depending on your machine.



Output explained

Before interpretting limiric's outputs, it may be helpful to understand a bit of what's going on in the background.

As you know, the magic of scRNA-seq revolves around the identification of distinct cell populations. In majority of scRNA-seq analysis workflows, distinct populations are identified according to the clusters that form when dimensionality reduction is performed on the most variable genes in a sample. For instance, the typical Seurat workflow uses the top two or three thousand variable genes.

In our case, instead of trying to isolate distinct cell-type associated populations, we are only interested in two broad populations: damaged and healthy cells. In general, most of what we know about how these two populations differ can be summarized by looking at the expression of mitochondrial and ribosomal genes. Thus, if we perform dimensionality reduction and clustering using only these genes, we expect to see a division of cells into two clusters associated with the cell's status as either damaged or healthy.

This is the basis of the low dimension mitochondrial & ribosomal clustering, aka limiric, algorithm.

Output directory

The limiric output will be created in the output_pathdirectory with the following structure

output_path/
|
├── RBCQC
|
├── CellQC
|
└── Filtered

RBCQC

Before damaged cells can be identified, limiric removes red blood cells, or cells that are highly contaminated with haemoglobin, from your data. This is done under the assumption that red blood cells will not be informative to your study. If this assumption should not be true, you can avoid this filtering using filter_rbc = FALSE.

NB Given many scRNA-seq protocols, such as the 10X Genomics protocol, advise for globin treatment, we hope for the contamination percentage to be as low as possible.


The output here comes in the form of a scatter plot where the red blood cells are coloured in blue. The percentage of the total cells that were removed is also given in the plot. Scatter plot


CellQC

This directory contains core diagnostic plots for the limiric damaged cell detection algorithm. As shown below, you will see four tSNE plots. In this example dataset, you can see the cells divide into two distinct clusters; but which contain damaged cells and which contain healthy cells?


From literature we know that damaged cells have high mitochondrial expression and low ribosomal expression. This is largely because damaged cells, such as those undergoing apoptosis, are characterised by compromised cell membranes. In this case, free cytoplasmic RNA- like ribosomal RNA- is more likely to have escaped before being captured and sequenced, meaning its expression will be low. While the RNA enclosed within the mitochondria, which has its own membrane, is retained, meaning its expression will be high. But many other factors likely contribute to this phenomenon, including impaired ribosomal biogenesis in damaged cells.


Now answering the question of which cluster is which is easy with the limiric tSNEs where the mitochondrial and ribosomal gene expression of each cell is made visible.



  • A Mitochondrial gene expression tSNE shows that cells in the smaller cluster express mitochondrial genes at a very high level, suggesting they are likely damaged.

  • The same conclusion can be reached looking at the Ribosomal gene expression tSNE where cells in the smaller cluster express ribosomal genes at a very low level.

  • The Complexity score measures the total number of mitochondrial and ribosomal genes expressed in a cell. Cells that express a large number of these genes are more likely to be metabolically active, undamaged cells while those that only express a few are likely not. In the Complexity score tSNE, the smaller cluster contains cells with very low complexities. This verifies by another metric that these cells are likely damaged.



Together, these metrics are used by the limiric algorithm to annotate the damaged cells, as shown in the fourth and final tSNE. In summary, by looking at the limiric tSNE plots, we can immediately tell that we have a small cluster of cells that are likley damaged and should be removed from the sample. This makes sense as, unless something went drastically wrong during sample preparations, you expect there to be far more healthy than damaged cells in your scRNA-seq data.



Filtered


The main output of the limiric function comes in the form of a barcodes.csv containing annotations for each cell barcode of the input data. This can easily be incorporated into an existing scRNA-seq analysis workflow for filtering (see example).

project_name_barcodes.csv

Barcode Limiric
AAACCTGAGATAGGAG cell
AAACCTGAGCTATGCT cell
AAACCTGAGCTGTTCA damaged
AAACCTGCACATTAGC damaged
... ...


Additionally, limiric will output a Seurat object with ready-filtered barcodes for seamless integration at the start of a Seurat pre-processing workflow.



Extended usage

Slightly-less Basic Usage

1. Perform ambient RNA correction

Standard good practice highly advises to correct for ambient RNA in your scRNAseq data. If you haven't already performed some kind of correction, the SoupX parameter allows you to do so. This will occur before anything else in the limiric workflow.

SRR1234567 <- limiric(
    project_name  = "SRR1234567",
    filtered_path = "/home/user/alignment/SRR1234567/filtered/",
    soupx        = TRUE,
    raw_path      = "/home/user/alignment/SRR1234567/raw/",
    output_path   = "/home/user/alignment/limiric/"
)

2. Isolate immune cells

If you have a sample where only the immune cells are of interest, include the following isolate_cd45 parameter. This will isolate the immune cells present in the sample, then identify damaged cells.

SRR1234567 <- limiric(
    project_name  = "SRR1234567",
    filtered_path = "/home/user/alignment/SRR1234567/filtered/",
    isolate_cd45  = TRUE,
    output_path   = "/home/user/alignment/limiric/"
)

NB This will change your output directory structure by adding a new IMCQC layer

output_path/
├── CellQC
|
├── IMCQC
|
├── RBCQC
|
└── Filtered
This output, like RBCQC, contains a scatter plot showing the removed cell in blue and the retained cells in grey. Scatter plot


3. Combine limiric annotations with DropletQC

Detect damaged cells and compare results with those from DropletQC.

SRR1234567 <- limiric(
    project_name  = "SRR1234567",
    filtered_path = "/home/user/alignment/SRR1234567/filtered/",
    droplet_qc    = TRUE,
    velocyto_path = "/home/user/alignment/velocyto/",
    output_path   = "/home/user/alignment/limiric/"
)

NB This will change your output directory structure by adding a new DropletQC layer

output_path/
├── CellQC
|
├── DropletQC
|
├── RBCQC
|
└── Filtered

This will output a scatter plot and tSNE showing the cells annotated as damaged by both limiric and DropletQC. Scatter plot

4. Combine previous condition

Perform ambient RNA correction with SoupX, filter red blood cells, isolate immune cells, detect damaged cells, and compare against DropletQC.

SRR1234567 <- limiric(
    project_name  = "SRR1234567",
    filtered_path = "/home/user/alignment/SRR1234567/filtered/",
    soupx         = TRUE,
    raw_path      = "/home/user/alignment/SRR1234567/raw/",
    droplet_qc    = TRUE,
    isolate_cd45  = TRUE,
    velocyto_path = "/home/user/alignment/velocyto/",
    output_path   = "/home/user/alignment/limiric/"
)

NB This will mean your output directory looks like this

output_path/

├── CellQC
|
├── DropletQC
|
├── IMCQC
| 
├── RBCQC
| 
└── Filtered


5. Process multiple samples with the same conditions as in Example 4
sample_list <- list(
    list(project_name  = "SRR1234567",
         filtered_path = "/home/user/alignment/SRR1234567/filtered/",
         soupx         = TRUE,
         raw_path      = "/home/user/alignment/SRR1234567/raw/",
         droplet_qc    = TRUE,
         isolate_cd45  = TRUE,
         velocyto_path = "/home/user/alignment/velocyto/",
         output_path   = "/home/user/alignment/limiric/"),
    
    list(project_name  = "SRR1234568",
         filtered_path = "/home/user/alignment/SRR1234568/filtered/",
         soupx         = TRUE,
         raw_path      = "/home/user/alignment/SRR1234568/raw/",
         droplet_qc    = TRUE,
         isolate_cd45  = TRUE,
         velocyto_path = "/home/user/alignment/velocyto/",
         output_path   = "/home/user/alignment/limiric/"),
    
    list(project_name  = "SRR1234569",
         filtered_path = "/home/user/alignment/SRR1234569/filtered/",
         soupx         = TRUE,
         raw_path      = "/home/user/alignment/SRR1234569/raw/",
         droplet_qc    = TRUE,
         isolate_cd45  = TRUE,
         velocyto_path = "/home/user/alignment/velocyto/",
         output_path   = "/home/user/alignment/limiric/")
)

GSE1234567 <- limiric(sample_list = sample_list)


Input file format

If your alignment output files are not zipped (end with a .gz extension), you will need to find a way to do this before using limiric. If you have a Linux or mac machine you can open the terminal in the directory where your files are stored and use gzip


cd path/to/directory
gzip *

This assumes that your files are the only items inside the directory. As with the standard output of many alignment algorithms such as STARsolo and CellRanger, each sample should be stored in its own directory with the following file naming convention :

path/
|
├── matrix.mtx
|
├── barcodes.tsv
|
└── features.tsv

If you have a Windows machine, the simplest solution is to install Windows Subsystem for Linux 🐧 which creates a new terminal environment for you with Linux capabilites. From there, you can do the same as above. Note that WSL uses a different path structure to access file systems where Windows directories are mounted under /mnt/

cd /mnt/c/Users/path/to/directory
gzip *


Downstream

Adding barcode annotations to pre-existing Seurat object

# Add output to pre-existing Seurat object
limiric_annotations <- read.csv2("path/to/project_name_barcodes.csv",
                      sep = ",",
                      col.names = c("barcodes", "limiric"))

seurat@meta.data$limiric <- limiric_annotations$limiric[match(rownames(seurat@meta.data), limiric_annotations$barcodes)]

# Using pre-existing dimensionality reductions, visualise the limiric annotations
limiric_visual <- DimPlot(seurat, group.by = limiric) 

# Filter cells marked as damaged by limiric
seurat_filtered <- subset(seurat, limiric == "cell")
seurat_filtered$limiric <- NULL # once used, remove the column