工具不好用你完全可以自己写一个

生信技能树 · 公众号 · · 2024-11-27 10:39

正文

理论上处理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 cell.barcodes                                                  header = FALSE))
cell.names cell.names feature.names 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
  。。。。。

也是可以很容易的降维聚类分群：

降维聚类分群

学徒作业

完成这个人类膀胱癌数据集GSE181294的两部分单细胞项目，分开做哦，然后对比文献 Dissecting the immune suppressive human prostate tumor microenvironment via integrated single-cell and spatial transcriptomic analyses. Nat Commun 2023 Feb ，PMID: 36750562 看看是否合理

这个时候如果你也想做单细胞转录组数据分析，最好是有自己的计算机资源哦，比如我们的 2024的共享服务器交个朋友福利价仍然是800 ，而且还需要有基本的生物信息学基础，也可以看看我们的生物信息学马拉松授课（买一得五），你的生物信息学入门课。

如果你已经熟悉了我们的课程，就联系我们报名吧~
（添加好友务必备注高校或者工作单位+姓名+马拉松，方便后续认识）

生信入门班：
学习以转录组数据为代表的组学数据分析，包括上游分析（从下机数据到表达矩阵）和下游分析（差异分析、富集分析等），无专业偏向性，顺带学习基因表达芯片。
R语言是为下游分析打基础，linux是为上游分析打基础。