单细胞多样本整合之harmony(seurat v5)

用户11414625

发布于 2024-12-20 15:24:36

47400

代码可运行

文章被收录于专栏：生信星球520生信星球520

运行总次数：0

代码可运行

0.背景知识

多样本整合这回事，其实有很多中方法，seurat5可以用一个参数支持5种整合算法。 Anchor-based CCA integration (method=CCAIntegration) Anchor-based RPCA integration (method=RPCAIntegration) Harmony (method=HarmonyIntegration) FastMNN (method= FastMNNIntegration) scVI (method=scVIIntegration)

不要选择困难症，我们一般就是用harmony。单纯harmony一种算法的话，之前V4版本的代码也仍然可以使用，这不是有新用法了再跟着跑一下嘛。

官方文档使用的示例数据和包下载起来一堆坑，让我们搞简单点，用GEO的数据GSE183904的两个样本。

为啥是两个样本呢，因为样本数量多了你电脑带不起来。。。

单细胞数据有很多种格式，csv是其中一种。

dir("GSE183904_RAW/")

## [1] "GSM5573467_sample2.csv.gz" "GSM5573472_sample7.csv.gz"

1.批量读取多个样本的数据

rm(list = ls())
if(!file.exists("f.Rdata")){
  fs = dir("GSE183904_RAW/")
  f = lapply(paste0("GSE183904_RAW/",fs),function(x){
    Matrix::Matrix(as.matrix(read.csv(x,row.names = 1)), sparse = T)
  })
  fs = stringr::str_split_i(fs,"_",1)
  names(f) = fs
  save(f,file = "f.Rdata")
}
load("f.Rdata")
str(f,max.level = 1)

## List of 2
##  $ GSM5573467:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##  $ GSM5573472:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots

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.GSM5573467" "counts.GSM5573472"

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
## AAACCTGAGGAGTAGA_9 SeuratProject       4237          691   2.666981   5.947604
## AAACCTGAGGGTGTTG_9 SeuratProject       6156         2448   6.270305  11.306043
## AAACCTGAGGTTCCTA_9 SeuratProject      16410         2905   3.266301   7.160268
##                     percent.hb
## AAACCTGAGGAGTAGA_9 0.000000000
## AAACCTGAGGGTGTTG_9 0.000000000
## AAACCTGAGGTTCCTA_9 0.006093845

咔，发现orig.ident 是”SeuratObject”,而不是样本名，所以给它手动改一下了。

这两种写法都可以得到两个数据分别多少列，即多少个细胞。

c(ncol(f[[1]]),ncol(f[[2]]))

## [1] 2769 2097

sapply(f, ncol)

## GSM5573467 GSM5573472 
##       2769       2097

obj@meta.data$orig.ident = rep(names(f),times = sapply(bj@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 < 5 &
#              nFeature_RNA < 4200 &
#              nCount_RNA < 18000 &
#              percent.rp <30 &
#              percent.hb <1
# )

这个数据已经被过滤过了，就不用再过滤了。

ok接下来是

5.降维聚类分群那一套

obj <- NormalizeData(obj) %>%
  FindVariableFeatures()%>%
  ScaleData(features = rownames(.)) %>%  
  RunPCA(features = VariableFeatures(.))  %>%
  IntegrateLayers(HarmonyIntegration)%>%
  FindNeighbors(reduction = 'harmony', dims = 1:15)%>%
  FindClusters(resolution = 0.5)%>%
  RunUMAP(reduction = "harmony", dims = 1:15)%>%
  RunTSNE(reduction = "harmony", dims = 1:15)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 4866
## Number of edges: 171641
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9137
## Number of communities: 13
## Elapsed time: 0 seconds

UMAPPlot(obj)+TSNEPlot(obj)

obj = JoinLayers(obj)
obj

## An object of class Seurat 
## 20010 features across 4866 samples within 1 assay 
## Active assay: RNA (20010 features, 2000 variable features)
##  3 layers present: data, counts, scale.data
##  4 dimensional reductions calculated: pca, harmony, 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] "CD8+ T-cells"      "B-cells"           "CD8+ T-cells"     
##  [4] "Epithelial cells"  "Monocytes"         "Epithelial cells" 
##  [7] "B-cells"           "Fibroblasts"       "HSC"              
## [10] "Epithelial cells"  "Endothelial cells" NA                 
## [13] "Epithelial 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"  "9"  "10" "11" "12"

scRNA <- RenameIdents(scRNA,new.cluster.ids)
levels(scRNA)

## [1] "CD8+ T-cells"      "B-cells"           "Epithelial cells" 
## [4] "Monocytes"         "Fibroblasts"       "HSC"              
## [7] "Endothelial cells"

DimPlot(scRNA, reduction = "tsne",label = T,pt.size = 0.5) + NoLegend()