Seurat analysis report of BD Rhapsody TAS-Seq WTA data

Demultiplexing cells by cell hashing

For assignment of each tags to each cell barcodes, read counts of each tag in each valid cell barcode, which defined by the cDNA matrix, were extracted from tag/cell barcode expression matrix. Unassigned cell barcodes were labeled as “not-detected” cells. Then, sum of the total read counts of each tags were normalized to the minimum sum count of each tags, and log2 fold-change between first most tag counts and second most tag counts within each cell barcode. Each cell barcode were ranked by the fold-change ascending order, and top 3.28% cells were identified as doublets, which were theoretically detectable doublets calculated by the Poisson’s distribution based on the number of loaded cells, total Rhapsody well number, and number of tags used.

Finally, remained cell barcodes were assigned to the first most counted tags. Each tag expression in each cell barcode were log2(x+1)-transformed, z-scaled by each cell barcode, and visualized by using pheatmap package in R version 4.1.2 (2021-11-01).

Figure 1: Heatmap representation of tag expression and assignment in each cell barcode

File locations

Demultiplexed plot was stored at ./result/UTWTA_results/plots/UTWTA_demultiplex.png. Demultiplexed expression matrix was stored at ./result/UTWTA_results/matrix/matrix_inflection_demulti_UTWTA.txt.gz.

Perform DBEC correction (mRNA leak-derived backgroud subtraction of gene expression) For DBEC correction, genes of which log2(x+1)-transformed maximum expression over 8 were selected, and biexponential transformation was performed by using FlowTrans package in R version 4.1.2 (2021-11-01). Next, normal mixture components (model E) was detected by using mclust package in R version 4.1.2 (2021-11-01), average expression of each component were calculated, and genes of which the maximum average expression over 5.5 were selected for background subtration. Then, if the components of which difference of the maximum average expression against maximal average expression under 5 were detected, the expression value of the components set as 0 for background subtracion.

File locations
DBEC threshold plots were stored at ./result/UTWTA_results/plots/DBEC_thresholds_UTWTA Demultiplexed/DBEC-corrected expression matrix was stored at ./result/UTWTA_results/matrix/matrix_inflection_demulti_DBEC_UTWTA.txt.gz

Figure 3: Heatmap representation of tag expression and assignment in each cell barcode (after DBEC correction)

Seurat analysis

Perform Seurat (version 4.0.4) analysis.
- https://satijalab.org/seurat/
Please see Seurat vignettes (PBMC3K tutorial) for detail explanation of the each step of Seurat analysis.
- https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
Entire Seurat object was saved by using qs package with .qs extension. Please use qs package to read .qs file with R.

We first calculating metrics (percent.mito, percent.ribosomal protein, percent.ribosomal RNA) and filtered out mitochondrial gene-high cells (over 0.25), doublets, and not-tag-assigned cells.
Because BD Rhapsody TAS-Seq data is non-UMI data (raw count data), we perform global normalization to 1M tags (according to the Tabula Muris (Nature 2019) parameter settings) and perform scaling by using Seurat ScaleData function with regressing out library size (total raw read count of each cell) as a confounding factor.

PCA and Jackstraw analysis

Next, PCA analysis was performed against 7618 of highly-variable genes identified by FindVariableFeatures (selection.method=mvp, mean.cutoff=c(0.1, Inf), dispersion.cutoff=c(0.5, Inf)) function in Seurat v4.0.4 in R version 4.1.2 (2021-11-01).

Figure 3-1. Jackstraw plot

1:71 PCs were selected (p<1e-05).

Jackstraw plot file was stored at ./result/UTWTA_results/Seurat/Seurat_plots/UTWTA_Jackstraw.png.

Clustering

Cell clustering was performed by using FindNeighbors and FindClusters function in Seurat v4.0.4 in R version 4.1.2 (2021-11-01).
1:71 PCs were used for clustering analysis.
Clustering resolution parameters were changed sequentially and evaluate clustering stability by using silhouette score against FIt-SNE space.

Seurat clustering silhouette score plot

We assumed that resolution 1.5 was optimal because mean silhouette score of resolution 1.5 was maximal among calculated resolutions.

Figure 4. Boxplot of silhouette scores of clustering in each resolution parameter.