工具不好用你完全可以自己写一个
工具不好用你完全可以自己写一个
理论上处理10x技术的单细胞转录组数据应该是毫无悬念了，因为太大众了所以到处是教程。但有些时候会被作者不小心埋下的雷给坑一下，比如这个人类膀胱癌数据集GSE181294明明写的是 hromium Single Cell 3’ Library and Gel Bead Kit v2 (10X Genomics). 但是读取就会失败，让我们拆解一下这个过程。
数据集是：GSE181294
451K 5 12 2022 GSM6133917_S1.barcode.csv.gz
40M 5 12 2022 GSM6133917_S1.counts.mtx.gz
72K 5 12 2022 GSM6133917_S1.genes.csv.gz
462K 5 12 2022 GSM6133918_S2.barcode.csv.gz
46M 5 12 2022 GSM6133918_S2.counts.mtx.gz
72K 5 12 2022 GSM6133918_S2.genes.csv.gz
。。。。
让人工智能大模型进行文件解释：
10x Genomics平台进行单细胞RNA测序时，会生成几种特定格式的文件，每种文件包含不同类型的数据。以下是您提供的文件名及其含义的解释：
- **Barcode文件 (
.barcode.csv.gz)**：- 这个文件包含了每个样本的条形码信息，通常用于追踪样本的身份和索引。
- 文件名示例：
GSM6133917_S1.barcode.csv.gz - 这里的
GSM6133917可能是一个样本的唯一标识符，S1表示样本的编号或批次。
- **Counts文件 (
.counts.mtx.gz)**：- 这个文件包含了基因表达的计数数据，通常是一个矩阵格式（Matrix Market格式），其中行代表基因，列代表单个细胞，元素值表示每个细胞中特定基因的读数或计数。
- 文件名示例：
GSM6133917_S1.counts.mtx.gz - 这个文件是基因表达分析中的核心数据，用于后续的数据处理和分析。
- **Genes文件 (
.genes.csv.gz)**：- 这个文件包含了基因的信息，通常包括基因标识符、基因名称和其他相关的基因注释信息。
- 文件名示例：
GSM6133917_S1.genes.csv.gz - 这个文件用于将计数数据与具体的基因实体关联起来，是进行基因表达分析的重要参考。
这些文件是单细胞RNA测序数据分析的基础，通常需要使用专门的软件和工具进行处理和分析。例如，.mtx.gz文件可以使用如Scanpy或Seurat等生物信息学软件包进行读取和分析，而.barcode.csv.gz和.genes.csv.gz文件则提供了必要的上下文信息，以确保数据的正确解释和分析。
因为后缀是.genes.csv.gz文件，所以可以比较肯定的是它属于早期10x技术的文件了。
使用Seurat官方Read10X函数报错：
Seurat官方Read10X函数可以读取一个文件夹，里面只需要有符合文件名规则的3个文件即可。但是读取这个数据集就报错了，信息如下所示：
Error in fixupDN.if.valid(value, x@Dim) :
length of Dimnames[[2]] (36979) is not equal to Dim[2] (36978)
其实打开这个Read10X函数可以看到它也是分开读取3个文件：
data <- readMM(file = matrix.loc)
cell.barcodes <- as.data.frame(data.table::fread(barcode.loc,
header = FALSE))
cell.names <- cell.barcodes[, cell.column]
cell.names <- readLines(con = barcode.loc)
feature.names <- as.data.frame(data.table::fread(ifelse(test = pre_ver_3,
yes = gene.loc, no = features.loc), header = FALSE))
如果Read10X函数报错，就应该是进去看具体是3个文件读取过程的哪个环节出问题。
解析函数以及对应的3个文件：
我们随意打开一个样品的3个文件，可以看到是 26369个基因 在 36978个细胞的表达量矩阵。
> d='10x/output/GSM6133917_S1'
> mtx=readMM( file.path(d,'matrix.mtx.gz') )
> mtx[1:4,1:4]
4 x 4 sparse Matrix of class "dgTMatrix"
[1,] . . . .
[2,] . . . .
[3,] . . . .
[4,] . . 3 .
> dim(mtx)
[1] 26369 36978
> cl=fread( file.path(d,'barcodes.tsv.gz') ,
+ header = T,data.table = F )
> dim(cl)
[1] 36978 3
> head(cl)
barcodes xcoord ycoord
1 HP1_TGCGGTGTCCCTGT 1135.50 1261.10
2 HP1_ACTTGACGGTAATA 3070.50 1922.00
3 HP1_AACGCMCTCACGTT 1099.40 1581.85
4 HP1_CCCTTCAGGGCGTS 3132.15 2846.30
5 HP1_TCCAACGGCCGTCC 2208.50 636.17
6 HP1_TTTATTAACCTCAT 3647.40 3919.50
> rl=fread( file.path(d,'features.tsv.gz') ,
+ header = T,data.table = F )
> head(rl)
x
1 A1BG
2 A1BG-AS1
3 A1CF
4 A2M
5 A2M-AS1
6 A2ML1
> dim(rl)
[1] 26369 1
其实没必要修改作者给出来的3个文件，因为我们可以自己模仿Seurat官方Read10X函数写一个实现类似的功能的工具：
构建“粗糙”的函数
需要简单的理解3个文件，以及Seurat的输入格式要求，我写的函数如下所示：
# 参考：https://mp.weixin.qq.com/s/tw7lygmGDAbpzMTx57VvFw
my_read10x <- function(d){
library(data.table)
library(Matrix)
# https://hbctraining.github.io/scRNA-seq/lessons/readMM_loadData.html
#
#d='10x/output/GSM6133917_S1'
mtx=readMM( file.path(d,'matrix.mtx.gz') )
mtx[1:4,1:4]
dim(mtx)
cl=fread( file.path(d,'barcodes.tsv.gz') ,
header = T,data.table = F )
head(cl)
rl=fread( file.path(d,'features.tsv.gz') ,
header = T,data.table = F )
head(cl)
dim(cl)
head(rl)
dim(mtx)
# length(unique(rl$V1))
# kp=!duplicated(rl$V1)
# table(kp)
# dim(mtx)
# mtx=mtx[kp,]
# rl=rl[kp,]
rownames(mtx)=rl[,1]
colnames(mtx)=cl[,1] #paste0('c',cl$V2)
mtx
}
接下来就可以使用我的my_read10x去读取每个样品的3个文件啦。然后就可以降维聚类分群，值得注意的是这个这个人类膀胱癌数据集GSE181294其实有两部分单细胞项目，另外一个项目给出来了的是csv文件格式， 也是可以自己解析它然后凑成为Seurat的输入即可构建对象：
ls -lh *gz|cut -d" " -f 7-
5.4M 8 2 2021 GSM5494342_HP1.count.csv.gz
5.3M 8 2 2021 GSM5494343_HP2.count.csv.gz
7.6M 8 2 2021 GSM5494344_HP3.count.csv.gz
5.0M 8 2 2021 GSM5494345_HP4.count.csv.gz
3.9M 8 2 2021 GSM5494346_HP0.count.csv.gz
8.7M 8 2 2021 GSM5494347_SCG-PCA21-N-LG.count.csv.gz
4.9M 8 2 2021 GSM5494348_SCG-PCA21-T-LG.count.csv.gz
。。。。。
也是可以很容易的降维聚类分群：
降维聚类分群
腾讯云开发者
