前言
看到一篇单细胞数据挖掘的文章，题为：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，野生博主，佛系回复，请礼貌提问，不要催我，也不要问我别人的代码错在了哪里哦。
