前言
看到一篇單細(xì)胞數(shù)據(jù)挖掘的文章屈呕，題為：Establishment of a Prognostic Model of Lung Adenocarcinoma Based on Tumor Heterogeneity
遂打算拿里面的數(shù)據(jù)跑一跑乐导，這個(gè)數(shù)據(jù)可以在GSE117570的補(bǔ)充文件里直接下載到喻犁。
1.批量讀取數(shù)據(jù)
雖然不是標(biāo)準(zhǔn)10X的三個(gè)文件碌补，但也可以搞，直接讀取為數(shù)據(jù)框，轉(zhuǎn)換為矩陣奈揍，自行創(chuàng)建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個(gè)文件來著赋续，這篇文章是只拿了其中3個(gè)男翰。
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
這個(gè)數(shù)據(jù)詭異，第一個(gè)樣本和第三個(gè)樣本里面有3個(gè)相同的barcode纽乱，需要處理掉蛾绎。所以加上下面一段，正常數(shù)據(jù)里不加哦
length(intersect(colnames(f$P1),colnames(f$P4)))
## [1] 3
f$P4 = f$P4[,!(colnames(f$P4) %in% colnames(f$P1))]
3.創(chuàng)建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是可以一次容納多個(gè)表達(dá)矩陣的鸦列，會(huì)存放在不同的layers
4.質(zhì)控
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
咔租冠，發(fā)現(xiàn)orig.ident 是”SeuratObject”,而不是樣本名，所以給它手動(dòng)改一下了薯嗤。
這兩種寫法都可以得到兩個(gè)數(shù)據(jù)分別多少列顽爹，即多少個(gè)細(xì)胞。
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"
#查看注釋準(zhǔn)確性
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()
搞掂~
