Data normalization using scTransform

Rationale

1. Per gene glm model parameters are learned. This can be done on a subset of genes/cells to make the run faster.
2. sctransform gives user an option to choose methods. For example, if set to ‘poisson’ regression is done; negative binomial theta parameter is estimated using the response residual. If method is set to ‘glmGamPo’, coefficient and theta are estimated by a call to glmGamPoi::glm_gp.
3. Regularize parameters, both offset and slope.
4. Transform UMI counts into person residuals.

# install required pacakges ;  glmGamPoi, sctransform,Seurat
#if (!requireNamespace("BiocManager", quietly = TRUE)) install.package ("BiocManager")
#BiocManager::install("glmGamPoi")
#install.packages("sctransform")
#install.packages("Seurat")
library(Seurat)

x10.dat <- Read10X("filtered_feature_bc_matrix")
x10 <- CreateSeuratObject(counts = x10.dat, project = "10x", min.cells = 3, min.features = 200)
dim(x10)
x10.df <- GetAssayData(x10,slot = "count")

# gene attributes
gene_attr <- data.frame(mean = rowMeans(x10.df),detection_rate = rowMeans(x10.df > 0), var=apply(x10.df,1,var))
gene_attr$log_mean <- log10(gene_attr$mean)
gene_attr$log_var <- log10(gene_attr$var)
rownames(gene_attr) <- rownames(x10.df)

# cell attributes 
cell_attr <- data.frame(umi_count = colSums(x10.df), n_gene = colSums(x10.df > 0))
rownames(cell_attr) <- colnames(x10.df)

head(gene_attr)

mean detection_rate          var   log_mean     log_var
Xkr4    0.0672385169   0.0536801328 9.760892e-02 -1.1723819 -1.01051048
Rp1     2.3583287216   0.4964028777 1.100220e+01  0.3726043  1.04147971
Sox17   0.0503597122   0.0058107360 8.942249e-01 -1.2979168 -0.04855325
Mrpl15  0.3303818484   0.2576092972 4.886588e-01 -0.4809838 -0.31099432
Lypla1  0.2814056447   0.2257885999 3.755357e-01 -0.5506672 -0.42534878
Gm37988 0.0008301051   0.0008301051 8.296456e-04 -3.0808669 -3.08110737

head(cell_attr)

                  umi_count n_gene
AAACCCAAGGGATGTC-1      3137   1436
AAACCCACATCCGGTG-1      3132   1432
AAACCCATCAAGCGTT-1      5195   2530
AAACCCATCAGAGTTC-1      2188   1172
AAACGAAAGACCATGG-1       548    358
AAACGAAAGATTAGCA-1      2287   1128

#plot mean vs variance
plot(gene_attr$log_mean, gene_attr$log_var)
abline(lm(log_var ~ log_mean, data = gene_attr), col="red")

# genes with non-zero expression in cells.
genes_cell_count <- rowSums(x10.df >= 0.01)  
umi <- x10.df

min_cells = 10
genes <- rownames(umi)[genes_cell_count >= min_cells]
umi <- umi[genes, ]
dim(umi)

[1] 14799  3614

latent_var = "umi_count"
model_str <- paste0('y ~ ', paste(latent_var, collapse = ' + '))
mat_data <- model.matrix(as.formula(gsub('^y', '', model_str)), cell_attr)


fitted_mod <- glmGamPoi::glm_gp(as.matrix(umi), design = mat_data, size_factors = FALSE)
attributes(fitted_mod)

$names
 [1] "Beta"                          "overdispersions"              
 [3] "overdispersion_shrinkage_list" "deviances"                    
 [5] "Mu"                            "size_factors"                 
 [7] "Offset"                        "ridge_penalty"                
 [9] "data"                          "model_matrix"                 

$class
[1] "glmGamPoi

Variance stabilizing transformation(seurat::vst)

1. Uses a regularized negative binomial regression model.
2. Returns corrected UMI matrix counts and pearson residual.
3. Regress out the unwanted effects of UMI counts.
4. Usage :

vst( umi, cell_attr = NULL, latent_var = c(“log_umi”), batch_var = NULL, latent_var_nonreg = NULL, n_genes = 2000, n_cells = NULL, method = “poisson”, do_regularize = TRUE, theta_regularization = “od_factor”, res_clip_range = c(-sqrt(ncol(umi)), sqrt(ncol(umi))), bin_size = 500, min_cells = 5, residual_type = “pearson”, vst 23 return_cell_attr = FALSE, return_gene_attr = TRUE, return_corrected_umi = FALSE, min_variance = -Inf, bw_adjust = 3, gmean_eps = 1, theta_estimation_fun = “theta.ml”, theta_given = NULL, exclude_poisson = FALSE, use_geometric_mean = TRUE, use_geometric_mean_offset = FALSE, fix_intercept = FALSE, fix_slope = FALSE, scale_factor = NA, vst.flavor = NULL, verbosity = 2, verbose = NULL, show_progress = NULL )

Get Ashwani Kumar’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

Returns:

1. umi_corrected: Matrix of corrected UMI counts.
2. model_str : Representation of model formula.
3. model_pars: Matrix of estimated model parameters per gene.
4. cell_attr: Data frame of cell meta data.
5. gene_attr: Data frame of gene attributes.

#Apply VST 
x10_vst <- sctransform::vst(x10.df,latent_var = c("log_umi"),return_cell_attr=TRUE, return_gene_attr=TRUE, return_corrected_umi=TRUE, verbosity = 1)

Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 15965 by 3614
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 3614 cells
Found 127 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 15965 genes
Computing corrected count matrix for 15965 genes
Calculating gene attributes
Wall clock passed: Time difference of 1.240291 mins

#Plot model parameters
sctransform::plot_model_pars(x10_vst,show_theta=TRUE)

#Inspect model using plot model function
sctransform::plot_model(x10_vst,x10.df,c('Malat1'), plot_residual = TRUE, show_fit = TRUE, show_nr = TRUE)

Reference:

1. Hafemeister, C., Satija, R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol 20, 296 (2019). https://doi.org/10.1186/s13059-019-1874-1
2. sctransform: Variance Stabilizing Transformations for Single Cell UMI Data (r-project.org)