單細(xì)胞分析的經(jīng)典包——Seurat包有自己的官方教程，跟著教程過一遍可以get大致過程（https://satijalab.org/seurat/archive/v3.1/pbmc3k_tutorial.html）

在我看來最難的部分其實(shí)是第一步的數(shù)據(jù)準(zhǔn)備階段，如果GEO數(shù)據(jù)庫可以下載到標(biāo)準(zhǔn)的10X數(shù)據(jù)（包含三個(gè)文件（barcodes.tsv/genes.tsv/matrix.mtx）），那直接用Seurat包自帶的Read10X函數(shù)就能一步解決。但是現(xiàn)實(shí)是GEO數(shù)據(jù)庫上傳的數(shù)據(jù)五花八門，或者我們需要整合各種數(shù)據(jù)進(jìn)行我們自己的生信再分析，就需要我們將數(shù)據(jù)整理成Seurat包能識(shí)別的標(biāo)準(zhǔn)格式。以下介紹幾種情況：
寫在最前面

setwd()

1.有直接的標(biāo)準(zhǔn)10X數(shù)據(jù)（喜大普奔）
解壓縮后可以得到三個(gè)文件（barcodes.tsv/genes.tsv/matrix.mtx），文件名修改到一模一樣
例如：GSE106273
下載后三個(gè)文件，解壓縮后文件名改為barcodes.tsv、genes.tsv、matrix.mtx（一個(gè)字也不差）

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.多個(gè)10X數(shù)據(jù)可以用merge函數(shù)合并
例如：GSE135927，只有一個(gè)raw data能下載

image.png

下載后整理成GSM4038043、GSM4038044兩個(gè)文件夾，分別含有barcodes.tsv、genes.tsv、matrix.mtx三個(gè)文件

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里只有矩陣數(shù)據(jù)
例如：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.數(shù)據(jù)需要篩選后再構(gòu)建Seurat矩陣
例如：GSE84465，需要篩選其中的Tumor細(xì)胞數(shù)據(jù)進(jìn)一步分析

image.png

下載第一個(gè)

image.png

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)

#篩選腫瘤細(xì)胞
index<-which(new.b$Tissue=="Tumor")
group<-new.b[index,]   #篩選的是行
a.filt<-a[,index]        #篩選的是列
dim(a.filt)
identical(colnames(a.filt),group$sample)

#構(gòu)建Seurat對(duì)象
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")  #用于后續(xù)分析