3.3 Vignette for creating pseudo-bulk counts matrix

3.3 Vignette for creating pseudo-bulk counts matrix#

In order to perform isoform switching analysis downstream, we would create psuedo samples by breaking/sub-sampling from each cluster into 3 to create replicates. By doing so, each cluster becomes ‘condition’ and each of the sub-clusters becomes a ‘replicate’

For easy sub-sampling, we will annotate sample names as ‘Sample_R+cluster.id+replicate{1-3}’, for ex: ‘Sample_R01’ corresponds to a sub-cluster within cluster index 0.

We follow the pseudo-bulk vignette here closely:

Pseudo-bulk ref : https://hbctraining.github.io/scRNA-seq_online/lessons/pseudobulk_DESeq2_scrnaseq.html

Other Refs:

Seurat Cheatsheet: https://hbctraining.github.io/scRNA-seq_online/lessons/seurat_cheatsheet.html
Seurat commands: https://satijalab.org/seurat/articles/essential_commands.html#fetchdata

Installing BiocStyle:

#{r}
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")

BiocManager::install("BiocStyle")

#{r}
remotes::install_github("cvarrichio/Matrix.utils")

loading required libraries:

#{r}
library(tidyverse)
library(DESeq2)
library(Seurat)
library(SingleCellExperiment)
library(BiocStyle)
library(data.table)
library(Matrix.utils)
knitr::opts_chunk$set(warning=FALSE, error=FALSE, message=FALSE)
library(dplyr)

Part1: Reading and Exploring Seurat object:#

kinnex_sc/complete_dataset/PBMC_BioIVT_10x3p_complete_refGuided.genes-sc_matrix_from_isoquant/PBMC_complete_scKinnex.genes-seurat_obj.rds

#{r}
seurat_readObj <- readRDS('kinnex_sc/complete_dataset/PBMC_BioIVT_10x3p_complete_refGuided.genes-sc_matrix_from_isoquant/PBMC_complete_scKinnex.genes-seurat_obj.rds')
seurat_readObj
length(unique(rownames(seurat_readObj@meta.data)))

Terminal Output:

An object of class Seurat: 30015 features across 12850 samples within 1 assay Active assay: RNA (30015 features, 2000 variable features) 3 layers present: counts, data, scale.data 2 dimensional reductions calculated: pca, umap

12850

#{r}
#library(dplyr)
Features(seurat_readObj) %>% head()
rownames(seurat_readObj) %>% head()
all(rownames(seurat_readObj@meta.data) == Cells(seurat_readObj))
Idents(seurat_readObj) %>% head()

Terminal Out:

[1] "novel-gene-chr10-13364^novel-gene-chr10-13364" "novel-gene-chr10-2287^novel-gene-chr10-2287"
[3] "novel-gene-chr10-26782^novel-gene-chr10-26782" "novel-gene-chr10-29791^novel-gene-chr10-29791"
[5] "novel-gene-chr10-30176^novel-gene-chr10-30176" "novel-gene-chr10-31726^novel-gene-chr10-31726"
[1] "novel-gene-chr10-13364^novel-gene-chr10-13364" "novel-gene-chr10-2287^novel-gene-chr10-2287"
[3] "novel-gene-chr10-26782^novel-gene-chr10-26782" "novel-gene-chr10-29791^novel-gene-chr10-29791"
[5] "novel-gene-chr10-30176^novel-gene-chr10-30176" "novel-gene-chr10-31726^novel-gene-chr10-31726"
[1] TRUE
AAACAACGAAAGAATC AAACAACGACAGTCTA AAACAACGAGTTAGAA AAACACCTGCTTCCAC AAACACCTGGAGGAGG
           0                3                0                1                1
AAACACCTGGTACCTA
           2
Levels: 0 1 2 3 4 5 6 7 8 9 10 11

Part2: Dividing clusters into sub-clusters#

seperating cluster info

#{r}
seurat_readObj@meta.data
seurat_clusters_info <- FetchData(object = seurat_readObj, vars = c("seurat_clusters"), layer = "meta.data")

randomly sampling indexes from each cluster:

#{r}
seurat_clusters_info$cell_bc <- rownames(seurat_clusters_info)

Dividing cluster into 3 subsets random sampling:

#{r}
for (id in unique(seurat_readObj$seurat_clusters)) {
    cluster_name <- paste("cluster_",id, sep="")
    print(cluster_name)
    df_temp <- seurat_clusters_info[seurat_clusters_info$seurat_clusters==id,]
    df_temp$sample_id <- sample(factor(rep(1:3, length.out=nrow(seurat_clusters_info[seurat_clusters_info$seurat_clusters==id,])),
                      labels=paste0("sample_R",id,1:3)))
    assign(cluster_name,df_temp)
}

adding samples ids to metadata objects:

#{r}
#summary(cluster_0$sample_id)
#summary(cluster_2$sample_id)
#summary(cluster_10$sample_id)
#summary(cluster_11$sample_id)
ss <- seurat_readObj
w_sample_ids <- rbind(cluster_0, cluster_1,cluster_2,cluster_3, cluster_4, cluster_5, cluster_6,
                  cluster_7,cluster_8,cluster_9,cluster_10,cluster_11)

#ss@meta.data
#w_sample_ids

Adding sample names in Sample:

#{r}
ss@meta.data <- merge(ss@meta.data, dplyr::select(w_sample_ids, sample_id), by=0, all=TRUE)
ss@meta.data

Assigning back to Seurat Object:

#{r}
seurat_readObj <- ss

Extracting metadata:

#{r}
metadata <- seurat_readObj@meta.data

Extract raw counts and metadata to create SingleCellExperiment object

#{r}
counts <- seurat_readObj@assays$RNA$counts

Set up metadata as desired for aggregation and DE analysis

#{r}
metadata$cluster_id <- factor(seurat_readObj@active.ident)

Part3 - Create single cell experiment object#

#{r}
sce <- SingleCellExperiment(assays = list(counts = counts),
                       colData = metadata)

Exploring the raw counts for the dataset Checking the assays present

#{r}
assays(sce)

Terminal Out:

List of length 1 names(1): counts

Check the counts matrix

#{r}
dim(counts(sce))
counts(sce)[1:6, 1:6]

Part4: Preparing the single-cell dataset for pseudobulk analysis#

Extracting necessary metrics for aggregation by cell type in a sample:

#{r}
# Extract unique names of clusters (= levels of cluster_id factor variable)
cluster_names <- levels(colData(sce)$cluster_id)
cluster_names

# Total number of clusters
length(cluster_names)

Number of cells in each cluster:

#{r}
for (i in cluster_names) {
 print(paste(i, length(colData(sce)$cluster_id[colData(sce)$cluster_id==i]), sep = ":"))
}

Terminal Out:

“0:2314” “1:2123” “2:1940” “3:1822” “4:1604” “5:1062” “6:1035” “7:281” “8:280” “9:235” “10:101” “11:53”

#{r}
# Extract unique names of samples (= levels of sample_id factor variable)
sample_names <- levels(colData(sce)$sample_id)
sample_names

# Total number of samples
length(sample_names)

Part5: Subset metadata#

Subset metadata to include only the variables you want to aggregate across (here, we want to aggregate by sample and by cluster)

#{r}
#colData(sce)
groups <- colData(sce)[, c("cluster_id", "sample_id")]
head(groups)

Aggregate across cluster-sample groups - transposing row/columns to have cell_ids as row names matching those of groups

#{r}
aggr_counts <- aggregate.Matrix(t(counts(sce)),
                            groupings = groups, fun = "sum")

Exploring aggregated output matrix

#{r}
class(aggr_counts)
dim(aggr_counts)
aggr_counts[1:6, 1:6]

Transpose aggregated matrix to have genes as rows and samples as columns

#{r}
aggr_counts <- t(aggr_counts)
aggr_counts[1:6, 1:6]

Understanding tstrsplit()

#{r}
## Exploring structure of function output (list)
tstrsplit(colnames(aggr_counts), "_") %>% str()

## Comparing the first 10 elements of our input and output strings
head(colnames(aggr_counts), n = 10)
head(tstrsplit(colnames(aggr_counts), "_")[[1]], n = 10)

aggr_counts

#{r}
# As a reminder, we stored our cell types in a vector called cluster_names
cluster_names


# Loop over all cell types to extract corresponding counts, and store information in a list

## Initiate empty list
counts_ls <- list()

for (i in 1:length(cluster_names)) {

    ## Extract indexes of columns in the global matrix that match a given cluster
    column_idx <- which(tstrsplit(colnames(aggr_counts), "_")[[1]] == cluster_names[i])

    ## Store corresponding sub-matrix as one element of a list
    counts_ls[[i]] <- aggr_counts[, column_idx]
    names(counts_ls)[i] <- cluster_names[i]

}

# Explore the different components of the list
str(counts_ls)

Part6: Generating matching metadata at the sample-level#

#{r}
# Reminder: explore structure of metadata
head(colData(sce))

# Extract sample-level variables
metadata <- colData(sce) %>%
as.data.frame() %>%
dplyr::select(seurat_clusters,sample_id)

dim(metadata)
head(metadata)

# Exclude duplicated rows
metadata <- metadata[!duplicated(metadata), ]

dim(metadata)
head(metadata)

Rename rows:

#{r}
rownames(metadata) <- metadata$sample_id
head(metadata)

Number of cells per sample and cluster

#{r}
t <- table(colData(sce)$sample_id,
       colData(sce)$cluster_id)
t

#{r}
temp <- '11_sample_R113'
tstrsplit(temp, "_")[[1]]
paste(tstrsplit(temp, "_")[[2]],tstrsplit(temp, "_")[[3]],sep='_')

Creating metadata list

#{r}
## Initiate empty list
metadata_ls <- list()

for (i in 1:length(counts_ls)) {

    ## Initiate a data frame for cluster i with one row per sample (matching column names in the counts matrix)
    df <- data.frame(cluster_sample_id = colnames(counts_ls[[i]]))
    head(df)
    ## Use tstrsplit() to separate cluster (cell type) and sample IDs
    df$cluster_id <- tstrsplit(df$cluster_sample_id, "_")[[1]]
    df$sample_id  <- paste(tstrsplit(temp, "_")[[2]],tstrsplit(temp, "_")[[3]],sep='_')


    ## Retrieve cell count information for this cluster from global cell count table
    idx <- which(colnames(t) == unique(df$cluster_id))
    cell_counts <- t[, idx]
    ## Remove samples with zero cell contributing to the cluster
    cell_counts <- cell_counts[cell_counts > 0]

    ## Match order of cell_counts and sample_ids
    sample_order <- match(df$sample_id, names(cell_counts))
    cell_counts <- cell_counts[sample_order]

    ## Append cell_counts to data frame
    df$cell_count <- cell_counts


    ## Join data frame (capturing metadata specific to cluster) to generic metadata
    df <- plyr::join(df, metadata,
                 by = intersect(names(df), names(metadata)))

    ## Update rownames of metadata to match colnames of count matrix, as needed later for DE
    rownames(df) <- df$cluster_sample_id

    ## Store complete metadata for cluster i in list
    metadata_ls[[i]] <- df
    names(metadata_ls)[i] <- unique(df$cluster_id)
    }

# Explore the different components of the list
str(metadata_ls)

we have matching lists of counts matrices and sample-level metadata for each cell type, and we are ready to proceed with pseudobulk differential expression analysis.

#{r}
# Double-check that both lists have same names
all(names(counts_ls) == names(metadata_ls))
#counts_ls$`0`
#counts_ls[[idx]]

In absence of ‘group_id’, we can assign cluster names as groups

#{r}
colnames(counts_ls[[1]])

merging the matrices - one for each cluster - corresponding to counts for 3 replicates - to get gene counts:

#{r}
merged_sm <- RowMergeSparseMatrices(counts_ls[[1]],counts_ls[[2]])

for (i in 3:length(counts_ls)) {
    print(i)
    print(colnames(counts_ls[[i]]))
    merged_sm <- RowMergeSparseMatrices(merged_sm, counts_ls[[i]])
}

colnames(merged_sm)

Writing combined counts to tsv: Note - the tedious code below, which can use an R proficient R, worksaround the structure of the dgCsparse matrix object to assign rownames as ‘isoquant_id’ to thee final counts table.

#{r}
#head(merged_sm, 3)
write.table(as.matrix(merged_sm),
        file ="kinnex_sc/complete_dataset/pseudo_bulk_counts.tsv",
        row.names=TRUE,
        sep="\t")

# temp <-
# read.table(file ="kinnex_sc/complete_dataset/pseudo_bulk_counts.tsv",
#            sep="\t")
#
#
# temp$isoform_id <- rownames(temp)
# head(temp)
#
# write.table(temp,
#             file ="kinnex_sc/complete_dataset/pseudo_bulk_counts.tsv",col.names = TRUE,
#             sep="\t")

3.3 Vignette for creating pseudo-bulk counts matrix

Contents

3.3 Vignette for creating pseudo-bulk counts matrix#

Part1: Reading and Exploring Seurat object:#

Part2: Dividing clusters into sub-clusters#

Part3 - Create single cell experiment object#

Part4: Preparing the single-cell dataset for pseudobulk analysis#

Part5: Subset metadata#

Part6: Generating matching metadata at the sample-level#