虽然现在(2023年12月30日)空间单细胞技术已经是百花齐放了,主要是基于测序的10x Genomics Visium平台和Slide-seq技术,另外基于成像的技术也有一点点市场,比如:seqFISH+、MERFISH、 NanoString GeoMx Digital Spatial Profiler (DSP),还有其它小众产品就不值得一提啦。
是否需要 spaceranger count 的定量过程
而且,通俗来说,单细胞转录组就是10x,空间单细胞也是10x,所以我们只能说是一起讨论10x技术的数据分析策略。
那么就不得不预先理解一下10x的空间单细胞文件格式,因为都是10x技术, 所以它也是正常走10x官网软件即可,比如前面我们介绍的单细胞转录组的cellranger的定量流程即可,代码我已经是多次分享了。参考:
- 10X单细胞转录组原始测序数据的Cell Ranger流程(仅需800元)
- 10X的单细胞转录组原始数据也可以在EBI下载
- 一个10x单细胞转录组项目从fastq到细胞亚群
- 一文打通单细胞上游:从软件部署到上游分析
- PRJNA713302这个10x单细胞fastq实战
- 一次曲折且昂贵的单细胞公共数据获取与上游处理
- 只能下载bam文件的10x单细胞转录组项目数据处理
- 不知道10x单细胞转录组样品和fastq文件的对应关系
- 10X单细胞转录组测序数据的 SRA转fastq踩坑那些事
- 10x的单细胞转录组fastq文件的R1和R2不能弄混哦
值得注意的是10x的空间单细胞使用的是Space Ranger,软件下载以及数据库文件压缩包下载:
虽然说Space Ranger跟cellranger软件名字不一样,但是很类似,感兴趣的可以试试看10x官网提供的案例数据:
基本上也是常规的fq文件命名格式即可,因为不同文件命名有不同含义,所以一定是不能弄混哦 :
- 28bp read 1 (16bp Visium spatial barcode, 12bp UMI),
- 50bp read2 (transcript)
- 10bp i7 sample barcode
- 10bp i5 sample barcode.
对上面的案例数据(Visium_FFPE_Human_Breast_Cancer_fastqs.tar),走 spaceranger 代码
raw=/home/x10/visium/visium_test/FFPE/rawdata
fq_dir=$raw/Visium_FFPE_Human_Breast_Cancer_fastqs
ref=/home/x10/pipeline/refdata-gex-GRCh38-2020-A/
spaceranger count --id=Visium_FFPE_Human_Breast_Cancer \
--transcriptome=$ref \
--fastqs=$fq_dir \
--sample=Visium_FFPE_Human_Breast_Cancer \
--image=$raw/Visium_FFPE_Human_Breast_Cancer_image.tif \
--slide=V11J26-008 \
--area=B1 \
--probe-set=$raw/Visium_FFPE_Human_Breast_Cancer_probe_set.csv \
--localcores=20 \
--localmem=80
20核用时半个小时即可,如果有服务器的话,这个流程非常简单!
不过绝大部分情况下目前的空间单细胞转录组文章都是提供了这个 spaceranger count 后的结果文件,我们要做的就是理解10x的空间单细胞文件格式详解即可,大概率是不需要跑这个 spaceranger count 本身。它输出文件是非常多的,详见官网介绍:
- https://software.10xgenomics.com/spatial-gene-expression/software/pipelines/latest/output/overview
- https://www.10xgenomics.com/support/software/space-ranger/analysis/outputs/space-ranger-count-outputs
但是我们大概是不需要全部都讲解,最重要的是首先是表达量矩阵和图片信息:
filtered_feature_bc_matrix/ # 文件夹
filtered_feature_bc_matrix.h5 #表达量矩阵 in HDF5 format.
spatial/ # 空间图片信息文件夹
web_summary.html # 网页报表
其中web_summary.html 这个 网页报表值得单独拿出来讲解因为它能告诉你你的空间单细胞数据的产出的质量,不过我们默认都是ok的,所以就不用看这个网页版报表文件啦。
然后我们来分开讲解表达量矩阵和空间图片信息。
首先是表达量矩阵
目前在单细胞转录组学中,表达量矩阵可以以不同的格式存储,其中 Market Exchange Format (MEX) 和 Hierarchical Data Format (HDF5) 是两种不同的格式。
其中Market Exchange Format (MEX)格式的表达量矩阵有3个文件,如下所示:
$ cd /home/jdoe/runs/sample345/outs
$ tree filtered_feature_bc_matrix
filtered_feature_bc_matrix
├── barcodes.tsv.gz
├── features.tsv.gz
└── matrix.mtx.gz
0 directories, 3 files
而且这3个文件必须是有固定的格式以及固定好的文件名字,不能修改!!!因为我们读取它的时候只需要文件夹的名字,文件夹里面的3个文件是一定要固定的!分别存储 列名(细胞barcode),行名(基因名字),表达量矩阵(稀疏矩阵格式)。
在R或者Python编程语言里面的,这3个文件都是可以分开独立读取的。使用 Python 的 scipy
库或 R 的 Matrix
和 data.table
等库来分别读取这三个单细胞转录组文件。
在Python里面是 csv
, os
, gzip
, and scipy.io
modules ,下面是 Python 代码示例:
import scipy.io
import pandas as pd
# 读取 matrix.mtx.gz 文件
matrix = scipy.io.mmread('matrix.mtx.gz')
# 读取 barcodes.tsv.gz 文件
barcodes_df = pd.read_csv('barcodes.tsv.gz', header=None, names=['barcodes'])
# 读取 features.tsv.gz 文件
features_df = pd.read_csv('features.tsv.gz', header=None, names=['features'])
# 打印示例数据
print("Matrix Shape:", matrix.shape)
print("Barcodes DataFrame:")
print(barcodes_df.head())
print("Features DataFrame:")
print(features_df.head())
下面是R 代码示例:
library(Matrix)
library(data.table)
# 读取 matrix.mtx.gz 文件
matrix <- readMM("matrix.mtx.gz")
# 读取 barcodes.tsv.gz 文件
barcodes_df <- fread("barcodes.tsv.gz", header = FALSE, col.names = c("barcodes"))
# 读取 features.tsv.gz 文件
features_df <- fread("features.tsv.gz", header = FALSE, col.names = c("features"))
# 打印示例数据
cat("Matrix Dimensions:", dim(matrix), "\n")
cat("Barcodes DataFrame:\n")
head(barcodes_df)
cat("Features DataFrame:\n")
head(features_df)
这些代码将帮助你分别读取单细胞转录组的矩阵、条形码和特征文件。请确保你的 Python 环境中已经安装了 scipy
和 pandas
库,而 R 环境中已经安装了 Matrix
和 data.table
包。如果你熟悉这两个编程语言,安装模块和包对你来说应该是手到擒来的事情。
而Hierarchical Data Format (HDF5) 格式的单细胞表达量矩阵就更简单了,简简单单的h5后缀的文件就存储上面的3个文件的全部信息,所以大概率是不需要大家分开读取上面的3个文件的,因为可以一次性读取。
然后是空间图像信息文件
同样的,这个spatial/ 空间图片信息文件夹里面的文件也并不是全部需要,首先是空间单细胞切片的图片本身:
tissue_hires_image.png
,分辨率很高,所以图片文件就很大。tissue_lowres_image.png
: 文件很小,因为图片分辨率低,官网说它 is always 600 pixels.
然后是 scalefactors_json.json 它关联了tissue_hires_image.png
and tissue_lowres_image.png
,案例内容如下所示:
$ cd /home/jdoe/runs/sampele345/outs/spatial
$ cat scalefactors_json.json
{"tissue_hires_scalef": 0.17011142,
"tissue_lowres_scalef": 0.051033426",
"fiducial_diameter_fullres": 144.4773339857,
"spot_diameter_fullres": 89.43834961021503}
感兴趣计算公式的可以自己看官网:https://www.10xgenomics.com/support/software/space-ranger/analysis/outputs/spatial-outputs
最后是 tissue_positions.csv
需要看看,它记录了对应的空间切片里面的每个barcode所在的切片的空间坐标,理论上有了这个csv文件其实就可以相当于有了空间切片的图片的细胞定位作用啦:
$ cd /home/jdoe/runs/sample345/outs/spatial
$ column -s, -t < tissue_positions.csv | head -n 10
barcode in_tissue array_row array_col pxl_row_in_fullres pxl_col_in_fullres
GTCACTTCCTTCTAGA-1 0 0 0 1832 11971
CACGGTCTCCTTACGA-1 0 0 2 1832 11833
ATAGCTGCGGATAAGA-1 0 0 4 1832 11695
GTCAGTATGTCCGGCG-1 0 0 6 1832 11556
ATGTACCAGTTACTCG-1 0 0 8 1831 11418
ACGCTCAGTGCACCGT-1 0 0 10 1831 11280
TCACTAACGTATAGTT-1 0 0 12 1831 11142
CGGTTAGGCCTGGACG-1 0 0 14 1831 11004
GATATCACCAGCATGG-1 0 0 16 1831 10866
简单的看几个案例
比如这个 2020-GSE111672-胰腺导管细胞的,其中一个样品GSM3036911_PDAC-A的文件如下所示 :
% ls -lh *|cut -d" " -f 7-
661B 9 22 2018 GSM3036911_PDAC-A-ST1-Cy3.jpg.gz
280M 8 26 2020 GSM3036911_PDAC-A-ST1-HE.jpg.gz
717K 11 4 2020 GSM3036911_PDAC-A-ST1-filtered.txt.gz
然后看2021-GSE158328-肠道发育的,自己下载 GSM4797916_A1.tar.gz 然后解压可以看到它每个样品其实有两个文件夹 :
% ls -lh *|cut -d" " -f 7-
raw_feature_bc_matrix:
28K 2 20 2020 barcodes.tsv.gz
298K 2 20 2020 features.tsv.gz
18M 2 20 2020 matrix.mtx.gz
spatial:
1.6M 2 20 2020 aligned_fiducials.jpg
1.8M 2 20 2020 detected_tissue_image.jpg
163B 2 20 2020 scalefactors_json.json
4.6M 2 20 2020 tissue_hires_image.png
496K 2 20 2020 tissue_lowres_image.png
180K 2 20 2020 tissue_positions_list.csv
接着看2022-复旦高强课题组-结直肠癌肝转移空间单细胞数据:
# 空间单细胞的文件在文献pdf里面的可以找到地址:ST-colon2
% ls -lh *|cut -d" " -f 7-
77K 12 15 2021 barcodes.tsv
162K 12 15 2021 coordinates.tsv
15M 12 19 2021 filtered_feature_bc_matrix.h5
601K 12 15 2021 genes.tsv
127M 12 15 2021 matrix.mtx
166B 12 19 2021 scalefactors_json.json
16M 12 19 2021 tissue_hires_image.png
1.4M 12 19 2021 tissue_lowres_image.png
184K 12 19 2021 tissue_positions_list.csv
这个文献对每个样品提供了Market Exchange Format (MEX) 和 Hierarchical Data Format (HDF5) 两种格式的表达量矩阵,这个时候我们一般来说会选择那个 filtered_feature_bc_matrix.h5 读取后分析即可。
再看看2022-GSE213216-子宫内膜-2个空间单细胞,我们挑选其中一个样品 GSM6690475 展示给大家看看 :
% ls -lh *|cut -d" " -f 7-
8.0K 10 27 2022 GSM6690475_BEME_346_barcodes.tsv.gz
326K 10 27 2022 GSM6690475_BEME_346_features.tsv.gz
12M 10 27 2022 GSM6690475_BEME_346_matrix.mtx.gz
175B 10 28 2022 GSM6690475_BEME_346_scalefactors_json.json.gz
4.7M 10 28 2022 GSM6690475_BEME_346_tissue_hires_image.png.gz
65K 10 28 2022 GSM6690475_BEME_346_tissue_positions_list.csv.gz
再看看2023-GSE217843-小鼠胰腺癌模型空间单细胞2个样品, 可以看到其中的GSM6727528_Tumor_DAY_30_ST_A1给出来的文件也很多:
% ls -lh *|cut -d" " -f 7-
1.2M 11 9 2022 GSM6727528_Tumor_DAY_30_ST_A1_aligned_fiducials.jpg.gz
19K 11 9 2022 GSM6727528_Tumor_DAY_30_ST_A1_barcodes.tsv.gz
1.2M 11 9 2022 GSM6727528_Tumor_DAY_30_ST_A1_detected_tissue_image.jpg.gz
284K 11 9 2022 GSM6727528_Tumor_DAY_30_ST_A1_features.tsv.gz
37M 11 10 2022 GSM6727528_Tumor_DAY_30_ST_A1_matrix.mtx.gz
150B 11 10 2022 GSM6727528_Tumor_DAY_30_ST_A1_scalefactors_json.json.gz
3.5M 11 10 2022 GSM6727528_Tumor_DAY_30_ST_A1_tissue_hires_image.png.gz
323K 11 10 2022 GSM6727528_Tumor_DAY_30_ST_A1_tissue_lowres_image.png.gz
65K 11 10 2022 GSM6727528_Tumor_DAY_30_ST_A1_tissue_positions_list.csv.gz
1.1M 11 10 2022 GSM6727529_Tumor_DAY_30_ST_B1_aligned_fiducials.jpg.gz
值得注意的是这个GSE217843数据集里面作者给出来了aligned_fiducials.jpg
: This image has the dimensions of tissue_hires_image.png
. Fiducial spots found by the fiducial alignment algorithm are highlighted in red. This file is useful to verify that fiducial alignment was successful.
上面的这些案例,空间单细胞文献,作者给大家的数据文件格式都不一样,但是都是可以分析的!怕的就是作者忘记了给任何空间信息文件给你,仅仅是提供了表达量矩阵, 这个就基本上无解了只能说是去写邮件给作者索取了。但是也有一些情况下是作者非常贴心的直接把空间单细胞数据使用R或者Python编程语言读取并且整理好存储为了编程语言里面的对象文件,你直接load即可,那就是最方便的。比如文章 :《Delineating the dynamic evolution from preneoplasia to invasive lung adenocarcinoma by integrating single-cell RNA sequencing and spatial transcriptomics》数据在GSE189487,每个病人有如下所示的4个数据文件:
GSM5702473_TD1_barcodes.tsv.gz 22.8 Kb
GSM5702473_TD1_features.tsv.gz 325.6 Kb
GSM5702473_TD1_matrix.mtx.gz 52.0 Mb
GSM5702473_TD1_tissue_positions_list.csv.gz 57.4 Kb
不过,它提供了在Dropbox的rdata文件可以下载,每个样品都是一个独立rdata文件,直接在r语言里面load即可。