分享是一种态度
最近在学单细胞测序,学习了B站生信技能树的单细胞教程,https://www.bilibili.com/video/BV1dt411Y7nn结合jimmy老师的代码,用GEO上的单细胞测序数据,做了一下分析。
我们选择的数据已经发表的文章题目是“Single cell RNA sequencing of human liver reveals distinct intrahepatic macrophage populations”,2018年发表在nature communications上,数据存在GSE115469。
在https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE115469下载数据文件”GSE115469_Data.csv.gz“
#一键清空
rm(list=ls())
options(stringsAsFactors = F)
#读入并处理数据
a=read.csv('GSE115469_Data.csv')
rownames(a)=a[,1] #将第一列基因名字作为列名
a=a[,-1]#去除第一列
dim(a)
[1] 20007 8445
#8445个细胞,总共有20007个基因
注意,因为这篇文章作者已经上传了过滤后的矩阵,所以在 “CreateSeuratObject”函数中我们不需要用“min.cells =” 和“min.features = ”来过滤基因和细胞,关于数据如何处理的,详见https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM3178782%5D
library(Seurat)
raw_sce <- CreateSeuratObject(counts = a)#
save(raw_sce,file="raw_sce.Rdata")
load("raw_sce.Rdata")
raw_sce
#An object of class Seurat
#20007 features across 8444 samples within 1 assay
#Active assay: RNA (20007 features, 0 variable features)
首先我们要找出线粒体基因和线粒体核糖体基因。
rownames(raw_sce)[grepl('^mt-',rownames(raw_sce),ignore.case = T)]
rownames(raw_sce)[grepl('^Rp[sl]',rownames(raw_sce),ignore.case = T)]
#其实在Seurat包中,PercentageFeatureSet函数可以用来提取并计算线粒体基因比例。
raw_sce[["percent.mt"]] <- PercentageFeatureSet(raw_sce, pattern = "^MT-")
#提取并计算线粒体核糖体基因比例,为raw_sce添加线粒体核糖体基因注释信息。
rb.genes <- rownames(raw_sce)[grep("^RP[SL]",rownames(raw_sce),ignore.case = T)]
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)
接下来我们可以通过subset函数过滤基因和样本,但是由于作者已经将过滤好了,所以我们不进行下一步subset操作。按照一般情况,都是设置一个阈值(10%)来过滤线粒体基因,但是由于肝细胞的特殊性,作者将阈值设置为50%,关于肝细胞阈值的选择,作者是这么解释的。
#过滤的代码如下:
raw_sce <- subset(raw_sce, subset = nFeature_RNA >100 & nCount_RNA > 1000 & percent.mt < 50)
#100,1000,50这3个数值要根据实际情况调整
sce=raw_sce
sce<- NormalizeData(sce, normalization.method = "LogNormalize",
scale.factor = 10000 )
GetAssay(sce,assay = "RNA")
#Assay data with 20007 features for 8444 cells
#First 10 features:RP11-34P13.7, FO538757.2, AP006222.2, RP4-669L17.10, RP5-857K21.4, RP11-206L10.9,LINC00115, FAM41C, RP11-54O7.1, RP11-54O7.3
nfeatures选择多少取决于实际情况,这里我们选择5000
sce <- FindVariableFeatures(sce,
selection.method = "vst", nfeatures = 5000)
由于保留的线粒体基因较多,因此我这里选择矫正的参数为"percent.mt","nCount_RNA"和"percent.ribo"试试看
sce <- ScaleData(sce,vars.to.regress=c("percent.mt","nCount_RNA","percent.ribo"))
#运行scaleData之前要进行NormalizeData
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce))
DimHeatmap(sce, dims = 1:12, cells = 100, balanced = TRUE)
#展示前12个主成分热图
判断选择的主成分数目,我选择19个
ElbowPlot(sce)
#调整参数resolution = 0.55,得到20个cluster
sce <- FindNeighbors(sce, dims = 1:19)
sce <- FindClusters(sce, resolution = 0.55)
#看看每个cluster的细胞数目
table(sce@meta.data$RNA_snn_res.0.55)
> 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
1313 1047 768 682 634 517 517 514 489 484 349 286 167 137 118 107 97 93
18 19
90 35
set.seed(123)
sce <- RunTSNE(object = sce, dims = 1:19, do.fast = TRUE)
DimPlot(sce,reduction = "tsne",label=T)
看看文献的tSNE结果
给图片添加meta信息,结果貌似跟文献的图片结果差不多。P3TLH样本的细胞(绿色部分)基本独立存在,由3个cluster组成
phe=data.frame(cell=rownames(sce@meta.data),
cluster =sce@meta.data$seurat_clusters)
head(phe)
table(phe$cluster)
tsne_pos=Embeddings(sce,'tsne')
DimPlot(sce,reduction = "tsne",group.by ='orig.ident')
再看看文献的图
DimPlot(sce,reduction = "tsne",label=T,split.by ='orig.ident')
#ggplot可视化
dat=cbind(tsne_pos,phe)
library(ggplot2)
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()
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=p+theme+guides(colour = guide_legend(override.aes = list(size=5)))
print(p)
#对marker基因进行可视化,这一步会把所有cluster的marker基因展示出来。这步骤运行时间较久
pro='first'
table(sce@meta.data$seurat_clusters)
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 = 5))
VlnPlot(object = sce, features =markers_genes,log =T )
ggsave(filename=paste0(pro,'_VlnPlot_subcluster_',i,'_sce.markers_heatmap.pdf'))
FeaturePlot(object = sce, features=markers_genes )
ggsave(filename=paste0(pro,'_FeaturePlot_subcluster_',i,'_sce.markers_heatmap.pdf'))
}
#提取marke基因,min.ct:设定研究的基因在两组细胞中的最小占比。
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE, min.pct = 0.25,
thresh.use = 0.25)#only.pos = TRUE:只返回positive的marker基因
print(x = head(sce.markers))
> p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
GSTA1 0 1.283148 0.865 0.199 0 0 GSTA1
APOM 0 1.222117 0.854 0.185 0 0 APOM
ANG 0 1.154127 0.986 0.374 0 0 ANG
TAT-AS1 0 1.143226 0.953 0.237 0 0 TAT-AS1
UGT2B7 0 1.130114 0.841 0.182 0 0 UGT2B7
CYP3A4 0 1.126983 0.776 0.225 0 0 CYP3A4
DT::datatable(sce.markers)
write.csv(sce.markers,file=paste0('_sce.markers_tsne.csv'))
#可视化marker基因,我们选择ALB基因可视化看看
markers_genes = rownames(sce.markers[c("ALB"),])
VlnPlot(object = sce, features =markers_genes,log=T,ncol=2)
FeaturePlot(object = sce,
features =markers_genes,
ncol=2)
#热图可视化marker基因的表达差异
library(dplyr)
top5 <- sce.markers %>% group_by(cluster) %>% top_n(5, avg_logFC)
DoHeatmap(sce,top5$gene,size=2)
ggsave(filename=paste0('_sce.markers_——stne_heatmap.pdf'))
save(sce,sce.markers,file = 'last_sce_tsne.Rdata')
sce_cc <- CellCycleScoring(sce, s.features = cc.genes$s.genes, g2m.features =cc.genes$ g2m.genes, set.ident = F)
head(sce_cc[[]])
DimPlot(sce_cc,reduction = "tsne",group.by ='Phase',cols=c("yellow","#458B74","#483D8B"))
再跟文献的图对比一下,可以看到中间那团大细胞群体基本都在G1期,跟我们的结果还是相似的。