使用merge整合多組單細(xì)胞數(shù)據(jù)

開始操作

第一步：準(zhǔn)備原始測序數(shù)據(jù)

我們下載第一個：GSE106273_RAW.tar（183.5 Mb）

下載后解壓，整個過程直接在Rstudio中的Terminal直接完成

image

第二步：整理數(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)前目錄下所有開頭是GSM的文件
fs=list.files('./','^GSM')
# 然后獲取四個樣本信息
library(stringr)
samples=str_split(fs,'_',simplify = T)[,1]
# 設(shè)置一個循環(huán)，對每個樣本信息做同樣的事：
#（1）找到包含這個樣本的文件(用grepl)
# (2)設(shè)置對應(yīng)的目錄名（str_split+paste）然后創(chuàng)建目錄（用dir.create）
# (3)將文件放到對應(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ù)需要讀取解壓后的文件，于是還要對所有的數(shù)據(jù)文件進(jìn)行解壓

find ./ -name "*gz" |xargs  gunzip

常見錯誤：

  • 說找不到Barcode文件，但明明存在Barcode：
```
Error in Read10X() : Barcode file missing

```





那很有可能是因為三個10X數(shù)據(jù)的命名出了問題，一定要命名成"barcodes.tsv" "genes.tsv""matrix.mtx"
【補充：**cellranger的V2版本**得到的結(jié)果分別是：barcodes.tsv、genes.tsv、matrix.mtx；**V3版本**得到的結(jié)果分別是：matrix.mtx.gz、features.tsv.gz、barcodes.tsv.gz】
  • 說找不到基因文件，那么就要看看測序數(shù)據(jù)是不是解壓后的
Error in Read10X() : Gene name or features file missing

第三步：批量讀取成10X對象

Read10X() + CreateSeuratObject()
# 因為Read10X函數(shù)需要對目錄進(jìn)行操作，所以先把目錄名提取出來
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的幫助文檔就知道，它是對列表或向量進(jìn)行循環(huán)，而apply是對數(shù)據(jù)框或矩陣操作)
library(Seurat)
sceList = lapply(folders,function(folder){ 
  CreateSeuratObject(counts = Read10X(folder), 
                     project = folder )
})
# 此時的sceList僅僅是一個堆砌了8個10X對象的集合，下一步就要真正合并起來
> 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

補充

官網(wǎng)的merge教程在：https://satijalab.org/seurat/v3.1/merge_vignette.html

描述了三種情況

第一種： merge兩個seurat對象（原始數(shù)據(jù)）

需要注意的是，組合數(shù)據(jù)時需要注明每個數(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ù)就會出現(xiàn)列名的標(biāo)識
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兩個以上（原始數(shù)據(jù)）

將參數(shù)y 設(shè)成一個向量，就可以指定其他的數(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)情況，只會組合原始數(shù)據(jù)，但如果有的數(shù)據(jù)時標(biāo)準(zhǔn)化之后的呢？

其實可以通過一個參數(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對比：

#################
# 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 . . .

上面的組合多個數(shù)據(jù)就結(jié)束了，接下來是檢查組合后的分群結(jié)果

首先檢查原樣本分群結(jié)果

# 歸一化+標(biāo)準(zhǔn)化（移除了不想要的差異來源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個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 )
image

然后鑒定亞群，看看它們的分群結(jié)果

ElbowPlot(sce.big)
# 官方建議，下游分析時可以多用幾個PCs試試
sce.big <- FindNeighbors(sce.big, dims = 1:20)
# 保持和原文一樣的15個亞群
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')
image
table(sce.big$orig.ident,sce.big@meta.data$RNA_snn_res.0.23)
image

看到，NP有三個主群、G有三個主群、L有三個主群、PI有兩個主群

對比原文的數(shù)據(jù)

它得到了3個NP、3個G、2個L、3個PI，其余的分給了Basal 4群

image

原文：劉小澤學(xué)習(xí)組合多個單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容