1.背景知識(shí)
眾所周知踢步,單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的輸入文件是多種多樣的癣亚,詳見:不同文件格式單細(xì)胞數(shù)據(jù)讀取流程 https://mp.weixin.qq.com/s/W7szy-Kg6G1N1ENHNRjGiw
除了標(biāo)準(zhǔn)三個(gè)文件的格式,就還有很多其他的获印!
今天來看一個(gè)小坑坑述雾。
這個(gè)數(shù)據(jù)是GSM4453576,給的是txt.gz格式,gz是壓縮格式玻孟,不用解壓可以正常讀取唆缴。
所以正常的讀取方式應(yīng)該是這樣:
a = data.table::fread("GSM4453576_P1_exp.txt.gz",data.table = F)
a = tibble::column_to_rownames(a,"V1")
a[1:4,1:4]
然后把這個(gè)表達(dá)矩陣傳遞給counts參數(shù)即可:
library(Seurat)
seu.obj <- CreateSeuratObject(counts = a,
min.cells = 3,
min.features = 200)
dim(seu.obj)
## [1] 20769 5317
2.發(fā)現(xiàn)bug并找到源頭
但是呢,后面做質(zhì)控小提琴圖時(shí)就會(huì)發(fā)現(xiàn)黍翎,好好的一個(gè)樣本居然被畫在兩個(gè)小提琴里了面徽?按理來說單樣本數(shù)據(jù)應(yīng)該是只有一個(gè)才對(duì)。
VlnPlot(seu.obj,"nCount_RNA")
[圖片上傳失敗...(image-3d653b-1721981237535)]
甚至有的數(shù)據(jù)是這樣(來自學(xué)生)
不仔細(xì)看都看不出來橫坐標(biāo)有3個(gè)匣掸,其中第三個(gè)只有一個(gè)點(diǎn)(可能是讀取的代碼有啥毛病吧趟紊,誰做錯(cuò)了參考我的作為答案就行,手動(dòng)狗頭)碰酝。
還給了個(gè)warning霎匈,讓設(shè)置drop=F,我起初以為這個(gè)參數(shù)就是解決問題的送爸,實(shí)際上它只能解決點(diǎn)的數(shù)量太少的”分組“被刪掉的問題唧躲,例如上圖中的”V5381”。并不能解決org.ident的問題碱璃。
VlnPlot不指定顏色的時(shí)候弄痹,就是按照orig.ident分配顏色的,所以問題出在orig.ident上面:
table(seu.obj$orig.ident)
##
## 1 2
## 1596 3721
果然嵌器,orig.ident分了1和2肛真。正常的單樣本數(shù)據(jù)應(yīng)該是同一個(gè)值,默認(rèn)是SeuratProject 爽航,如果CreateSeuratObject時(shí)設(shè)置了project參數(shù)的話蚓让,那么orig.ident的值就會(huì)是設(shè)置的值,沒啥用讥珍,設(shè)不設(shè)置都一樣用历极。
name為啥會(huì)這樣呢,是因?yàn)镃reateSeuratObject默認(rèn)的樣本識(shí)別機(jī)制:
幫助文檔里面有兩個(gè)參數(shù):
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是分隔符衷佃,上面說索的”第幾部分“指的是按照某分隔符分割后的”第幾部分“趟卸,默認(rèn)值是_。
所以按照幫助文檔里的例子就很好理解了:假如列名是BARCODE_CLUSTER_CELLTYPE的格式氏义,默認(rèn)值就是把按照”“為分隔符(names.delim的默認(rèn)值為”)拆分出來的第1部分(names.field默認(rèn)值為1)作為orig.ident锄列!
這樣也就理解了它是怎么給我們把一個(gè)樣本拆成兩個(gè)的:
library(stringr)
colnames(a) %>%
str_split_i("_",1) %>%
table()
## .
## 1 2
## 1596 3721
和上面table(seu.obj$orig.ident)的結(jié)果是呼應(yīng)的。
所以問題就是列名第一部分不統(tǒng)一惯悠。
3.解決辦法
有很多種辦法可以解決:
例如names.field = 2邻邮,忽略第一部分,按照第二部分克婶,但第二部分是真正的barcode筒严,每個(gè)細(xì)胞都不一樣丹泉,所以就和正常的、沒有前綴的數(shù)據(jù)會(huì)同等處理鸭蛙,都設(shè)為一個(gè)值嘀掸。
seu.obj <- CreateSeuratObject(counts = a,
min.cells = 3,
min.features = 200,
names.field = 2)
table(seu.obj$orig.ident)
##
## SeuratProject
## 5317
再例如names.delim=“,”,逗號(hào)可以換成任何一個(gè)在列名里不存在的符號(hào)规惰,這樣就會(huì)把全部列名識(shí)別為第一個(gè)部分睬塌。
seu.obj <- CreateSeuratObject(counts = a,
min.cells = 3,
min.features = 200,
names.delim = ",")
table(seu.obj$orig.ident)
##
## SeuratProject
## 5317
再例如改列名,刪掉前綴歇万,不過這種方法有點(diǎn)風(fēng)險(xiǎn)揩晴,某些抽風(fēng)數(shù)據(jù)會(huì)有重復(fù)的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
其實(shí)還有別的方法贪磺,比如RenameIdents(),就不一一列舉了硫兰。
4.總結(jié)
需要了解單細(xì)胞對(duì)象,多多積累經(jīng)驗(yàn)寒锚,尋找bug源頭并杜絕它劫映。
我發(fā)現(xiàn)后時(shí)候大家是被動(dòng)的去解決問題,而不是尋找源頭來杜絕刹前,比如某位學(xué)員的思維是把1和2合并然后去批次泳赋!
事實(shí)是這個(gè)數(shù)據(jù)是一個(gè)樣本,被分成兩個(gè)僅僅是函數(shù)的默認(rèn)參數(shù)與數(shù)據(jù)實(shí)際情況不符喇喉,這不是真的兩個(gè)樣本祖今,既然被錯(cuò)誤的拆分,我們就應(yīng)該避免這個(gè)錯(cuò)誤拣技,而不是在錯(cuò)誤的基礎(chǔ)上繼續(xù)往前走千诬,都沒有批次效應(yīng),何來去批次之說呢膏斤?
越發(fā)感覺到培養(yǎng)解決問題的思維比解決問題本身更重要徐绑。