手把手教你单细胞测序分析实战

1、数据读取：

以GSE106273为例：

image

下载三个combined文件，修改文件名：Matrix文件改为“matrix.mtx.gz”（这个文件是对基因的表达量），genes文件改为“features.tsv.gz”（这个文件是对基因进行注释），barcode改为“barcode.tsv.gz”，三个文件放到一个文件夹（这里是“SCRNA”文件夹）进行Seurat文件读取（这三个是cellranger count之后的数据格式，很多测序服务公司，一般都会给到cellranger跑完之后的表达矩阵数据，同样在GEO数据上其他单细胞项目提供的数据，所以可以从这三个文件之后开始分析。

library(Seurat)
library(ggplot2)
library(clustree)
library (cowplot)
library (dplyr)
pbmc.data <- Read10x(data.dir = "./SCRNA")
#读取之后创建Seurat对象：
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc", 
                           min.cells = 3, min.features = 200)

#还有另一种数据读入的方法（有时GSE上传的是这种），读取比较简单，但这样读入需要再创建Seurat对象：
a<-read.table("GSE84465_GBM_All_data.csv.gz")
head(rownames(a))
tail(rownames(a),10) #检查一下最后几行基因名称有没有问题

a=a[1:(nrow(a)-5),]#最后5行行名异常，剔除
             
#完了之后创建Seurat对象：
pbmc <- CreateSeuratObject(counts = a, #meta.data = sce_meta, #临床数据读取，这里先不读取
        min.cells = 3, min.features = 50)#最少细胞数不能少于3个，最少基因数不能低于50个）
            #这里没有读入meta.data，关于meta.data的搜索如下：
image.png
image.png
#meta.data文件读取：
b <-read.table("sraRunTable.txt",sep = ",", header = T) #样本信息
table(b$PATIENT_ID)
table(b$Tissue)
table(b$Tissue , b$PATIENT_ID)

#提取b的子集（从meta.data提取需要的列，并使meta.data中包含表达矩阵a文件中的样本名）：
new.b <- b[,c("Plate.ID","Well" ,"Tissue" , "PATIENT_ID")]
new.b$sample <- paste0("X",b$Plate.ID, ".", b$Well)
head(new.b)
identical(colnames(a) ,new.b$sample)#看a的列名跟b的行名是否一致,返回TRUE


#接着还可以继续筛选肿瘤样本
index = which(new.b$Tissue=='Tumor')
group=new.b [index,]
head(group)#筛选肿瘤样本
a.f=a[,index]


#重新创建Seurat对象：
sce_meta=data.frame(PATIENT_ID=new.b$PATIENT_ID, 
                    row.names=new.b$sample)#创建一个新的data.frame
pbmc <- CreateSeuratObject(counts = a, meta.data = sce_meta,
                           min.ce1ls = 3, min.features = 50)
#现在meta.data多了一行PATIENT_ID信息，后续可以用来分组
head(pbmc@meta.data)
table(pbmc@meta.data$PATIENT_ID)

2、Seurat质控+降维+细胞聚类分亚群

读入Seurat对象后，首先进行质控： 查看线粒体百分比，如果都是0，那就不需要过滤

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern ="MT-")
table(pbmc@meta.data$percent.mt)

VlnPlot(object = pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
        group.by= "PATIENT_ID")#画质控图，
image.png

去除表达基因数量少于500或多于6000的细胞，去除线粒体RNA比例高于40%的细胞，如下：

              subset = nFeature_RNA > 500 & nFeature_RNA < 6000 & percent.mt < 20)
                                                                     
#之后可以画相关性散点图：
plot1 <- FeatureScatter(pbmc,feature = "nCount_RNA",
                        feature2 = "percent.mt",
                        group.by = "PATIENT_ID", pt.size = 1.5)#mt百分比很低，所以可以不用画这个

plot2 <- FeatureScatter(pbmc,feature = "nCount_RNA",
                        feature2 = "nFeature_RNA",
                        group.by= "PATIENT_ID", pt.size=1.5)

plot1+plot2 #完成质控，下面开始发现高变异基因。
image.png

下一步鉴定高变异基因，以便用高变异基因进行聚类(基本上不需要改动代码)：

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize",  scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 1500) #变异最大的1500个基因

#输出特征方差图：
top10 <- head(VariableFeatures(pbmc), 10)
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = c(top10,"HLA-DPB1","CCL4"),repel = TRUE)
plot1+plot2
image.png

查看前五个主成分：

a11.genes <- rownames(pbmc)
pbmc<- ScaleData(pbmc, features = a11.genes)
pbmc <- RunPCA(pbmc, features =VariableFeatures(object = pbmc))#这一步比较关键
print(pbmc[["pca"]], dims = 1:5,  nfeatures = 5)

#绘制一些PCA相关的图形：
VizDimLoadings(object = pbmc, dims = 1:2, reduction = "pca", nfeatures = 20) #4个PC 20个基因
DimPlot(pbmc, reduction = "pca",  group.by="PATIENT_ID")
DimHeatmap(pbmc, dims = 1:10, cells = 500, balanced = TRUE, nfeatures = 30, ncol=5)#修改dims =10试试

image.png

接着画关键的2张图：

pbmc <- JackStraw(object = pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(object = pbmc, dims = 1:20)

JackStrawPlot(object = pbmc, dims = 1:20, reduction = "pca")#根据PC的p-value变化最大的PC值，存在主观性
ElbowPlot(pbmc, reduction="pca")#根据E1bowP1ot确定PC数量

JackStrawPlot(object = pbmc, dims = 1:20, reduction = "pca")+ElbowPlot(pbmc, reduction="pca")
image.png

TSNE聚类分析:

pbmc <- FindNeighbors (pbmc,dims = 1:20)
pbmc <- FindClusters(pbmc, resolution = 0.5)#大约3K细胞的单细胞数据集，将resolution参数设置在0.4-1.2之间，数据集增加resolution 值对应增加。不同的resolution值可以获得不同的cluster数目，值越大cluster数目越多，默认值是0.5.
#这一句运行完之后，会自动选择PC数量，本例选择16，及16个细胞cluster。

#tSNE降维
pbmc <- RunTSNE(object = pbmc, dims = 1:20)
DimPlot(pbmc, reduction = "tsne", label = TRUE, pt.size = 1)

image.png

ggplot作图及UMAP聚类：

table(pbmc@meta.data$seurat_clusters)

mydata<- FetchData(pbmc,vars = c("tSNE_1","tSNE_2","PATIENT_ID"))
ggplot(mydata,aes(x= tSNE_1 , y = tSNE_2 ,color = PATIENT_ID)) +geom_point(size = 1 , alpha =1 )+ 
  geom_point(size = 1.5, alpha = 1) + 
  scale_color_manual(values=c("brown3", "chartreuse3", "khaki1", "plum1", "lightskyblue", "darkorange", "hotpink",  "seashell3", "lightblue2"))+
  guides(colour=guide_legend(override.aes=list(size=5)))+ 
  theme(axis.text.y = element_blank(),   
        axis.ticks.y = element_blank(),   
        axis.text.x = element_blank(),  
        axis.ticks.x = element_blank())+  
  theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),
        panel.background = element_rect(fill = "white"))

mydata<- FetchData(pbmc,vars = c("tSNE_1","tSNE_2","seurat_clusters"))
ggplot(mydata,aes(x= tSNE_1 , y = tSNE_2 ,color = seurat_clusters)) +geom_point(size = 1 , alpha =1 )+ 
  geom_point(size = 1.5, alpha = 1) + 
  guides(colour=guide_legend(override.aes=list(size=5)))+ 
  theme(axis.text.y = element_blank(),   
        axis.ticks.y = element_blank(),   
        axis.text.x = element_blank(),   
        axis.ticks.x = element_blank())+  
  theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),
        panel.background = element_rect(fill = "white"))

pbmc <- RunUMAP(pbmc, dims = 1:20)
# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(pbmc, reduction = "umap",
        #cols = c('#FF0033','#CC3399','#993366','#009933', '#0066CC','#003399','#9966CC','#99CC33','#003399'),
        #label = T,
        #label.box = T,
        pt.size = 1
)
image.png

3、singleR细胞族群注释

#细胞注释的准备工作：利用celldex准备注释数据库并加装SingleR包：
#BiocManager: :install("sCRNAseq")
library(celldex)
hpca.se <- HumanPrimaryCellAtlasData() #celldex包里的参考数据集
hpca.se

table(hpca.se@colData@listData$label.main) #列出可注释的细胞数量和种类
table(hpca.se@colData@listData$label.fine) #运行这一句会发现T细胞也有很多种

#SingleR 这个包本身并不会自带数据库啦，而是专门的把数据库文件丢给了celldex包。
#通过celldex得到注释数据库之后就可以加载SingleR：

library(SingleR)
#对数据做简单的转换，将Seurat对象转换成SingleCell支持的对象，这两行都可以
datas= GetAssayData(pbmc, slot = 'data' ) 
datas=as.SingleCellExperiment(pbmc)

cluster_pbmc=pbmc@meta.data$seurat_clusters #这一句比较重要，是上一步做完聚类之后自动生成的(包含那16个聚类)
table(pbmc@meta.data$seurat_clusters)#这一句可以查看每个cluster包含多少个细胞，可以用堆叠柱状图来展示


#接下来通过关键基因将每个cluster进行细胞注释：
pred.hesc<-SingleR(test = datas,ref = hpca.se,labels = hpca.se$label.main,
                   clusters = cluster_pbmc,assay.type.test = "logcounts",assay.type.ref = "logcounts")

table(pred.hesc$labels)

celltype = pred.hesc$labels
celltype #会把每个cluster注释出来
                     
#创建一个data.frame并运行for循环，将其添加到meta信息中，作图：
celltype = data.frame(cluster=rownames(pred.hesc), celltype=pred.hesc$labels)

pbmc@meta.data$celltype= "NA"
                     
for(i in 1: 3511){
                  index=pbmc@meta.data$seurat_clusters[i]
                  pbmc@meta.data$celltype[i]= celltype[index, 2]
                  } #3588个barcode？
#理解一下这个循环：
pbmc@meta.data$seurat_clusters[10]#第十个barcode属于哪个cluster,返回7说明是第7个cluster
celltype[7, 2]#第二列第7行的细胞类型

#作图：
DimPlot(pbmc, group.by = 'celltype', reduction = "tsne", label = TRUE, pt.size = 1)
image.png

4、Seurat细胞亚群差异基因鉴定(即marker基因)：

有的人把4放在3之前，这里放在3之后。

差异表达分析（得到这些基因可以做GO/KEGG/GSEA/PPI/预后/免疫表达等分析）

#找到每一个cluster当中的marker，并且只展示阳性的marker 。
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, 
                               min.pct=0.25,logfc.threshold = 0.25) #计算当前的cluster跟其他所有cluster的差异。这一步会比较慢。
head(pbmc.markers, n = 10) #返回差异基因和属于哪个细胞的cluster，可挑选p值和logfc进一步作为marker基因
write.table(pbmc.markers ,file = 'pbmc.markers.txt', sep = '\t' , row.nameS=F , quote=F)

#选出差异显著的基因：
library("tidyverse")
sig.markers<-pbmc.markers %>% select(gene, everything())%>%subset(p_val_adj<0.05 & abs(pbmc.markers$avg_log2FC)>1)
dim(sig.markers)
head(sig.markers)

上面是每个cluster的DEGs，需要进一步注释到哪一类细胞（有些cluster可能对应同一种细胞）：

  clusterID = as.numeric(sig.markers[i,]$cluster)
  celltypes = celltype[clusterID,2]
  sig.markers$celltype[i] = celltypes
}

#还可以继续提取具体细胞（Astrocyte）的marker基因
ascell = subset(sig.markers, sig.markers$celltype =="Astrocyte")

#运行完之后就多了一列celltype
write.table(sig.markers, file="sigmarkers.xls", sep="\t", row.names=F, quote=F)

top10 <- pbmc.markers %>% group_by(cluster) %>%  top_n(n = 10, wt = avg_log2FC)#提取差异大的基因
top10

#差异基因可视化：
VlnPlot(object = pbmc, features = c("EGFR", "CCL3", "AGXT2L1", "DCN", "GPR17", "MoG" , "ST1"), group.by = "celltype", pt.size = 1)
image.png
#FeaturePlot：在tSNE 或 PCA图中展示重点基因的表达:
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", 
                               "LYZ", "PPBP",  "CD8A"))

FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", 
                               "LYZ", "PPBP", "CD8A"), cols = c("gray", red"),
reduction = "tsne")
image.png

image.png
#细胞注释的可靠性,plotScoreHeatmap显示所有引用细胞类型上所有细胞的得分，这允许用户检查数据集中预测细胞类型的可信度。
plotScoreHeatmap(pred.hesc)
#plotScoreHeatmap显示所有引用细胞类型上所有细胞的得分，这允许用户检查数据集中预测细胞类型的可信度。每个类群/细胞的实际分配标签显示在顶部的颜色栏中。
#我们前面用到了HP（**HumanPrimaryCellAtlasData**）和BP（**BlueprintEncodeData**）两个数据集，所以结果会有2+1个，最上面那个是结合两个注释集的Combined Scores。
#对于下图，关键点是检查分数(scores)在每个类群/细胞中的分布情况。理想情况下，每个类群/细胞(即热图的列)应该有一个明显大于其他细胞的分数，这表明它明确地分配给了单个标签。
image.png

5、细胞拟时序分析（细胞发育轨迹）

根据细胞的高可变基因的变化模拟细胞出现的变化（当前的状态以及发育路线），
要考虑生物学知识，可针对某些亚群进行分析（再做一次分群）,因为有些cluster明显不能从一个发育树分出来。高、低分化状态？

library(monocle)
#先提取细胞亚群（monocle的idents）-这一步也可以不做，对所有细胞分析发育轨迹，参考文件“！！！【生信补给站】跟着NC学pseudotime monocle2 拟时序分析 + 树形图.md”
celltype
Monocyte<-subset(pbmc, idents = c(0,1,8,10))#先提取Monocyte细胞大类，数字0表示cluster的对应数字，从前面celltype数据框中得到。这里的pbmc是注释后的seurat结果（可以从上一步结果中导入）。
table(Monocyte@meta.data$celltype)
Macrophage<-subset(pbmc, idents = c(2,7))#提取Macrophage细胞大类
table(Macrophage@meta.data$celltype)  #发现1,2,7的cluster都是macrophage

data <- as(as.matrix(Macrophage@assays$RNA@counts),"sparseMatrix")#用到counts矩阵

#创建对象准备做拟时序分析：
pd <- new ( 'AnnotatedDataFrame', data = Macrophage@meta.data)
fData <- data.frame(gene_short_name = row.names (data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)

#构建S4对象，CellDataSet
HSMM <- newCellDataSet(data, phenoData = pd, featureData = fd, 
                       #lowerDetectionLimit = 0.5,
                       expressionFamily = negbinomial.size())
#注意expressionFamily的选择：
#单细胞稀疏矩阵的话用negbinomial.size()，如果是UMI的话 不要标准化；
#FPKM值用tobit()；
#logFPKM值用gaussianff()

HSMM<- estimateSizeFactors(HSMM) #估计size factors
HSMM<- estimateDispersions(HSMM) #估计dispersions 
HSMM<- detectGenes(HSMM, min_expr = 0.1) #过滤低质量细胞
#print(head(fData(HSMM)))

disp_table <- dispersionTable(HSMM)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
HSMM_myo <- setOrderingFilter(HSMM, unsup_clustering_genes$gene_id)
plot_ordering_genes(HSMM)
#上面这六句不用改，官方网站提供的
#这里选择基因的方式有很多，说明文档中建议以下4种选择基因的方式(此处大概是第三种)
#（1）选择发育差异表达基因
#（2）选择clusters差异表达基因
#（3）选择离散程度高的基因
#（4）自定义发育marker基因
image.png
#首先降维处理
HSMM_myo <- reduceDimension(HSMM_myo, max_components= 2, num_dim = 20, #画在平面上，所以填2
                            method  = 'DORTree', verbose = T) #这一步非常慢，而且只能用单核GPU

HSMM_myo <-orderCells(HSMM_myo)
                            
plot_cell_trajectory(HSMM_myo)
image.png
##以state进行着色
colour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00",  
         "#808000","#FF00FF","#FA8072","#7B68EE","#9400D3","#800080","#A0522D","#D2B48C","#D2691E","#87CEEB","#40E0D0","#5F9EA0",
         "#FF1493","#0000CD","#008B8B","#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4","#32CD32","#F0E68C","#FFFFE0","#EE82EE",
         "#FF6347","#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887")

plot_cell_trajectory(HSMM_myo, color_by = "State") + scale_color_manual(values = colour)
plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime", ce1l_size = 0.75)#根据拟时间值着色
                   ###绘制cluster的分面图
plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters", ce1l_size = 0.75)#如果有Seurat生的rds文件
image.png

image.png

image.png
#可以使用分面单独查看各单一celltype的时序状态                      
plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters") +
  facet_wrap(~seurat_clusters, nrow = 2) #设置几行几列展示
image.png
#添加“树形图”
p1 <- plot_cell_trajectory(HSMM_myo, x = 1, y = 2, color_by = "seurat_clusters") + 
  theme(legend.position='none',panel.border = element_blank()) + #去掉第一个的legend
  scale_color_manual(values = colour) 
p2 <- plot_complex_cell_trajectory(HSMM_myo, x = 1, y = 2,
                                   color_by = "seurat_clusters")+
  scale_color_manual(values = colour) +
  theme(legend.title = element_blank()) 

p1 / p2
image.png

更多期待待复现：

CellChat细胞通讯分析（一） (qq.com)

文献数据复现——原发性和转移性肝细胞癌多细胞生态系统的单细胞图谱 (qq.com)

使用cytoTRACE评估不同单细胞亚群的分化潜能 - 知乎 (zhihu.com)

scRNA单细胞拟时分析工具介绍——CytoTRACE (qq.com)

scRNA单细胞拟时分析工具介绍——CytoTRACE (qq.com)

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 212,542评论 6赞 493
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,596评论 3赞 385
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 158,021评论 0赞 348
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,682评论 1赞 284
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 65,792评论 6赞 386
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,985评论 1赞 291
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,107评论 3赞 410
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,845评论 0赞 268
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,299评论 1赞 303
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,612评论 2赞 327
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,747评论 1赞 341
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,441评论 4赞 333
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 40,072评论 3赞 317
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,828评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,069评论 1赞 267
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,545评论 2赞 362
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,658评论 2赞 350

推荐阅读更多精彩内容