单细胞多样本整合之CCA(Seuratv5)

前言

看到一篇单细胞数据挖掘的文章，题为：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()

搞掂~