单细胞分析的经典包——Seurat包有自己的官方教程，跟着教程过一遍可以get大致过程（https://satijalab.org/seurat/archive/v3.1/pbmc3k_tutorial.html）
在我看来最难的部分其实是第一步的数据准备阶段，如果GEO数据库可以下载到标准的10X数据（包含三个文件（barcodes.tsv/genes.tsv/matrix.mtx）），那直接用Seurat包自带的Read10X函数就能一步解决。但是现实是GEO数据库上传的数据五花八门，或者我们需要整合各种数据进行我们自己的生信再分析，就需要我们将数据整理成Seurat包能识别的标准格式。以下介绍几种情况：
写在最前面
setwd()
1.有直接的标准10X数据（喜大普奔）
解压缩后可以得到三个文件（barcodes.tsv/genes.tsv/matrix.mtx），文件名修改到一模一样
例如：GSE106273
下载后三个文件，解压缩后文件名改为barcodes.tsv、genes.tsv、matrix.mtx（一个字也不差）
image
image
pbmc.data <- Read10X(data.dir = "C:/Users/fhche/Desktop/GSE106273")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc", min.cells = 3, min.features = 200)
head(pbmc@meta.data)
2.多个10X数据可以用merge函数合并
例如：GSE135927，只有一个raw data能下载
image.png
下载后整理成GSM4038043、GSM4038044两个文件夹，分别含有barcodes.tsv、genes.tsv、matrix.mtx三个文件
image.png
image.png
GSM4038043<- Read10X(data.dir = "C:/Users/fhche/Desktop/GSE135927/GSM4038043")
pbmc1 <- CreateSeuratObject(counts = GSM4038043,
min.cells = 3,
min.features = 200)
GSM4038044<- Read10X(data.dir = "C:/Users/fhche/Desktop/GSE135927/GSM4038044")
pbmc2 <- CreateSeuratObject(counts = GSM4038044,
min.cells = 3,
min.features = 200)
head(pbmc2@meta.data)
pbmc = merge(pbmc1, pbmc2,
add.cell.ids = c("GSM4038043", "GSM4038044"),
merge.data = TRUE)
as.data.frame(pbmc@assays$RNA@counts[1:10, 1:2])
head(pbmc@meta.data)
3.GEO里只有矩阵数据
例如：GSE157703
image.png
解压后
image.png
library(data.table)
library("R.utils")
pca1 <- fread("GSM4773521_PCa1_gene_counts_matrix.txt.gz",
data.table = F)
pca1[1:4,1:4]
d1=pca1[,-1]
rownames(d1)=pca1[,1]
pca2 <- fread("GSM4773522_PCa2_gene_counts_matrix.txt.gz",
data.table = F)
pca2[1:4,1:4]
d2=pca2[,-1]
rownames(d2)=pca2[,1]
pbmc1 <- CreateSeuratObject(counts = d1,
min.cells = 3,
min.features = 200,
project = "pca1")
pbmc2 <- CreateSeuratObject(counts = d2,
min.cells = 3,
min.features = 200,
project = "pca2")
pbmc = merge(pbmc1, pbmc2, add.cell.ids = c("pca1", "pca2"),
project = "PCA", merge.data = TRUE)
as.data.frame(pbmc@assays$RNA@counts[1:10, 1:2])
head(pbmc@meta.data)
4.数据需要筛选后再构建Seurat矩阵
例如：GSE84465，需要筛选其中的Tumor细胞数据进一步分析
image.png
下载第一个
image.png
下载Metadata的SraRunTable.txt
a<-read.table("GSE84465_GBM_All_data.csv.gz")
head(rownames(a))
tail(rownames(a),10)
a=a[1:(nrow(a)-5),] #最后5行行名异常，剔除
b<-read.table("SraRunTable.txt",sep = ",",header = T) #样本信息
table(b$patient_id)
table(b$Tissue)
table(b$Tissue, b$patient_id)
new.b <- b[,c("plate_id","Well","Tissue","patient_id")]
new.b$sample <- paste0("X",b.group$plate_id,".",b.group$Well)
head(new.b)
identical(colnames(a),new.b$sample)
#筛选肿瘤细胞
index<-which(new.b$Tissue=="Tumor")
group<-new.b[index,] #筛选的是行
a.filt<-a[,index] #筛选的是列
dim(a.filt)
identical(colnames(a.filt),group$sample)
#构建Seurat对象
library("Seurat")
datagroup=data.frame(Patient_ID=group$patient_id,
group=new.b[index,3],
row.names = group$sample)
pbmc <- CreateSeuratObject(counts = a.filt,
meta.data = datagroup,
min.cells = 3,
min.features = 50)
head(pbmc@meta.data)
最后记得
saveRDS(pbmc,file="pbmc_raw.rds") #用于后续分析
