劉小澤寫(xiě)于19.10.8
前幾天單細(xì)胞天地推送了一篇整合scRNA數(shù)據(jù)的文章：使用seurat3的merge功能整合8個(gè)10X單細(xì)胞轉(zhuǎn)錄組樣本
這次根據(jù)推送佑惠，再結(jié)合自己的理解寫(xiě)一寫(xiě)
前言
單細(xì)胞數(shù)據(jù)未來(lái)會(huì)朝著多樣本發(fā)展电媳，因此數(shù)據(jù)整合是一項(xiàng)必備技能。cellranger中自帶了aggr的整合功能但绕，而這篇文章（Differentiation dynamics of mammary epithelial cells revealed by single-cell RNA-sequencing）的作者也正是這么做得到的組合后的表達(dá)矩陣，然后用Read10X讀入
關(guān)于文章
這是發(fā)表在2017年10月的NC文章。
作者的論點(diǎn)是：乳腺上皮細(xì)胞對(duì)研究乳腺癌的發(fā)展很重要症副，但目前只有很少的marker可以追蹤這群細(xì)胞，因此有必要探索乳腺發(fā)育不同階段的乳腺上皮細(xì)胞變化政基。
實(shí)驗(yàn)涉及了四個(gè)時(shí)期：8 weeks virgin =》nulliparous (NP) 未懷孕時(shí)期贞铣、14.5d gestation (G) 妊娠期第14.5天、6d lactation (L) 哺乳期第6天沮明、11d post involution (PI) 完全退化第11天辕坝。其中每個(gè)時(shí)期都采集兩只老鼠的組織細(xì)胞，所以一共8個(gè)樣本荐健，然后使用10X建庫(kù)酱畅，那么最后的測(cè)序文件就是8*3 = 24個(gè)
如果要對(duì)這24個(gè)文件分別去整合琳袄，使用seurat的merge函數(shù)即可，不過(guò)問(wèn)題的關(guān)鍵是：如何用代碼將這些樣本區(qū)分開(kāi)纺酸，然后分別構(gòu)建對(duì)象窖逗，最后merge這些對(duì)象
開(kāi)始操作
第一步：準(zhǔn)備原始測(cè)序數(shù)據(jù)
我們下載第一個(gè)：GSE106273_RAW.tar（183.5 Mb） https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE106273&format=file
感覺(jué)手機(jī)下載速度就是比電腦快，這個(gè)文件在手機(jī)上下載3分鐘餐蔬，電腦預(yù)計(jì)時(shí)間1小時(shí)
下載后解壓碎紊，整個(gè)過(guò)程直接在Rstudio中的Terminal直接完成
第二步：整理數(shù)據(jù)
思路：根據(jù)中間的分組信息（NP、G）將包含的文件放到相應(yīng)的文件夾中
方法一：shell腳本
# 將同一組數(shù)據(jù)放在同一目錄下
ls GSM* | awk -F '_' '{print $2"_"$3}'| uniq | while read i;do mkdir $i;mv *$i*gz $i;done
# 各自重命名
find -name "*barcodes.tsv.gz" | while read i;do mv $i $(dirname $i)/barcodes.tsv.gz;done
find -name "*genes.tsv.gz" | while read i;do mv $i $(dirname $i)/genes.tsv.gz;done
find -name "*matrix.mtx.gz" | while read i;do mv $i $(dirname $i)/matrix.mtx.gz;done
方法二：R腳本
# 列出當(dāng)前目錄下所有開(kāi)頭是GSM的文件
fs=list.files('./','^GSM')
# 然后獲取四個(gè)樣本信息
library(stringr)
samples=str_split(fs,'_',simplify = T)[,1]
# 設(shè)置一個(gè)循環(huán)樊诺，對(duì)每個(gè)樣本信息做同樣的事：
#（1）找到包含這個(gè)樣本的文件(用grepl)
# (2)設(shè)置對(duì)應(yīng)的目錄名（str_split+paste）然后創(chuàng)建目錄（用dir.create）
# (3)將文件放到對(duì)應(yīng)目錄(采用的是file.rename)并重命名文件
lapply(unique(samples),function(x){
y=fs[grepl(x,fs)]
folder=paste(str_split(y[1],'_',simplify = T)[,2:3],
collapse = '')
dir.create(folder,recursive = T)
file.rename(y[1],file.path(folder,"barcodes.tsv.gz"))
file.rename(y[2],file.path(folder,"genes.tsv.gz"))
file.rename(y[3],file.path(folder,"matrix.mtx.gz"))
})
需要注意的是仗考，Read10X函數(shù)需要讀取解壓后的文件，于是還要對(duì)所有的數(shù)據(jù)文件進(jìn)行解壓
find ./ -name "*gz" |xargs gunzip
常見(jiàn)錯(cuò)誤：
說(shuō)找不到Barcode文件词爬，但明明存在Barcode：
Error in Read10X() : Barcode file missing那很有可能是因?yàn)槿齻€(gè)10X數(shù)據(jù)的命名出了問(wèn)題痴鳄，一定要命名成"barcodes.tsv" "genes.tsv""matrix.mtx"
【補(bǔ)充：cellranger的V2版本得到的結(jié)果分別是：barcodes.tsv、genes.tsv缸夹、matrix.mtx痪寻；V3版本得到的結(jié)果分別是：matrix.mtx.gz、features.tsv.gz虽惭、barcodes.tsv.gz】說(shuō)找不到基因文件橡类，那么就要看看測(cè)序數(shù)據(jù)是不是解壓后的
Error in Read10X() : Gene name or features file missing
第三步：批量讀取成10X對(duì)象
Read10X() + CreateSeuratObject()
# 因?yàn)镽ead10X函數(shù)需要對(duì)目錄進(jìn)行操作，所以先把目錄名提取出來(lái)
folders=list.files('./',pattern='[12]$')
> folders
[1] "G_1" "G_2" "L_1" "L_2"
[5] "NP_1" "NP_2" "PI_1" "PI_2"
# 然后使用lapply進(jìn)行循環(huán)(看下lapply的幫助文檔就知道芽唇，它是對(duì)列表或向量進(jìn)行循環(huán)顾画，而apply是對(duì)數(shù)據(jù)框或矩陣操作)
library(Seurat)
sceList = lapply(folders,function(folder){
CreateSeuratObject(counts = Read10X(folder),
project = folder )
})
# 此時(shí)的sceList僅僅是一個(gè)堆砌了8個(gè)10X對(duì)象的集合，下一步就要真正合并起來(lái)
> sceList
[[1]]
An object of class Seurat
27998 features across 2915 samples within 1 assay
Active assay: RNA (27998 features)
[[2]]
An object of class Seurat
27998 features across 3106 samples within 1 assay
Active assay: RNA (27998 features)
第四步：組合
sce.big <- merge(sceList[[1]],
y = c(sceList[[2]],sceList[[3]],sceList[[4]],
sceList[[5]],sceList[[6]],
sceList[[7]],sceList[[8]]),
add.cell.ids = folders,
project = "mouse8")
> table(sce.big$orig.ident)
G_1 G_2 L_1 L_2 NP_1 NP_2 PI_1 PI_2
2915 3106 5906 3697 2249 2127 1500 4306
save(sce.big,file = 'sce.big.merge.mouse8.Rdata') # 保存的數(shù)據(jù)是1.4G
補(bǔ)充
官網(wǎng)的merge教程在：https://satijalab.org/seurat/v3.1/merge_vignette.html
描述了三種情況
第一種： merge兩個(gè)seurat對(duì)象（原始數(shù)據(jù)）
需要注意的是匆笤，組合數(shù)據(jù)時(shí)需要注明每個(gè)數(shù)據(jù)的名稱研侣，使用add.cell.ids參數(shù)指定
pbmc.combined <- merge(pbmc4k, y = pbmc8k, add.cell.ids = c("4K", "8K"), project = "PBMC12K")
pbmc.combined
## An object of class Seurat
## 33694 features across 12721 samples within 1 assay
## Active assay: RNA (33694 features)
# 之后的組合數(shù)據(jù)就會(huì)出現(xiàn)列名的標(biāo)識(shí)
head(colnames(pbmc.combined))
## [1] "4K_AAACCTGAGAAGGCCT" "4K_AAACCTGAGACAGACC" "4K_AAACCTGAGATAGTCA"
## [4] "4K_AAACCTGAGCGCCTCA" "4K_AAACCTGAGGCATGGT" "4K_AAACCTGCAAGGTTCT"
table(pbmc.combined$orig.ident)
##
## PBMC4K PBMC8K
## 4340 8381
第二種：merge兩個(gè)以上（原始數(shù)據(jù)）
將參數(shù)y 設(shè)成一個(gè)向量，就可以指定其他的數(shù)據(jù)
pbmc.big <- merge(pbmc3k,
y = c(pbmc4k, pbmc8k),
add.cell.ids = c("3K", "4K", "8K"),
project = "PBMC15K")
unique(sapply(X = strsplit(colnames(pbmc.big), split = "_"), FUN = "[", 1))
## [1] "3K" "4K" "8K"
table(pbmc.big$orig.ident)
##
## pbmc3k PBMC4K PBMC8K
## 2638 4340 8381
第三種：merge歸一化炮捧、標(biāo)準(zhǔn)化數(shù)據(jù)
默認(rèn)情況庶诡，只會(huì)組合原始數(shù)據(jù)，但如果有的數(shù)據(jù)時(shí)標(biāo)準(zhǔn)化之后的呢咆课？
其實(shí)可以通過(guò)一個(gè)參數(shù)merge.data = TRUE指定
pbmc4k <- NormalizeData(pbmc4k)
pbmc8k <- NormalizeData(pbmc8k)
pbmc.normalized <- merge(pbmc4k,
y = pbmc8k,
add.cell.ids = c("4K", "8K"),
project = "PBMC12K",
merge.data = TRUE)
看看第一種組合raw data和第三種組合normalized data對(duì)比：
#################
# raw data
#################
GetAssayData(pbmc.combined)[1:10, 1:15]
## 10 x 15 sparse Matrix of class "dgCMatrix"
##
## RP11-34P13.3 . . . . . . . . . . . . . . .
## FAM138A . . . . . . . . . . . . . . .
## OR4F5 . . . . . . . . . . . . . . .
## RP11-34P13.7 . . . . . . . . . . . . . . .
## RP11-34P13.8 . . . . . . . . . . . . . . .
## RP11-34P13.14 . . . . . . . . . . . . . . .
## RP11-34P13.9 . . . . . . . . . . . . . . .
## FO538757.3 . . . . . . . . . . . . . . .
## FO538757.2 . . . . . . . . . 1 . . . . .
## AP006222.2 . . . . . . . . . . . 1 . . .
#################
# normalized data
#################
GetAssayData(pbmc.normalized)[1:10, 1:15]
## 10 x 15 sparse Matrix of class "dgCMatrix"
##
## RP11-34P13.3 . . . . . . . . . . . . . . .
## FAM138A . . . . . . . . . . . . . . .
## OR4F5 . . . . . . . . . . . . . . .
## RP11-34P13.7 . . . . . . . . . . . . . . .
## RP11-34P13.8 . . . . . . . . . . . . . . .
## RP11-34P13.14 . . . . . . . . . . . . . . .
## RP11-34P13.9 . . . . . . . . . . . . . . .
## FO538757.3 . . . . . . . . . . . . . . .
## FO538757.2 . . . . . . . . . 0.7721503 . . . . .
## AP006222.2 . . . . . . . . . . . 1.087928 . . .
上面的組合多個(gè)數(shù)據(jù)就結(jié)束了末誓，接下來(lái)是檢查組合后的分群結(jié)果
首先檢查原樣本分群結(jié)果
# 歸一化+標(biāo)準(zhǔn)化（移除了不想要的差異來(lái)源nCount_RNA）
sce.big <- NormalizeData(sce.big)
sce.big <- ScaleData(sce.big, vars.to.regress = c('nCount_RNA'),
model.use = 'linear', use.umi = FALSE)
# 默認(rèn)選2000個(gè)HVGs
sce.big <- FindVariableFeatures(object = sce.big,
mean.function = ExpMean,
dispersion.function = LogVMR,
mean.cutoff = c(0.0125,4),
dispersion.cutoff = c(0.5,Inf))
# 降維（PCA+tSNE）
sce.big <- RunPCA(object = sce.big, pc.genes = VariableFeatures(sce.big))
sce.big <- RunTSNE(object = sce.big, dims.use = 1:10)
DimPlot(object = sce.big, reduction = "tsne")
# 當(dāng)然也有ICA的選擇
# sce.big <- RunICA(sce.big )
然后鑒定亞群，看看它們的分群結(jié)果
ElbowPlot(sce.big)
# 官方建議书蚪，下游分析時(shí)可以多用幾個(gè)PCs試試
sce.big <- FindNeighbors(sce.big, dims = 1:20)
# 保持和原文一樣的15個(gè)亞群
sce.big <- FindClusters(sce.big, resolution = 0.23)
head(Idents(sce.big), 5)
# 新的亞群結(jié)果
DimPlot(object = sce.big, reduction = "tsne",
group.by = 'RNA_snn_res.0.23')
# 原樣本分群結(jié)果
DimPlot(object = sce.big, reduction = "tsne",
group.by = 'orig.ident')
table(sce.big$orig.ident,sce.big@meta.data$RNA_snn_res.0.23)
看到喇澡，NP有三個(gè)主群、G有三個(gè)主群殊校、L有三個(gè)主群晴玖、PI有兩個(gè)主群
對(duì)比原文的數(shù)據(jù)
它得到了3個(gè)NP、3個(gè)G、2個(gè)L呕屎、3個(gè)PI宪萄，其余的分給了Basal 4群
歡迎關(guān)注我們的公眾號(hào)～_～
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩，打造生信星球榨惰，想讓它成為一個(gè)不拽術(shù)語(yǔ)拜英、通俗易懂的生信知識(shí)平臺(tái)。需要幫助或提出意見(jiàn)請(qǐng)后臺(tái)留言或發(fā)送郵件到jieandze1314@gmail.com
