SCTransform：单细胞样本的标准化

每次我在国内直接用这种方法下载数据集都没有成功，如上，下载安装包，本地安装：

install.packages(“H:/singlecell/Seurat/stxBrain.SeuratData_0.1.1.tar.gz”, repos = NULL, type = “source”)
brain <- LoadData(“stxBrain”, type = “anterior1”)
当然，我们不推荐这种读取数据的方法，毕竟没有人把我们的数据也打包成一个R包的样子，我们拿到的是10 X Space Ranger的输出结果：

├── analysis
│?? ├── clustering
│?? ├── diffexp
│?? ├── pca
│?? ├── tsne
│?? └── umap
├── cloupe.cloupe
├── filtered_feature_bc_matrix
│?? ├── barcodes.tsv.gz
│?? ├── features.tsv.gz
│?? └── matrix.mtx.gz
├── filtered_feature_bc_matrix.h5
├── metrics_summary.csv
├── molecule_info.h5
├── possorted_genome_bam.bam
├── possorted_genome_bam.bam.bai
├── raw_feature_bc_matrix
│?? ├── barcodes.tsv.gz
│?? ├── features.tsv.gz
│?? └── matrix.mtx.gz
├── raw_feature_bc_matrix.h5
├── spatial # 空间信息全在这 :这些文件是用户提供的原始全分辨率brightfield图像的下采样版本。
下采样是通过box滤波实现的，它对全分辨率图像中像素块的RGB值进行平均，得到下采样图像中一个像素点的RGB值。
│?? ├── aligned_fiducials.jpg 这个图像的尺寸是tissue_hires_image.png。由基准对齐算法发现的基准点用红色高亮显示。此文件对于验证基准对齐是否成功非常有用。
│?? ├── detected_tissue_image.jpg
│?? ├── scalefactors_json.json
│?? ├── tissue_hires_image.png 图像的最大尺寸为2,000像素
│?? ├── tissue_lowres_image.png 图像的最大尺寸为600像素。
│?? └── tissue_positions_list.csv
└── web_summary.html
其实只要 Space Ranger的输出结果 filtered_feature_bc_matrix.h5文件和一个spatial文件夹就可以了，一如？stxBrain中给出的示例：

setwd(“H:\singlecell\spaceranger\V1_Mouse_Brain_Sagittal_Anterio/”)

Load the expression data

expr.url <- ‘H:\singlecell\spaceranger\V1_Mouse_Brain_Sagittal_Anterio/V1_Mouse_Brain_Sagittal_Anterior_filtered_feature_bc_matrix.h5’
expr.data <- Seurat::Read10X_h5(filename = expr.url )
anterior1 <- Seurat::CreateSeuratObject(counts = expr.data, project = ‘anterior1’, assay = ‘Spatial’)
anterior1$slice<-1anterior1$region <- ‘anterior’

Load the image data

img.url <- ‘H:\singlecell\spaceranger\V1_Mouse_Brain_Sagittal_Anterio/V1_Mouse_Brain_Sagittal_Anterior_spatial.tar.gz’
untar(tarfile = img.url)
img <- Seurat::Read10X_Image(image.dir = ‘H:\singlecell\spaceranger\V1_Mouse_Brain_Sagittal_Anterio/spatial’)
Seurat::DefaultAssay(object = img) <- ‘Spatial’
img <- img[colnames(x = anterior1)]
anterior1[[‘image’]] <- img

anterior1

An object of class Seurat
31053 features across 2696 samples within 1 assay
Active assay: Spatial (31053 features)

}
如果这两个文件在同一个文件下的话，也可以这样读入：

list.files(“H:\singlecell\spaceranger\V1_Mouse_Brain_Sagittal_Anterio”) # 注意命名要正确

[1] “filtered_feature_bc_matrix.h5” “spatial”

[3] “V1_Mouse_Brain_Sagittal_Anterior_filtered_feature_bc_matrix.h5” “V1_Mouse_Brain_Sagittal_Anterior_spatial.tar.gz”

brain<-Seurat::Load10X_Spatial(“H:\singlecell\spaceranger\V1_Mouse_Brain_Sagittal_Anterio”)
？Load10X_Spatial
brain

An object of class Seurat
31053 features across 2696 samples within 1 assay
Active assay: Spatial (31053 features)
空间数据如何存储在Seurat中?
来自10x的visium数据包括以下数据类型:

通过基因表达矩阵得到一个点(spot )
组织切片图像(采集数据时H&E染色)
用于显示的原始高分辨率图像与低分辨率图像之间的比例因子。
在Seurat对象中，spot by基因表达矩阵与典型的“RNA”分析类似，但包含spot水平，而不是单细胞水平的数据。图像本身存储在Seurat对象中的一个images 槽（slot ）中。图像槽还存储必要的信息，以将斑点与其在组织图像上的物理位置相关联。

数据预处理
在spot中基因表达数据进行初始的预处理步骤与典型的scRNA-seq相似。我们首先需要对数据进行规范化，以考虑数据点之间测序深度的差异。我们注意到，对于空间数据集来说，分子数/点的差异可能是巨大的，特别是在整个组织的细胞密度存在差异的情况下。我们在这里看到大量的异质性，这需要有效的标准化。

熟悉的函数名~

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

这些图表明，分子计数（molecular counts）在点间的差异不仅是技术性的，而且还依赖于组织解剖。例如，组织中神经元消耗的区域(如皮层白质)，可再生地显示出较低的分子计数。因此，标准方法(如LogNormalize函数)可能会有问题，因为它会强制每个数据点在标准化之后具有相同的底层“大小”。

作为一种替代方法，我们推荐使用sctransform (Hafemeister和Satija，已出版)，它构建了基因表达的正则化负二项模型，以便在保留生物差异的同时考虑技术因素。有关sctransform的更多信息，请参见 here的预印和here的Seurat教程。sctransform将数据归一化，检测高方差特征，并将数据存储在SCT分析中。

brain <- SCTransform(brain, assay = “Spatial”, return.only.var.genes = FALSE, verbose = FALSE)
sct 与log-归一化相比，结果如何?
为了探究规范化方法的不同，我们研究了sctransform和log规范化结果如何与UMIs的数量相关。为了进行比较，我们首先重新运行sctransform来存储所有基因的值(这将会比较慢)，并通过NormalizeData运行一个log规范化过程。

作者：周运来就是我
链接：https://www.jianshu.com/p/f6da86489784

Seurat从3.0版本引进了SCTransform这个函数用来对数据做标准化，并且这一个函数可以代替三个函数（NormalizeData, ScaleData, FindVariableFeatures）的运行。 且其对测序深度的校正效果要好于log标准化。（10万以内的细胞都建议使用SCT标准化）
️SCTransform对测序深度的校正效果很好，也可用于矫正线粒体等因素的影响，但不能用于批次矫正。

sctransform方法学文章
官网：
https://cran.r-project.org/web/packages/sctransform/index.html
https://satijalab.org/seurat/archive/v3.0/sctransform_vignette.html

官网对这个函数的描述：A normalization method for single-cell UMI count data using a variance stabilizing transformation. The transformation is based on a negative binomial regression model with regularized parameters. As part of the same regression framework, this package also provides functions for batch correction, and data correction. See Hafemeister and Satija 2019 doi:10.1186/s13059-019-1874-1 for more details. （这是一个用方差稳定变换对单细胞UMI count 数据标准化的方法，方差稳定变换是基于负二项回归。这个函数在对数据进行均一化的同时还可以去除线粒体红细胞等混杂因素的影响。）

1. NormalizeData, ScaleData, FindVariableFeatures和SCTransform结果的存放位置
   示例数据下载：pbmc3k

#读入数据，创建seurat对象
pbmc <- Read10X(‘./filtered_gene_bc_matrices/hg19/’)
pbmc <- CreateSeuratObject(pbmc,project = ‘pbmc3k’,min.cells = 3,min.features = 200)

质控

pbmc[[“percent.mt”]] <- PercentageFeatureSet(pbmc, pattern = “^MT-”)
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
常规NormalizeData, ScaleData, FindVariableFeatures操作

pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc, selection.method = “vst”, nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
View(pbmc)
这一步的结果都储存在pbmc@assays[[“RNA”]]中。

再使用SCTransform进行标准化

pbmc <- SCTransform(pbmc, vars.to.regress = “percent.mt”, verbose = FALSE)
View(pbmc)
再次查看pbmc的数据结构，新出现了一个pbmc@assays[[“SCT”]]的槽，并将SCTransform的结果储存在中，结果中也有对应的原始矩阵，Normalize矩阵，scale后的矩阵和找到的var.features（SCTransform默认是找3000个var.features，上面的方法 默认是2000个）。
在常规分析中，使用少量的PC既能关注到关键的生物学差异，又能够不引入更多的技术差异，相当于一种保守性的做法。它会失去一些生物差异信息但是同时又在常规手段中比较安全。
但sctransform的归一化、标准化都做得不错，多输入一些PCs能提取更多的生物差异，并且兼顾不引入技术误差。sctransform认为：新增加的这1000个基因就包含了之前没有检测到的微弱的生物学差异。而且，即使使用全部的全部的基因去做下游分析，得到的结果也是和sctransform这3000个基因的结果相似

因此NormalizeData, ScaleData, FindVariableFeatures和SCTransform的结果存放在不同的slot中，不必担心进行了Normalize和scale再进行SCTransform会对数据过度处理，在后续处理数据的时候注意选择需要的矩阵即可。

**

SCTransform得到的结果都保存在pbmc[[“SCT”]]中： 值得注意的是，SCTransform分析结果是@scale.data，@counts和@data都是用@scale.data反向计算出来的。

**

pbmc[[“SCT”]]@scale.data ：包含了残差数据，用作PCA的输入。这个数据不是稀疏矩阵，因此会占用大量内存。不过SCTransform函数计算的时候，为了节省内存，默认使用了return.only.var.genes = TRUE ，只对差异基因进行scale的结果。设置为FALSE相当于三步法的ScaleData(features=all.genes)。
关于为什么SCT运行完基因数会减少

Christoph作了回复：sctransform的min_cells参数会忽略一部分细胞。如果想要所有细胞，调整参数即可。

pbmc[[“SCT”]]@counts ：相当于细胞数据量相同情况下的值（也就是均一化了的UMI count值）。
pbmc[[“SCT”]]@data：包含了上面count值的log-normalized结果，有利于后面可视化
目前可以使用pbmc[[“SCT”]]@data结果进行差异分析和数据整合，实际上，官方更推荐直接使用残差值pbmc[[“SCT”]]@scale.data。 但是seurat目前不支持，因此还得用@data的数据。
⚠️在运行完三步法标准化再运行SCT标准化，后续分析默认使用SCT矩阵。可以使用DefalutAssay函数可以查看目前默认用的是哪个矩阵。

2. 两种方法比较
   library(VennDiagram)
   library(gplots)
   v1 <- pbmc@assays[[‘RNA’]]@var.features
   v2 <- pbmc@assays[[‘SCT’]]@var.features
   p <- venn.diagram(list(v1,v2), col=c(“red”,“green”), alpha = c(0.5, 0.5),
   category.names=c(‘pbmc’,‘pbmc.sc’), filename=NULL, margin=0.15)
   grid.draw§

可以看到两种方法找到的高变基因还是有一些差别

3. SCTransform做了些什么
   SCTransform利用了正则化负二项分布（regularized negative binomial regression）计算了技术噪音模型，得到的残差是归一化值，有正有负。正值表示：考虑到细胞群体中基因的平均表达量和细胞测序深度，某个细胞的某个基因所包含的UMIs比预测值要高。

首先提取出表达矩阵并对每个基因做简单的统计包括表达量均值，方差，以及可检测率（detection rate）。从一张均值和方差的散点图可以看到在均值较小的一定范围内，均值和方差有线性关系。均值再变大时，方差呈现过度分散（overdispersion）。

作者：Hayley笔记
链接：https://www.jianshu.com/p/21b2f08652ed