SPARK识别空间转录组数据显著高变基因R包学习
·
SPARK
SPARK用于在空间转录组学研究中识别具有空间表达模式的基因。其采用广义空间线性模型，直接对各种空间转录组技术生成的计数数据进行建模。它依托惩罚准似然算法实现大规模的可扩展计算，并利用最新开发的统计公式进行假设检验；这不仅能有效控制I类错误（假阳性），同时还具备极高的统计效能。推荐应用场景：样本量小于3,000，且数据稀疏性相对较低。
SPARK分析流程
1.导入
rm(list=ls())
library(SPARK)
library(Seurat)
library(ggplot2)
library(hdf5r)
library(qs)
# 示例数据
load("./Layer2_BC_Count.rds")
2.数据预处理
将编码在列名中的连续型空间坐标解析出来，并关联每个spot的测序深度，形成结构化的样本注释表info，为后续空间定位、质控和可视化奠定基础。
info <- cbind.data.frame(x=as.numeric(sapply(strsplit(colnames(rawcount),split="x"),"[",1)),
y=as.numeric(sapply(strsplit(colnames(rawcount),split="x"),"[",2)),
total_counts=apply(rawcount,2,sum))
rownames(info) <- colnames(rawcount)
3.创建SPARK对象
# 创建用于分析的SPARK对象。此步骤排除低表达的基因。
# 过滤基因和细胞/spots
spark <- CreateSPARKObject(counts=rawcount,
location=info[,1:2],
percentage = 0.1,
min_total_counts = 10)
4.拟合统计模型
# 计算每个细胞/spot的总计数
spark@lib_size <- apply(spark@counts, 2, sum)
# 以前十个基因为例。
# spark@counts <- spark@counts[1:10,]
# 在零假设下拟合统计模型
spark <- spark.vc(spark,
covariates = NULL,
lib_size = spark@lib_size,
num_core = 5,
verbose = F)
5.空间表达模式基因
# 对具有空间表达模式的基因进行检验。默认情况下，核矩阵会根据坐标自动计算，并检查核矩阵的正定性。此外，也支持用户自行提供核矩阵。
# 计算pval
options(mc.cores = 1)
spark <- spark.test(spark,
check_positive = T,
verbose = F)
# 输出最后结果
head(spark@res_mtest[,c("combined_pvalue","adjusted_pvalue")])
# combined_pvalue adjusted_pvalue
# GAPDH 7.477432e-09 4.233461e-06
# MAPKAPK2 1.016092e-01 1.000000e+00
# MCL1 1.149522e-08 6.079079e-06
# TMEM109 4.304043e-01 1.000000e+00
# TMEM189 6.189066e-01 1.000000e+00
# ITPK1 7.213287e-01 1.000000e+00
可以按照p值进行排序，然后选择靠前的基因并根据表达量综合判断即可
SPARK-X
SPARK-X基于一个稳健的协方差检验框架，能够对通过不同技术平台获得的多种空间转录组学数据进行建模。它依托代数层面的创新实现可扩展的高效计算，并结合新开发的统计公式进行假设检验，从而生成校准良好的p值，并具备较高的统计功效。SPARK-X具有极高的计算效率，是目前唯一可扩展应用于HDST 数据的空间表达分析方法。推荐应用场景：样本量大于3,000，无论数据稀疏性结构如何，均表现良好。
SPARK-X分析流程
1.导入
rm(list=ls())
library(SPARK)
library(Seurat)
library(ggplot2)
library(hdf5r)
library(qs)
data <- Read10X_h5("./0-GSM8633896/GSM8633896/filtered_feature_bc_matrix.h5")
dim(data)
object <- CreateSeuratObject(counts = data,
assay = "Spatial",
min.cells = 3, #过滤在少于3个细胞中表达的基因，以减少低表达基因的干扰
project = "GSM8633896")
object
# 再读取
image <- Read10X_Image(image.dir = "./0-GSM8633896/GSM8633896/",
image.name = "tissue_hires_image.png",
filter.matrix = TRUE)
image
dim(image)
image <- image[Cells(x = object)]# 筛选图像对象中包含的SPOT
DefaultAssay(object = image) <- "Spatial" # 设置图像对象的默认assay为"Spatial"
object@images[["GSM8633896"]]@scale.factors
# 添加图片到object中
object[["GSM8633896"]] <- image
object@images[["GSM8633896"]]@scale.factors$lowres <- object@images[["GSM8633896"]]@scale.factors$hires
sce.all <- object
## 简单探索一下数据结构
as.data.frame(sce.all[["Spatial"]]$counts[1:4,1:4])
as.data.frame(LayerData(sce.all, assay = "Spatial", layer = "counts")[1:5,1:5])
head(sce.all@meta.data)
table(sce.all$orig.ident)
Layers(sce.all)
Assays(sce.all)
qsave(sce.all, file="sce.all.qs")
2.数据预处理
## 标准化
#options(future.globals.maxSize= 20*1024*1024^2)
sce.all <- SCTransform(sce.all, assay ="Spatial", verbose = T)
sce.all <- RunPCA(sce.all, assay ="SCT", verbose = FALSE)
sce.all <- FindNeighbors(sce.all, reduction ="pca", dims = 1:30)
sce.all <- FindClusters(sce.all, verbose = FALSE,resolution = 0.3)
sce.all <- RunUMAP(sce.all, reduction ="pca", dims = 1:30)
p1 <- DimPlot(sce.all,
reduction ="umap",
label = TRUE,
label.size = 7);p1
ggsave(filename ="DimPlot.pdf", width = 9,height = 6, plot = p1)
# check
SpatialFeaturePlot(sce.all,
features = c("FOXP2", "CD4"),
pt.size = 3)
sp_count <- GetAssayData(sce.all,layer = "counts")
sp_count[1:5,1:5]
# 5 x 5 sparse Matrix of class "dgCMatrix"
# AAACAAGTATCTCCCA-1 AAACAGAGCGACTCCT-1 AAACATTTCCCGGATT-1 AAACCCGAACGAAATC-1 AAACCGTTCGTCCAGG-1
# AL627309.1 . . . . .
# AL627309.5 . . . . .
# AP006222.2 . . . . .
# LINC01409 . . . . .
# LINC01128 . 1 . . .
# 提取每个样本的注释信息，即位置或坐标。
# 从原始数据中提取坐标。
site <- GetTissueCoordinates(sce.all)
info <- cbind.data.frame(x=site$x,
y=site$y,
total_counts=apply(rawcount,2,sum))
rownames(info) <- colnames(rawcount)
location <- as.matrix(info)
head(info)
# x y total_counts
# AAACAAGTATCTCCCA-1 19800 16897 10596
# AAACAGAGCGACTCCT-1 18399 6025 11989
# AAACATTTCCCGGATT-1 18936 20220 9681
# AAACCCGAACGAAATC-1 22055 15384 11687
# AAACCGTTCGTCCAGG-1 9386 17512 9857
# AAACCTAAGCAGCCGG-1 16507 21431 11800
sp_count[26:30,1:5]
# 5 x 5 sparse Matrix of class "dgCMatrix"
# 1000x100 1000x103 1000x113 1000x114 1000x116
# Gm42418 1 . 2 . .
# Gm10925 . . . . .
# Gm7135 . . . . .
# Atrx . . . . .
# Celf2 . . . . .
# 移除线粒体基因
mt_idx <- grep("mt-",rownames(sp_count))
if(length(mt_idx)!=0){
sp_count <- sp_count[-mt_idx,]
}
3.SPARK-X分析数据
sparkX <- sparkx(sp_count,location,
numCores=1,option="mixture")
head(sparkX$res_mtest)
# combinedPval adjustedPval
# AP006222.2 7.736043e-01 1.000000000
# LINC01409 1.944115e-02 0.263246626
# LINC01128 5.347086e-05 0.001186625
# LINC00115 2.873962e-01 1.000000000
# FAM41C 6.868424e-02 0.848029374
# TMEM53 5.566834e-03 8.294263e-02
同样也能够得到需要差异基因，接下来就可以挑选相应的基因进行研究及可视化。
4.可视化
SpatialFeaturePlot(sce.all,
features = c("TMEM53"),
pt.size = 3)
参考资料
-
SPARK/SPARK-X github：https://xzhoulab.github.io/SPARK/02_SPARK_Example/
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台。更多相关内容可关注公众号：生信方舟 。
更多推荐
8
0- 0
已为社区贡献1条内容
所有评论(0)