什么？一个单细胞样本变成俩了？

1.背景知识

众所周知，单细胞转录组数据的输入文件是多种多样的，详见:不同文件格式单细胞数据读取流程 https://mp.weixin.qq.com/s/W7szy-Kg6G1N1ENHNRjGiw

除了标准三个文件的格式，就还有很多其他的！

今天来看一个小坑坑。

这个数据是GSM4453576，给的是txt.gz格式，gz是压缩格式，不用解压可以正常读取。

所以正常的读取方式应该是这样：

a = data.table::fread("GSM4453576_P1_exp.txt.gz",data.table = F)
a = tibble::column_to_rownames(a,"V1")
a[1:4,1:4]

然后把这个表达矩阵传递给counts参数即可：

library(Seurat)
seu.obj <- CreateSeuratObject(counts = a, 
                              min.cells = 3, 
                              min.features = 200)
dim(seu.obj)

## [1] 20769  5317

2.发现bug并找到源头

但是呢，后面做质控小提琴图时就会发现，好好的一个样本居然被画在两个小提琴里了？按理来说单样本数据应该是只有一个才对。

VlnPlot(seu.obj,"nCount_RNA")

[图片上传失败...(image-3d653b-1721981237535)]

甚至有的数据是这样（来自学生）

image.jpeg

不仔细看都看不出来横坐标有3个，其中第三个只有一个点(可能是读取的代码有啥毛病吧，谁做错了参考我的作为答案就行，手动狗头)。

image.jpeg

还给了个warning，让设置drop=F，我起初以为这个参数就是解决问题的，实际上它只能解决点的数量太少的”分组“被删掉的问题，例如上图中的”V5381”。并不能解决org.ident的问题。

VlnPlot不指定颜色的时候，就是按照orig.ident分配颜色的，所以问题出在orig.ident上面：

table(seu.obj$orig.ident)

## 
##    1    2 
## 1596 3721

果然，orig.ident分了1和2。正常的单样本数据应该是同一个值，默认是SeuratProject ，如果CreateSeuratObject时设置了project参数的话，那么orig.ident的值就会是设置的值，没啥用，设不设置都一样用。

name为啥会这样呢，是因为CreateSeuratObject默认的样本识别机制：

帮助文档里面有两个参数：

names.field:For the initial identity class for each cell, choose this field from the cell’s name. E.g. If your cells are named as BARCODE_CLUSTER_CELLTYPE in the input matrix, set names.field to 3 to set the initial identities to CELLTYPE.

names.delim:or the initial identity class for each cell, choose this delimiter from the cell’s column name. E.g. If your cells are named as BARCODE-CLUSTER-CELLTYPE, set this to “-” to separate the cell name into its component parts for picking the relevant field.

其中names.field是指定列名里面的第几部分是初始分类(初始分类就是样本) names.delim是分隔符，上面说索的”第几部分“指的是按照某分隔符分割后的”第几部分“，默认值是_。

所以按照帮助文档里的例子就很好理解了：假如列名是BARCODE_CLUSTER_CELLTYPE的格式，默认值就是把按照”“为分隔符(names.delim的默认值为”)拆分出来的第1部分(names.field默认值为1)作为orig.ident！

这样也就理解了它是怎么给我们把一个样本拆成两个的：

library(stringr)
colnames(a) %>%
  str_split_i("_",1) %>%
  table()

## .
##    1    2 
## 1596 3721

和上面table(seu.obj$orig.ident)的结果是呼应的。

所以问题就是列名第一部分不统一。

3.解决办法

有很多种办法可以解决：

例如names.field = 2，忽略第一部分，按照第二部分，但第二部分是真正的barcode，每个细胞都不一样，所以就和正常的、没有前缀的数据会同等处理，都设为一个值。

seu.obj <- CreateSeuratObject(counts = a, 
                              min.cells = 3, 
                              min.features = 200,
                              names.field = 2)
table(seu.obj$orig.ident)

## 
## SeuratProject 
##          5317

再例如names.delim=“,”，逗号可以换成任何一个在列名里不存在的符号，这样就会把全部列名识别为第一个部分。

seu.obj <- CreateSeuratObject(counts = a, 
                              min.cells = 3, 
                              min.features = 200,
                              names.delim = ",")
table(seu.obj$orig.ident)

## 
## SeuratProject 
##          5317

再例如改列名，删掉前缀，不过这种方法有点风险，某些抽风数据会有重复的barcode。

colnames(a) = str_split_i(colnames(a) ,"_",2) 
seu.obj <- CreateSeuratObject(counts = a, 
                              min.cells = 3, 
                              min.features = 200)
table(seu.obj$orig.ident)

## 
## SeuratProject 
##          5317

其实还有别的方法，比如RenameIdents(),就不一一列举了。

4.总结

需要了解单细胞对象，多多积累经验，寻找bug源头并杜绝它。

我发现后时候大家是被动的去解决问题，而不是寻找源头来杜绝，比如某位学员的思维是把1和2合并然后去批次！

事实是这个数据是一个样本，被分成两个仅仅是函数的默认参数与数据实际情况不符，这不是真的两个样本，既然被错误的拆分，我们就应该避免这个错误，而不是在错误的基础上继续往前走，都没有批次效应，何来去批次之说呢？

越发感觉到培养解决问题的思维比解决问题本身更重要。