前言
看到一篇单细胞数据挖掘的文章，题为：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)
fs = dir("GSE117570_RAW/");fs
## [1] "GSM3304007_P1_Tumor_processed_data.txt.gz"
## [2] "GSM3304011_P3_Tumor_processed_data.txt.gz"
## [3] "GSM3304013_P4_Tumor_processed_data.txt.gz"
fs2 = str_split(fs,"_",simplify = T)[,2];fs2
## [1] "P1" "P3" "P4"
原本是8个文件来着，这篇文章是只拿了其中3个。
rm(list = ls())
if(!file.exists("f.Rdata")){
fs = dir("GSE117570_RAW/")
f = lapply(paste0("GSE117570_RAW/",fs),function(x){
Matrix::Matrix(as.matrix(read.table(x,check.names = F)), sparse = T)
})
names(f) = fs2
save(f,file = "f.Rdata")
}
load("f.Rdata")
str(f,max.level = 1)
## List of 3
## $ P1:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## $ P3:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## $ P4:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
这个数据诡异，第一个样本和第三个样本里面有3个相同的barcode，需要处理掉。所以加上下面一段，正常数据里不加哦
length(intersect(colnames(f$P1),colnames(f$P4)))
## [1] 3
f$P4 = f$P4[,!(colnames(f$P4) %in% colnames(f$P1))]
3.创建Seurat对象
library(Seurat)
library(tidyverse)
library(patchwork)
obj = CreateSeuratObject(counts = f,min.cells = 3,min.features = 200)
names(obj@assays$RNA@layers)
## [1] "counts.P1" "counts.P3" "counts.P4"
CreateSeuratObject是可以一次容纳多个表达矩阵的，会存放在不同的layers
4.质控
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")
obj[["percent.rp"]] <- PercentageFeatureSet(obj, pattern = "^RP[SL]")
obj[["percent.hb"]] <- PercentageFeatureSet(obj, pattern = "^HB[^(P)]")
head(obj@meta.data, 3)
## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rp
## AAACCTGGTACAGACG-1 SeuratProject 4338 1224 2.512679 18.71830
## AAACGGGGTAGCGCTC-1 SeuratProject 11724 2456 2.021494 28.59092
## AAACGGGGTCCTCTTG-1 SeuratProject 3353 726 2.117507 55.26394
## percent.hb
## AAACCTGGTACAGACG-1 0
## AAACGGGGTAGCGCTC-1 0
## AAACGGGGTCCTCTTG-1 0
咔，发现orig.ident 是”SeuratObject”,而不是样本名，所以给它手动改一下了。
这两种写法都可以得到两个数据分别多少列，即多少个细胞。
c(ncol(f[[1]]),ncol(f[[2]]))
## [1] 1832 328
sapply(f, ncol)
## P1 P3 P4
## 1832 328 1420
obj@meta.data$orig.ident = rep(names(f),times = sapply(obj@assays$RNA@layers, ncol))
VlnPlot(obj,
features = c("nFeature_RNA",
"nCount_RNA",
"percent.mt",
"percent.rp",
"percent.hb"),
ncol = 3,pt.size = 0.1, group.by = "orig.ident")
obj = subset(obj,
percent.mt < 20 &
#nFeature_RNA < 4200 &
#nCount_RNA < 18000 &
percent.rp <50 #&
#percent.hb <1
)
ok接下来是
5.降维聚类分群那一套
obj <- NormalizeData(obj) %>%
FindVariableFeatures()%>%
ScaleData(features = rownames(.)) %>%
RunPCA(features = VariableFeatures(.)) %>%
IntegrateLayers(CCAIntegration)%>%
FindNeighbors(reduction = 'integrated.dr', dims = 1:15)%>%
FindClusters(resolution = 0.5)%>%
RunUMAP(reduction = "integrated.dr", dims = 1:15)%>%
RunTSNE(reduction = "integrated.dr", dims = 1:15)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 3319
## Number of edges: 119327
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8788
## Number of communities: 9
## Elapsed time: 0 seconds
UMAPPlot(obj)+TSNEPlot(obj)
obj = JoinLayers(obj)
obj
## An object of class Seurat
## 8013 features across 3319 samples within 1 assay
## Active assay: RNA (8013 features, 2000 variable features)
## 3 layers present: data, counts, scale.data
## 4 dimensional reductions calculated: pca, integrated.dr, umap, tsne
6.SingleR注释
library(celldex)
library(SingleR)
ls("package:celldex")
## [1] "BlueprintEncodeData" "DatabaseImmuneCellExpressionData"
## [3] "HumanPrimaryCellAtlasData" "ImmGenData"
## [5] "MonacoImmuneData" "MouseRNAseqData"
## [7] "NovershternHematopoieticData"
f = "../supp/single_ref/ref_BlueprintEncode.RData"
if(!file.exists(f)){
ref <- celldex::BlueprintEncodeData()
save(ref,file = f)
}
ref <- get(load(f))
library(BiocParallel)
scRNA = obj
test = scRNA@assays$RNA$data
pred.scRNA <- SingleR(test = test,
ref = ref,
labels = ref$label.main,
clusters = scRNA@active.ident)
pred.scRNA$pruned.labels
## [1] "Monocytes" "Epithelial cells" "CD8+ T-cells" "Epithelial cells"
## [5] "Macrophages" "Macrophages" "Mesangial cells" "B-cells"
## [9] "B-cells"
#查看注释准确性
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(scRNA)
levels(scRNA)
## [1] "0" "1" "2" "3" "4" "5" "6" "7" "8"
scRNA <- RenameIdents(scRNA,new.cluster.ids)
levels(scRNA)
## [1] "Monocytes" "Epithelial cells" "CD8+ T-cells" "Macrophages"
## [5] "Mesangial cells" "B-cells"
DimPlot(scRNA, reduction = "tsne",label = T,pt.size = 0.5) + NoLegend()
搞掂~
