一篇单细胞数据挖掘文章的图表复现

前言

看到一篇单细胞数据挖掘的文章，题为：Establishment of a Prognostic Model of Lung Adenocarcinoma Based on Tumor Heterogeneity

遂打算拿里面的数据跑一跑，这个数据可以在GSE117570的补充文件里直接下载到。

1.批量读取数据

虽然不是标准10X的三个文件，但也可以搞，直接读取为数据框，转换为矩阵，自行创建Seurat对象就可以啦。

rm(list = ls())
library(stringr)
library(Seurat)
library(dplyr)
f = dir("GSE117570_RAW/");f
## [1] "GSM3304007_P1_Tumor_processed_data.txt.gz"
## [2] "GSM3304011_P3_Tumor_processed_data.txt.gz"
## [3] "GSM3304013_P4_Tumor_processed_data.txt.gz"
s = str_split(f,"_",simplify = T)[,1];s
## [1] "GSM3304007" "GSM3304011" "GSM3304013"

原本是8个文件来着，这篇文章是只拿了其中3个。

批量读取，顺便把线粒体基因超过5%的细胞过滤掉了。

scelist = list()
for(i in 1:length(f)){
  tmp = read.table(paste0("GSE117570_RAW/",f[[i]]),
                 check.names = F)
  tmp = as.matrix(tmp)

  scelist[[i]] <- CreateSeuratObject(counts = tmp, 
                           project = s[[i]], 
                           min.cells = 3, 
                           min.features = 50)
  print(dim(scelist[[i]]))
  scelist[[i]][["percent.mt"]] <- PercentageFeatureSet(scelist[[i]], pattern = "^MT-")
  scelist[[i]] <- subset(scelist[[i]], subset = percent.mt < 5)
  print(dim(scelist[[i]]))
}
## [1] 6611 1832
## [1] 6611 1466
## [1] 3211  328
## [1] 3211  116
## [1] 7492 1423
## [1] 7492  281
names(scelist)  = s

过滤的还挺狠的奥

2.多样本的整合

然后用CCA方法完成多个样本的整合，我把nfeatures参数设置的大了一些，不然整合完了只剩下2000个基因。

for (i in 1:length(scelist)) {
  scelist[[i]] <- NormalizeData(scelist[[i]], verbose = FALSE)
  scelist[[i]] <- FindVariableFeatures(scelist[[i]],verbose = FALSE,nfeatures = 8000)
}
features <- SelectIntegrationFeatures(object.list = scelist,nfeatures = 8000)
sce.anchors <- FindIntegrationAnchors(object.list = scelist,anchor.features = features)
sce.integrated <- IntegrateData(anchorset = sce.anchors)

DefaultAssay(sce.integrated) <- "integrated"
sce.all = sce.integrated

看看三个经典指标

VlnPlot(sce.all, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

FeatureScatter(sce.all, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

sce.all <- FindVariableFeatures(sce.all, nfeatures = 1500)
top10 <- head(VariableFeatures(sce.all), 10)
plot1 <- VariableFeaturePlot(sce.all)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot2

关于高变基因这一步呢，还是有点不得劲的。虽然原文里有这个图，作者也说是用vst方法挑出来的，但是seurat官方教程里整合数据后是不做normalize和找高变化基因的，直接就是scaleData了。而且整合后的数据不是count矩阵，用不了vst。细节没有太多描述，就这样继续做呗。

3. 降维聚类分群啦

跑PCA，选择多少个主成分用于后续分析

all.genes <- rownames(sce.all)
sce.all <- ScaleData(sce.all, features = all.genes)
sce.all[["integrated"]]@scale.data[30:34,1:3]
##        AAACCTGGTACAGACG-1_1 AAACGGGGTAGCGCTC-1_1 AAACGGGGTCCTCTTG-1_1
## CTSD              0.4832165            0.1260508           -0.8700229
## MT1X              1.0439902            0.7763393           -0.5553776
## FCGR3A            0.8605570           -0.5830022           -0.5830022
## TIMP1             0.3118005            2.1084136           -0.7909911
## KRT8             -0.2545172           -0.2545172           -0.2545172
sce.all <- RunPCA(sce.all, npcs = 30)
DimPlot(sce.all, reduction = "pca")

ElbowPlot(sce.all)

sce.all <- JackStraw(sce.all, num.replicate = 100)
sce.all <- ScoreJackStraw(sce.all, dims = 1:20)
JackStrawPlot(sce.all, dims = 1:20)

跑tsne

sce.all <- RunTSNE(sce.all,  dims = 1:15)
sce.all <- FindNeighbors(sce.all,dims = 1:15)
sce.all <- FindClusters(sce.all, resolution = 0.6) #分辨率
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1863
## Number of edges: 62779
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8709
## Number of communities: 11
## Elapsed time: 0 seconds
length(levels(Idents(sce.all)))
## [1] 11
table(sce.all@active.ident) 
## 
##   0   1   2   3   4   5   6   7   8   9  10 
## 384 377 207 184 169 149 146  94  72  50  31
p1 = DimPlot(sce.all,reduction = "tsne",group.by =  "orig.ident") 
p2 = DimPlot(sce.all,reduction = "tsne",label=T ) 
p1+p2

细胞不按照样本聚集，看起来糊在一起就对了。

SingleR注释

celldex的数据还要下载，碗束十分感人，还是不要下载了。我随便搜了一下：

如此敷衍的关键词也没问题哈哈哈哈。于是技能树的讲师通过搜索，获得了技能树本树的帮助。你点进去看看就知道ref_Hematopoietic.RData这个数据是怎么来的啦！

library(celldex)
library(SingleR)
#ref <- celldex::HumanPrimaryCellAtlasData()
ref <- get(load("ref_Hematopoietic.RData"))
library(BiocParallel)
pred.scRNA <- SingleR(test = sce.all@assays$integrated@data, 
                      ref = ref,
                      labels = ref$label.main, 
                      clusters = sce.all@active.ident)
pred.scRNA$pruned.labels
##  [1] "Monocytes"       "CD4+ T cells"    "CD8+ T cells"    "Dendritic cells"
##  [5] "Monocytes"       "Erythroid cells" "B cells"         "B cells"        
##  [9] "Monocytes"       "HSCs"            "HSCs"
plotScoreHeatmap(pred.scRNA, clusters=pred.scRNA@rownames, fontsize.row = 9,show_colnames = T)

new.cluster.ids <- pred.scRNA$pruned.labels
names(new.cluster.ids) <- levels(sce.all)
levels(sce.all)
##  [1] "0"  "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10"
sce.all <- RenameIdents(sce.all,new.cluster.ids)
levels(sce.all)
## [1] "Monocytes"       "CD4+ T cells"    "CD8+ T cells"    "Dendritic cells"
## [5] "Erythroid cells" "B cells"         "HSCs"
TSNEPlot(object = sce.all, pt.size = 0.5, label = TRUE)

搞掂~

我的本职工作是生信入门和数据挖掘线上直播课程讲师，如果想要系统学习搞定生信数据分析，可以发邮件或者来生信星球公号找到我，有缘分的话你迟早会来的~

其实我转移到了别的平台，简书很久才看一次，评论看到的时候基本已经过去很多天，另外平时工作繁忙，很少有功夫回复。如果是跟我的教程学习，有卡住的情况，迫切需要我帮助的问题，请参考生信星球公号的答疑公告(因简书不允许外链，所以无法放链接)，把问题描述的图文并茂、表达清楚自己的意思，发邮件到我的邮箱xjsun1221@163.com，野生博主，佛系回复，请礼貌提问，不要催我，也不要问我别人的代码错在了哪里哦。