My Profile Photo

alevin-fry-tutorials


Tutorials for using the alevin-fry single-cell RNA-seq pipeline


Spatial single-cell quantification with alevin-fry

In this tutorial we will look at how to process a spatial gene expression (Visium) experiment using an alevin-fry based pipeline. Note : This tutorial is meant to mimic the original tutorial for spatial expression analysis with alevin written by Avi Srivastava. Thus, most of the descriptive text and commands are taken directly from that tutorial. However, here we will be analyzing the data using the alevin-fry pipeline instead of alevin.

af_spatial.utf8

Introduction

Slide-seq V1-V2 (Rodriques, Samuel G., et al., Stickels, Robert R., et al., 2020) and 10x Genomics’ Visium Gene Expression has enabled the transcriptome-wide measurements of the molecular signals in a tissue with spatial localization at single-cell level. From the perspective of quantification, these spatial single-cell technologies require the generation of a gene-count matrix for each spot on the tissue, which are identified by spatially marked cellular barcodes. Based on the size of a spot & other technological limitations, each spot potentially contains multiple cells and, unlike single-cell sequencing, a single barcode can represent multiple cells. In this tutorial we generate the spatially-resolved gene-count matrix for each spot using alevin-fry and we visualize it using Seurat.

In this tutorial, we describe how to preprocess this data using alevin-fry, and how to load the data into R using Seurat and perform some basic analyses. Alevin-fry makes it easy to get from the raw FASTQ files to an accurate count matrix fast. For example, processing the data for this tutorial took less than 15 minutes (using 16 threads, and excluding index construction). The sections below contain code chunks executed in different languages, indicated with different background colors. R code is represented by grey boxes and shell commands by blue boxes.

Index the reference sequences

In order to quantify the abundances with alevin-fry, we need to generate the index of the reference sequences. We already saw how to do this for the reference transcriptome in earlier tutorials, and here we will assume that you’ve built the splici transcriptome reference for mouse. Specifically, you can obtain the ref-gex-mm10-2020-A reference from the 10x website, and use the make_splici_txome function from usefulaf utilities to generate the transcriptome and 2 and 3 column transcript-to-gene mappings. Then, index the resulting transcriptome_splici.fa with salmon to obtain the index we’ll use for mapping. In subsequent steps, we’ll assume the index resides in refdata-gex-mm10-2020-A/sal_index, but just change the index path below if you have it somewhere different.

Download the raw spatial sequencing data

We download the FASTQ files for the Mouse Brain Serial Section 1 Visium experiment. The raw reads files are available on the website with the hyperlinked FASTQs. However, since there is a registration requirements from the source, we are not providing the direct download link. Once downloaded the reads would look like the following:

tar -xvf V1_Mouse_Brain_Sagittal_Anterior_fastqs.tar

We also need the spatial imaging data which contains the spatial coordinates for each cellular barcode on the tissue slide (similar data can be found in the Slides-seq V1-V2). The data can be found on the same web link under the name Spatial imaging data. Once downloaded the content of the file looks as follows:

tar -xvzf V1_Mouse_Brain_Sagittal_Anterior_spatial.tar.gz

Map the reads with alevin

Next we’ll map the reads to the splici mouse index using alevin. Note: Visium uses 16 length cellular barcode (cb) and 12 length UMI, so we can just use the --chromiumV3 flag. The mapping command we will use looks like this:

salmon alevin -l ISR -i refdata-gex-mm10-2020-A/sal_index \
  -1 V1_Mouse_Brain_Sagittal_Anterior_Section_1_fastqs/V1_Mouse_Brain_Sagittal_Anterior_Section_1_S7_L001_R1_001.fastq.gz \
  V1_Mouse_Brain_Sagittal_Anterior_Section_1_fastqs/V1_Mouse_Brain_Sagittal_Anterior_Section_1_S7_L002_R1_001.fastq.gz \
  -2 V1_Mouse_Brain_Sagittal_Anterior_Section_1_fastqs/V1_Mouse_Brain_Sagittal_Anterior_Section_1_S7_L001_R2_001.fastq.gz \
  V1_Mouse_Brain_Sagittal_Anterior_Section_1_fastqs/V1_Mouse_Brain_Sagittal_Anterior_Section_1_S7_L002_R2_001.fastq.gz \
  --chromiumV3 \
  -o alevin_out -p 16 --sketch

Once this is completed, the mapping information that we need will exist in the directory alevin_out.

Quantify with alevin-fry

Next, we’ll follow the familiar steps to quantify the mapped reads with alevin-fry. We’ll use the knee method for generating the permit list and, since we are using USA mode quantifications to appropriately tag reads of intronic origin, we’ll use the cr-like UMI resolution scheme.

alevin-fry generate-permit-list -k -d fw -i alevin_out -o alevin_out/permitlist_knee/ 

alevin-fry collate -i alevin_out/permitlist_knee/ -r alevin_out -t 16 

alevin-fry quant -r cr-like --use-mtx -m mouse_ref/t2g_3col.tsv -i alevin_out/permitlist_knee/ -o alevin_out/quant_knee_cr-like.usa -t 16 

where, above refdata-gex-mm10-2020-A is the path to the directory where you generated the splici reference transcriptome, including the 2 and 3-column transcript to gene mappings. Note : Since we are using the splici transcriptome, we make use here of the 3-column format transcript to gene tsv file (passed via the -m option to quant) to allow us to differentiate between UMIs that likely originate from spliced or unspliced molecules.

At this point, we have mapped and quantified the data. The quantification output will live in alevin_out/quant_knee_cr-like.usa.

Loading data into R

We start by loading the required R packages.

suppressPackageStartupMessages({
    library(devtools)
    library(ggplot2)
    library(SingleCellExperiment)
    library(Seurat)
    library(patchwork)
})

# set the seed 
set.seed(271828)

This is the function we’ll use to load our quantification matrices into R.

#' Read alevin-fry quantifications into a SingleCellExperiment object
load_fry <- function(frydir, which_counts = c('S', 'A'), verbose = FALSE) {
  suppressPackageStartupMessages({
    library(rjson)
    library(Matrix)
    library(SingleCellExperiment)
  })
  
  # read in metadata
  qfile <- file.path(frydir, "quant.json")
  if (!file.exists(qfile)) {
    qfile <- file.path(frydir, "meta_info.json")
  }
  meta_info <- fromJSON(file = qfile)
  ng <- meta_info$num_genes
  usa_mode <- meta_info$usa_mode
  
  if (usa_mode) {
    if (length(which_counts) == 0) {
      stop("Please at least provide one status in 'U' 'S' 'A' ")
    }
    if (verbose) {
      message("processing input in USA mode, will return ", paste(which_counts, collapse = '+'))
    }
  } else if (verbose) {
    message("processing input in standard mode, will return spliced count")
  }

  # read in count matrix
  af_raw <- readMM(file = file.path(frydir, "alevin", "quants_mat.mtx"))
  # if usa mode, each gene gets 3 rows, so the actual number of genes is ng/3
  if (usa_mode) {
    if (ng %% 3 != 0) {
      stop("The number of quantified targets is not a multiple of 3")
    }
    ng <- as.integer(ng/3)
  }
  
  # read in gene name file and cell barcode file
  afg <- read.csv(file.path(frydir, "alevin", "quants_mat_cols.txt"), 
                  strip.white = TRUE, header = FALSE, nrows = ng, 
                  col.names = c("gene_ids"), row.names = 1)
  afc <- read.csv(file.path(frydir, "alevin", "quants_mat_rows.txt"), 
                  strip.white = TRUE, header = FALSE,
                  col.names = c("barcodes"), row.names = 1)

  # if in usa_mode, sum up counts in different status according to which_counts
  if (usa_mode) {
    rd <- list("S" = seq(1, ng), "U" =  seq(ng + 1, 2 * ng),
               "A" =  seq(2 * ng + 1, 3 * ng))
    o <- af_raw[, rd[[which_counts[1]]], drop = FALSE]
    for (wc in which_counts[-1]) {
      o <- o + af_raw[, rd[[wc]], drop = FALSE]
    }
  } else {
    o <- af_raw
  }
  
  # create SingleCellExperiment object
  sce <- SingleCellExperiment(list(counts = t(o)),
                              colData = afc,
                              rowData = afg
  )
  sce
}

Load the quantifications

Next, we’ll load the quantifications into a SingleCellExperiemnt object. The quantifications will be loaded from the USA-mode matrix, and we’ll keep the S (confidently-spliced) and A (ambiguous) UMIs for each gene.

quant_sce <- load_fry('alevin_out/quant_knee_cr-like.usa', verbose=TRUE)
## processing input in USA mode, will return S+A

Because we’ve named each gene with its Entrez ID, but for later analysis we’ll want the gene names, we perform the mapping here (before we put the data into a Seurat object). The gene id to name mapping was obtained directly from our GTF file, but you can find the tsv file that results here.

gid_to_gname <- read.table('mouse_gid_to_gname.tsv')
rownames(quant_sce) <- gid_to_gname$V2[match(rownames(quant_sce), gid_to_gname$V1)]

Now that everything is prepared as we need, we create the Seurat object.

assay <- "Spatial"
brain <- CreateSeuratObject(counts = counts(quant_sce), project = "SPATIAL", assay = assay)
## Warning: Non-unique features (rownames) present in the input matrix, making
## unique
## Warning: The following arguments are not used: row.names
brain
## An object of class Seurat 
## 32085 features across 2907 samples within 1 assay 
## Active assay: Spatial (32085 features, 0 variable features)

Next, we load the spatial image data to extract the spatial 2D coordinates of each cellular barcode, and add the metadata to the Seurat object.

image.data <- Read10X_Image("spatial/")
# Since the names of alevin cb is different from 10x
# we rename the cells and filter the image data
# to have the metadata for only quantified cells
rownames(image.data@coordinates) <- gsub("-1", "", rownames(image.data@coordinates))
common.cells <- intersect(Cells(x = brain), rownames(image.data@coordinates))

image.data <- image.data[common.cells]
brain <- subset(brain, cells = common.cells)

# adding image data to Seurat object
DefaultAssay(object = image.data) <- "Spatial"
brain[['slice']] <- image.data

Visualization

Once we have the Seurat object we can use the exciting features made available through the Seurat package. We follow the workflow of the spatial vignette provided by the Seurat package for visualizing the Spatial data.

plot1 <- VlnPlot(brain, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(brain, features = "nCount_Spatial") + theme(legend.position = "right")
wrap_plots(plot1, plot2)