学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！
下面是《上海中医药大学研究生》的分享
前面jimmy老师分享了两个祖传的单细胞转录组数据分析代码，非常给力，是标准流程：
- 祖传的单个10x样本的seurat标准代码
- 祖传的单个10x样本的seurat标准代码（人和鼠需要区别对待）
其中有一个环节是需要比较seurat分群以及singleR的分群，这样就可以合理的命名啦。在jimmy老师的督促下，我使用老师的代码处理了GSE135927数据集，直接套用了jimmy老师的标准代码，希望对所有的初学者有帮助！
首先进入GEO可以看到是两个10X的样本：
教程目录大纲如下：
- 1、准备原始分析数据
- 2、创建Seurat对象
- 3、过滤质控
- 4、降维聚类
- 5、clusters细胞类型注释
1、准备原始分析数据
代码语言：javascript复制# 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] TRUEsamples=list.files("GSE135927_RAW/")
samples## [1] "GSM4038043" "GSM4038044"# 是两个文件夹的名字哦2、创建Seurat对象
代码语言：javascript复制# 循环读取两个文件夹下面的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、过滤质控
代码语言：javascript复制#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(C)*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、降维聚类
代码语言：javascript复制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 secondstable(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$seurat_clusters,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 30DimPlot(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 1p=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(p)
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 80p <- 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细胞类型注释
代码语言：javascript复制# 这里的代码是针对上一步分的9个cluster注释celltype，并不是jimmy老师那样的针对每个样本独立注释哦，而且呢，因为我电脑网络问题，采取了云盘下载singleR数据库的方式，如果你的网络好，就直接看jimmy老师的教程哈！refdata <- get(load("ref_Monaco_114s.RData"))
sce_for_SingleR <- GetAssayData(sce, slot="data")
clusters <- sce@meta.data$seurat_clusters
pred.hesc <- SingleR(test = sce_for_SingleR, 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
## 1plotScoreHeatmap(pred.hesc)
#tSNE可视化
celltype = data.frame(ClusterID=rownames(pred.hesc), celltype=pred.hesc$labels, stringsAsFactors = F)
#如下为sce对象注释细胞cluster鉴定结果。
sce@meta.data$celltype = "NA"
#先新增列celltype，值均为NA，然后利用下一行代码循环填充
for(i in 1:nrow(celltype)){
sce@meta.data[which(sce@meta.data$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- 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 cellsp=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(p)
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项目数据哦， 大家如果有自己的项目，赶快用起来吧！
- 祖传的单个10x样本的seurat标准代码
- 祖传的单个10x样本的seurat标准代码（人和鼠需要区别对待）
