seurat 与geo geo单细胞导入seurat实战 数据下载 数据分析 scrna-seq

前面jimmy老师分享了两个祖传的单细胞转录组数据分析代码，非常给力，是标准流程：

祖传的单个10x样本的seurat标准代码
祖传的单个10x样本的seurat标准代码（人和鼠需要区别对待）
其中有一个环节是需要比较seurat分群以及singleR的分群，这样就可以合理的命名啦。在jimmy老师的督促下，我使用老师的代码处理了GSE135927数据集，直接套用了jimmy老师的标准代码，希望对所有的初学者有帮助！

首先进入GEO可以看到是两个10X的样本：

教程目录大纲如下：

1、准备原始分析数据
2、创建Seurat对象
3、过滤质控
4、降维聚类
5、clusters细胞类型注释
1、准备原始分析数据

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE135927

两个样本，共6个文件。需要分到两个文件夹里，并且重命名

fs=list.files(’./GSE135927_RAW/’,’^GSM’)
fs

自行下载GSE135927数据集的GSE135927_RAW压缩包并且解压哦，这样上面的代码就可以运行啦

然后获取两个样本信息，因为是批量，所以下面的代码可能不好理解，需要熟练掌握R语言哦

samples=str_split(fs,’_’,simplify = T)[,1]

lapply(unique(samples),function(x){
y=fs[grepl(x,fs)]
folder=paste0(“GSE135927_RAW/”, str_split(y[1],’_’,simplify = T)[,1])
dir.create(folder,recursive = T)
#为每个样本创建子文件夹
file.rename(paste0(“GSE135927_RAW/”,y[1]),file.path(folder,“barcodes.tsv.gz”))
#重命名文件，并移动到相应的子文件夹里
file.rename(paste0(“GSE135927_RAW/”,y[2]),file.path(folder,“features.tsv.gz”))
file.rename(paste0(“GSE135927_RAW/”,y[3]),file.path(folder,“matrix.mtx.gz”))
})

[[1]]

[1] TRUE

[[2]]

[1] TRUE

samples=list.files(“GSE135927_RAW/”)
samples

[1] “GSM4038043” “GSM4038044”

是两个文件夹的名字哦

2、创建Seurat对象

循环读取两个文件夹下面的10x的的3个文件

sceList = lapply(samples,function(pro){
folder=file.path(“GSE135927_RAW/”,pro)
CreateSeuratObject(counts = Read10X(folder),
project = pro )
})
sce.big <- merge(sceList[[1]],
y = c(sceList[[2]]),
add.cell.ids = samples,
project = “ls_12”)
sce.big

An object of class Seurat

33538 features across 2042 samples within 1 assay

Active assay: RNA (33538 features, 0 variable features)

table(sce.big$orig.ident)

GSM4038043 GSM4038044

1359 683

#save(sce.big,file = ‘sce.big.merge.ls_12.Rdata’)
3、过滤质控
#raw_sce=get(load(file = ‘sce.big.merge.ls_12.Rdata’))
#线粒体基因
raw_sce=sce.big
rownames(raw_sce)[grepl(’^mt-’,rownames(raw_sce),ignore.case = T)]

[1] “MT-ND1” “MT-ND2” “MT-CO1” “MT-CO2” “MT-ATP8” “MT-ATP6”

#核糖体蛋白基因
rownames(raw_sce)[grepl(’^Rp[sl]’,rownames(raw_sce),ignore.case = T)]

[1] “RPL22” “RPL11” “RPS6KA1” “RPS8”

[5] “RPL5” “RPS27” “RPS6KC1” “RPS7”

[9] “RPS27A” “RPL31” “RPL37A” “RPL32”

[13] “RPL15” “RPSA” “RPL14” “RPL29”

#计算上述两种基因在所有细胞中的比例
raw_sce[[“percent.mt”]] <- PercentageFeatureSet(raw_sce, pattern = “^MT-”)

rb.genes <- rownames(raw_sce)[grep("^RP[SL]",rownames(raw_sce))]
C<-GetAssayData(object = raw_sce, slot = “counts”)
percent.ribo <- Matrix::colSums(C[rb.genes,])/Matrix::colSums©*100
raw_sce <- AddMetaData(raw_sce, percent.ribo, col.name = “percent.ribo”)

plot1 <- FeatureScatter(raw_sce, feature1 = “nCount_RNA”, feature2 = “percent.mt”)
plot2 <- FeatureScatter(raw_sce, feature1 = “nCount_RNA”, feature2 = “nFeature_RNA”)
CombinePlots(plots = list(plot1, plot2))

VlnPlot(raw_sce, features = c(“percent.ribo”, “percent.mt”), ncol = 2)

VlnPlot(raw_sce, features = c(“nFeature_RNA”, “nCount_RNA”), ncol = 2)

VlnPlot(raw_sce, features = c(“percent.ribo”, “nCount_RNA”), ncol = 2)

raw_sce #2042个

An object of class Seurat

33538 features across 2042 samples within 1 assay

Active assay: RNA (33538 features, 0 variable features)

#按照三个指标过滤细胞
raw_sce1 <- subset(raw_sce, subset = nFeature_RNA > 200 & nCount_RNA > 1000 & percent.mt < 20)
raw_sce1 #1887个

An object of class Seurat

33538 features across 1887 samples within 1 assay

Active assay: RNA (33538 features, 0 variable features)

4、降维聚类
sce=raw_sce1
#counts归一化
sce <- NormalizeData(sce, normalization.method = “LogNormalize”,
scale.factor = 10000)
GetAssay(sce,assay = “RNA”)

Assay data with 33538 features for 1887 cells

First 10 features:

MIR1302-2HG, FAM138A, OR4F5, AL627309.1, AL627309.3, AL627309.2,

AL627309.4, AL732372.1, OR4F29, AC114498.1

#前2000个高变feature RNA
sce <- FindVariableFeatures(sce,
selection.method = “vst”, nfeatures = 2000)

中心化，为下一步PCA做准备

sce <- ScaleData(sce)
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce))
#看看前12个主成分
DimHeatmap(sce, dims = 1:12, cells = 100, balanced = TRUE)

#判断最终选取的主成分数，这里我判断18个
ElbowPlot(sce)

sce <- FindNeighbors(sce, dims = 1:18)
sce <- FindClusters(sce, resolution = 0.9)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 1887

Number of edges: 67622

Running Louvain algorithm…

Maximum modularity in 10 random starts: 0.8140

Number of communities: 9

Elapsed time: 0 seconds

table(sce@meta.data$RNA_snn_res.0.9)

0 1 2 3 4 5 6 7 8

325 303 275 273 214 145 139 133 80

分成9个cluster

#tSNE可视化
set.seed(123)
sce <- RunTSNE(object = sce, dims = 1:18, do.fast = TRUE)
DimPlot(sce,reduction = “tsne”,label=T)

在两个样本里9个cluster的分布情况

table(sce@meta.data$seura{t}_{c}lusters,sce@meta.data$orig.ident)

GSM4038043 GSM4038044

0 108 217

1 75 228

2 237 38

3 271 2

4 201 13

5 92 53

6 137 2

7 74 59

8 50 30

DimPlot(sce,reduction = “tsne”,label=T,split.by =‘orig.ident’)

#save(dat,file=‘merge_for_tSNE.pos.Rdata’)
#load(file='merge_for_tSNE.pos.Rdata)

#再尝试用ggplot2可视化上图
phe=data.frame(cell=rownames(sce@meta.data),
cluster =sce@meta.data$seurat_clusters)
tsne_pos=Embeddings(sce,‘tsne’)
dat=cbind(tsne_pos,phe)
head(dat)

tSNE_1 tSNE_2

GSM4038043_AAACCTGAGGGCTCTC-1 -13.779792 1.304764

GSM4038043_AAACCTGAGTGAACGC-1 20.375308 5.068255

GSM4038043_AAACCTGGTATGAATG-1 -25.223356 -3.429902

GSM4038043_AAACCTGTCTGCGACG-1 6.194084 14.146905

GSM4038043_AAACGGGAGGTGCAAC-1 5.695203 -28.687411

GSM4038043_AAACGGGTCGCTAGCG-1 5.739440 -10.725417

cell cluster

GSM4038043_AAACCTGAGGGCTCTC-1 GSM4038043_AAACCTGAGGGCTCTC-1 0

GSM4038043_AAACCTGAGTGAACGC-1 GSM4038043_AAACCTGAGTGAACGC-1 2

GSM4038043_AAACCTGGTATGAATG-1 GSM4038043_AAACCTGGTATGAATG-1 0

GSM4038043_AAACCTGTCTGCGACG-1 GSM4038043_AAACCTGTCTGCGACG-1 3

GSM4038043_AAACGGGAGGTGCAAC-1 GSM4038043_AAACGGGAGGTGCAAC-1 5

GSM4038043_AAACGGGTCGCTAGCG-1 GSM4038043_AAACGGGTCGCTAGCG-1 1

p=ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=cluster))+geom_point(size=0.95)
p=p+stat_ellipse(data=dat,aes(x=tSNE_1,y=tSNE_2,fill=cluster,color=cluster),
geom = “polygon”,alpha=0.2,level=0.9,type=“t”,linetype =2,show.legend = F)+coord_fixed()
print§

theme= theme(panel.grid =element_blank()) + ## 删去网格
theme(panel.border = element_blank(),panel.background = element_blank()) + ## 删去外层边框
theme(axis.line = element_line(size=1, colour = “black”))
p+theme+guides(colour = guide_legend(override.aes = list(size=5)))

#寻找每个cluster的高变代表基因，并选取前5个，进行可视化
table(sce@meta.data$seurat_clusters)

0 1 2 3 4 5 6 7 8

325 303 275 273 214 145 139 133 80

p <- list()
for( i in unique(sce@meta.data$seurat_clusters) ){
markers_df <- FindMarkers(object = sce, ident.1 = i, min.pct = 0.25)
print(x = head(markers_df))
markers_genes = rownames(head(x = markers_df, n = 4))
p1 <- VlnPlot(object = sce, features =markers_genes,log =T,ncol = 2)
p[[i]][[1]] <- p1
p2 <- FeaturePlot(object = sce, features=markers_genes,ncol = 2)
p[[i]][[2]] <- p2
}

p_val avg_logFC pct.1 pct.2 p_val_adj

NRP2 1.051573e-129 1.1153915 0.585 0.059 3.526767e-125

RPL32 1.930788e-125 0.7562735 1.000 0.976 6.475476e-121

RPS3A 2.407233e-125 0.7527394 1.000 0.985 8.073378e-121

EEF1A1 1.509065e-122 0.7042242 1.000 0.997 5.061102e-118

RPL18A 6.076286e-119 0.7034375 1.000 0.983 2.037865e-114

RPS14 8.819300e-115 0.5999948 1.000 0.981 2.957817e-110

p[[1]][[2]]

p[[2]][[1]]

#查阅所有的marker基因，可用于人工注释cell type
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE, min.pct = 0.25,
thresh.use = 0.25)
DT::datatable(sce.markers)

#热图可视化每个cluster的marker基因表达差异
top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_logFC)
DoHeatmap(sce,top10$gene,size=3)

#save(sce,sce.markers,file = ‘sce.output.merge.ls_12.Rdata’)
5、clusters细胞类型注释

这里的代码是针对上一步分的9个cluster注释celltype，并不是jimmy老师那样的针对每个样本独立注释哦，而且呢，因为我电脑网络问题，采取了云盘下载singleR数据库的方式，如果你的网络好，就直接看jimmy老师的教程哈！

refdata <- get(load(“ref_Monaco_114s.RData”))
sce_for_SingleR <- GetAssayData(sce, slot=“data”)
clusters <- sce@meta.data$seura{t}_{c}lusterspred.hesc<-SingleR(test=sc{e}_{f}o{r}_{S}ingleR,ref=refdata,labels=refdata$label.fine,
#因为样本主要为免疫细胞（而不是全部细胞），因此设置为label.fine
method = “cluster”, clusters = clusters,
#这里我们为上一步分的9个cluster注释celltype
assay.type.test = “logcounts”, assay.type.ref = “logcounts”)
table(pred.hesc$labels)

Central memory CD8 T cells Effector memory CD8 T cells

3 2

T regulatory cells Th17 cells

2 1

Th2 cells

1

plotScoreHeatmap(pred.hesc)

#tSNE可视化
celltype = data.frame(ClusterID=rownames(pred.hesc), celltype=pred.hesc$labels, stringsAsFactors = F)
#如下为sce对象注释细胞cluster鉴定结果。

sce@meta.dataKaTeX parse error: Expected 'EOF', got '#' at position 17: …elltype = "NA" #̲先新增列celltype，值均…seurat_clusters == celltype$ClusterID[i]){,}^{\prime }celltyp{e}^{\prime }]<-celltype$celltype[i]}

DimPlot(sce, group.by=“celltype”, label=T, label.size=5, reduction=‘tsne’)

#利用ggplot2再可视化上图
phe=data.frame(cell=rownames(sce@meta.data),
cluster =clusters,
celltype=sce@meta.data$celltype)
tsne_pos=Embeddings(sce,‘tsne’)
dat=cbind(tsne_pos,phe)
head(dat)

tSNE_1 tSNE_2

GSM4038043_AAACCTGAGGGCTCTC-1 -13.779792 1.304764

GSM4038043_AAACCTGAGTGAACGC-1 20.375308 5.068255

GSM4038043_AAACCTGGTATGAATG-1 -25.223356 -3.429902

GSM4038043_AAACCTGTCTGCGACG-1 6.194084 14.146905

GSM4038043_AAACGGGAGGTGCAAC-1 5.695203 -28.687411

GSM4038043_AAACGGGTCGCTAGCG-1 5.739440 -10.725417

cell cluster

GSM4038043_AAACCTGAGGGCTCTC-1 GSM4038043_AAACCTGAGGGCTCTC-1 0

GSM4038043_AAACCTGAGTGAACGC-1 GSM4038043_AAACCTGAGTGAACGC-1 2

GSM4038043_AAACCTGGTATGAATG-1 GSM4038043_AAACCTGGTATGAATG-1 0

GSM4038043_AAACCTGTCTGCGACG-1 GSM4038043_AAACCTGTCTGCGACG-1 3

GSM4038043_AAACGGGAGGTGCAAC-1 GSM4038043_AAACGGGAGGTGCAAC-1 5

GSM4038043_AAACGGGTCGCTAGCG-1 GSM4038043_AAACGGGTCGCTAGCG-1 1

celltype

GSM4038043_AAACCTGAGGGCTCTC-1 Th2 cells

GSM4038043_AAACCTGAGTGAACGC-1 Effector memory CD8 T cells

GSM4038043_AAACCTGGTATGAATG-1 Th2 cells

GSM4038043_AAACCTGTCTGCGACG-1 T regulatory cells

GSM4038043_AAACGGGAGGTGCAAC-1 Central memory CD8 T cells

GSM4038043_AAACGGGTCGCTAGCG-1 Central memory CD8 T cells

p=ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=celltype))+geom_point(size=0.95)
p=p+stat_ellipse(data=dat,aes(x=tSNE_1,y=tSNE_2,fill=celltype,color=celltype),
geom = “polygon”,alpha=0.2,level=0.9,type=“t”,linetype = 2,show.legend = F)+coord_fixed()
print§

theme= theme(panel.grid =element_blank()) + ## 删去网格
theme(panel.border = element_blank(),panel.background = element_blank()) + ## 删去外层边框
theme(axis.line = element_line(size=1, colour = “black”))
p+theme+guides(colour = guide_legend(override.aes = list(size=5)))

理论上jimmy老师的代码是可以适用于绝大部分10x项目数据哦， 大家如果有自己的项目，赶快用起来吧！