空间转录组|数据读入,标准数据形式外,还有哪些"天残地缺"可以读取

简介: 空间转录组|数据读入,标准数据形式外,还有哪些"天残地缺"可以读取

本文首发于“生信补给站”公众号  https://mp.weixin.qq.com/s/9DJGGU1BBOFO5GMAl9ypFA


空间转录组测序可以同时获得细胞的空间位置信息和基因表达数据,虽然囿于当前单个spot的精度问题,但是在组织细胞功能,肿瘤生物学、发育过程等需要空间位置的研究领域仍然可以提供很多非常有价值的东西。

本节会在spaceranger处理fastq,Seurat处理标准数据(h5 + spatial文件夹)外,额外提供以下几种情况如何处理

(1)raw/filtered_feature_bc_matrix 的三个文件 + spatial 文件夹

(2)矩阵文件 + spatial 文件夹

(3)spatial文件夹如果只提供了tissue_hires_image.png

(4)未提供spatial文件夹,只提供了位置信息

数据来源于2022年CELL 文章 A human fetal lung cell atlas uncovers proximal-distal gradients of differentiation and key regulators of epithelial fates 中的空转Fastq数据(E-MTAB-11265)。

一 fastq数据,spaceranger分析


如若初始是fastq数据,首先在10X官网(https://www.10xgenomics.com/cn/support/software/space-ranger/downloads)下载spaceranger软件然后安装,仍然是使用count参数处理fastq 数据即可。


/path/spaceranger count --id=6332STDY10289523 \
--transcriptome=/path/refdata-gex-GRCh38-2020-A/ \
--fastqs=/path/ST_data/ \
--sample=6332STDY10289523 \
--image=/path/V10S24-031_D1.jpg \
--slide=V10S24-031 \
--area=D1

注意slide 和 area的参数 ,可能来源于图片的命名(本示例中的V10S24-031_D1.jpg) ,可能来自于method 或者code的链接中。

如果没有slide信息的话,可以设置为--unknown-slide

二 Seurat 数据读入 


上面spaceranger count后会得到outs文件夹,里面有网页版的质控报告,loupe文件,spatial文件夹 和 h5文件等,用于后续分析。

1, Load10X_Spatial 最优输入

outs文件夹中的 filtered_feature_bc_matrix.h5 + spatial文件夹 ,注意这两部分结果需要在同一级。


library(Seurat)
library(hdf5r)
spe = Load10X_Spatial(data.dir = "./outs/",
                            filename = "filtered_feature_bc_matrix.h5",
                            assay = "Spatial", 
                            slice = "Lung")
head(spe@meta.data,2)
#                      orig.ident nCount_Spatial nFeature_Spatial
#AAACAAGTATCTCCCA-1 SeuratProject           7232             2774
#AAACACCAATAACTGC-1 SeuratProject          20915             6318
#空转的话,就是SpatialFeaturePlot ,相较于单细胞转录组的FeaturePlot
SpatialFeaturePlot(spe, features = "nFeature_Spatial")

额,报错了,如下:

Error in FUN(left, right) : non-numeric argument to binary operator

推荐google搜索,解决方法如下(将spe@images$Lung@coordinates的信息改为):

for (i in colnames((spe@images$Lung@coordinates))) {
  spe@images$Lung@coordinates[[i]] <- as.integer(spe@images$Lung@coordinates[[i]])
}

再次尝试,问题解决,注意如果你的数据没有该报错可以不用as.integer 这部分!!!

SpatialFeaturePlot(spe, features = "nFeature_Spatial")
p0 <- SpatialDimPlot(spe,alpha = 0)
p1 <- SpatialDimPlot(spe,alpha = 1)
p0 + p1

image.png

非常建议通过 str(spe)查看一下空转的数据结构 ,有助于后续更好的分析!

2,10X的三个文件 + spatial文件夹

默认的Load10X_Spatial函数不能像单细胞转录组一样直接读取3个文件,通过View(Load10X_Spatial) 发现该函数写的比较死。其实我们只需要提取前面读取单细胞数据部分,使用Read10X即可


spe2 = Read10X("./outs/filtered_feature_bc_matrix/")
image2 <- Read10X_Image(image.dir = file.path("./outs/", 
                                              "spatial"), filter.matrix = TRUE)
spe2 <- CreateSeuratObject(counts = spe2, assay = "Spatial")
image2 <- image2[Cells(x = spe2)]
DefaultAssay(spe2 = image2) <- "Spatial"
spe2[["slice1"]] <- image2
#没有报错,无需转化
for (i in colnames((spe2@images$slice1@coordinates))) {
  spe2@images$slice1@coordinates[[i]] <- as.integer(spe2@images$slice1@coordinates[[i]])
}
SpatialFeaturePlot(spe2, features = "nFeature_Spatial")

Load10X_Spatial函数中的spe2[["slice1"]] 部分,也就是为什么在默认情况下image的slide 名字是slice1 ,所以如果多样本分析的时候记得改名字

3,读取rds文件

部分文献会提供rds文件,但是因为版本不一致常会遇到问题




spe3 <- readRDS("./STData/P1N_Spatial.rds")
SpatialFeaturePlot(spe3 , features = "MS4A1")

google检索解决方案:https://github.com/satijalab/seurat/issues/6312 ,可以UpdateSeuratObject 一下



spe3<- UpdateSeuratObject(spe3)
SpatialFeaturePlot(spe3, features = "MS4A1")

4,只提供 tissue_hires_image.png

spatial文件夹中会有tissue_hires_image.png和 tissue_lowres_image.png 两种tiff图,而Read10X_Image 函数中默认读取 tissue_lowres_image.png图,

因此如果(1)文献只提供了tissue_hires_image.png 或者(2)想绘制高清一些的图,可以使用如下的方法

4.1 尚无Seurat object

文献只提供hires图时候直接使用Load10X_Spatial会报错,可以先指定使用tissue_hires_image图片 ,如下

img = Read10X_Image("./outs/spatial/", 
                    image.name = "tissue_hires_image.png")
spe4 = Load10X_Spatial("./outs/", image = img)
spe4@images$slice1@scale.factors$lowres = spe4@images$slice1@scale.factors$hires
#没有报错,无需转化
for (i in colnames((spe4@images$slice1@coordinates))) {
  spe4@images$slice1@coordinates[[i]] <- as.integer(spe4@images$slice1@coordinates[[i]])
}
#对比一下两张图的差异
p1 <- SpatialFeaturePlot(spe2, features = "nFeature_Spatial")
p2 <- SpatialFeaturePlot(spe4, features = "nFeature_Spatial")
p1 + p2

image.png

其他的图可能会差异更明显。

4.2  已有Seurat object

如果前面已经创建了seurat object ,不想从头分析的话 也可以直接替换的,如下先重新定义一个spe5


spe5 <- spe2
lowres = spe5@images$slice1
hires = Read10X_Image("./outs/spatial/", 
                      image.name = "tissue_hires_image.png")
hires@scale.factors$lowres = hires@scale.factors$hires
#Set the assay and key to the same as the lowres (I'm not sure if this is necessary or not).
hires@assay = lowres@assay
hires@key = lowres@key
#Finally, add the image into the object and everything should be good!
spe5@images$slice1 = hires
#没有报错,无需转化
for (i in colnames((spe5@images$slice1@coordinates))) {
  spe5@images$slice1@coordinates[[i]] <- as.integer(spe5@images$slice1@coordinates[[i]])
}
p1 <- SpatialFeaturePlot(spe2, features = "nFeature_Spatial")
p2 <- SpatialFeaturePlot(spe5, features = "nFeature_Spatial")
p1 + p2

解答来源于https://github.com/satijalab/seurat/discussions/4833

5,只提供位置信息

没有提供spatial文件夹,只提供了空转spot的位置信息,可以按照单细胞的方式读取,但是无法绘制如上的空转切片为背景的图。

以上是总结的一些空转数据的读入方式,请查收!

相关文章
|
9月前
|
运维 Shell Python
【运维知识高级篇】超详细的Shell编程讲解2(变量切片+统计变量长度+字串删除+字串替换+七种方法进行数值运算+整数比较+多整数比较+文件判断+字符串比对+正则比对+配合三剑客的高阶用法)(一)
【运维知识高级篇】超详细的Shell编程讲解2(变量切片+统计变量长度+字串删除+字串替换+七种方法进行数值运算+整数比较+多整数比较+文件判断+字符串比对+正则比对+配合三剑客的高阶用法)
93 0
|
8月前
|
机器学习/深度学习 存储 C语言
二进制优化的快读模板,以及常用的读入形式。
二进制优化的快读模板,以及常用的读入形式。
42 0
|
9月前
|
算法 Linux Python
SGAT丨hapmap 格式hmp.txt文件转换,基因型和表型文件样品关联筛选提取的快速方法
SGAT丨hapmap 格式hmp.txt文件转换,基因型和表型文件样品关联筛选提取的快速方法
|
9月前
|
运维 Shell Perl
【运维知识高级篇】超详细的Shell编程讲解2(变量切片+统计变量长度+字串删除+字串替换+七种方法进行数值运算+整数比较+多整数比较+文件判断+字符串比对+正则比对+配合三剑客的高阶用法)(二)
【运维知识高级篇】超详细的Shell编程讲解2(变量切片+统计变量长度+字串删除+字串替换+七种方法进行数值运算+整数比较+多整数比较+文件判断+字符串比对+正则比对+配合三剑客的高阶用法)(二)
94 0
|
11月前
|
数据采集 移动开发 数据可视化
空间转录组|Load10X_Spatial函数修改适配多形式数据 + 空转标准流程
空间转录组|Load10X_Spatial函数修改适配多形式数据 + 空转标准流程
376 0
输入一个整形数(最多可以到亿位),然后按汉语的习惯,将其读出来并输出。如1052,读作:一千零五十二。 程序运行示例: 1052 一千零五十二
输入一个整形数(最多可以到亿位),然后按汉语的习惯,将其读出来并输出。如1052,读作:一千零五十二。 程序运行示例: 1052 一千零五十二
140 0
|
算法
改错题:用户从键盘任意输入一个数字表示月份值n,程序显示该月份对应的英文表示,若n不在1~12之间,则输出“Illegal month”。 注意: (1)请将修改正确后的完整源程序拷贝粘贴到答题区内。
改错题:用户从键盘任意输入一个数字表示月份值n,程序显示该月份对应的英文表示,若n不在1~12之间,则输出“Illegal month”。 注意: (1)请将修改正确后的完整源程序拷贝粘贴到答题区内。
259 0
读取文件结束的判定的概念,使用方法和文件缓冲区的位置
读取文件结束的判定的概念,使用方法和文件缓冲区的位置
105 0
|
SQL 数据采集 NoSQL
数据预处理-航线类型操作类型-读取规则到程序|学习笔记
快速学习数据预处理-航线类型操作类型-读取规则到程序
274 0
数据预处理-航线类型操作类型-读取规则到程序|学习笔记
键盘录入一组学生信息,保存到文件中,保存形式为每个学生名字占一行
键盘录入一组学生信息,保存到文件中,保存形式为每个学生名字占一行