【R code】Seurat CCAを用いてreference dataを元に自分のdataをmappingする方法

タイトルにあるように、自分のdataがreference dataに対してどこのclusterに属しているかを知りたい、という状況になることがあると思います。

しかし、reference dataと自分のdataにbatch effectが存在する場合やそもそもreference data自体にもbatch effectがありそうなdataである場合にはbatch effect除去を行う必要が生じます。

Seuratではその様な方法を提示してくれています；https://satijalab.org/seurat/articles/integration_mapping.html

Seuratのデモンストレーションではない場合、私の場合のコードを紹介します。

Reference dataの読み込み、clusterの抽出

# reference dataの入力
A.data, B.data, C.dataをCreateSeuratObjectで読み込み
your.data <- merge(A.data, y= c(B.data, C.data), add.cell.ids = c("A.data","B.data","C.data"), project = "original reference")

## your.dataを解析して、必要なclusterのみを抽出する場合 
#### cluster 0-13, 15, 17を抽出したい場合
Cell.idents <- c(rep("OK",14), "NG", "OK", "NG", "OK")
names(Cell.idents) <- levels(your.data)  
your_sub.data <- your.data
your_sub.data <- RenameIdents(your_sub.data, Cell.idents) #clusterのannotationを変更
your_sub.data@meta.data$cell_type <- your_sub.data@active.ident
your_sub.data <- subset(x=your_sub.data, subset=cell_type=="OK")

# CCAを行うためにSeuratObjectをlist化する
## batch差がA.data, B.data, C.dataのサンプル間に存在しているため、orig.identでdataを分ける。
your_cca.list <- SplitObject(your_sub.data, split.by ="orig.ident")

Reference dataに対してCCAでbatch effect除去を行い、解析方法として問題ないかを確認する

# ここからはseuratのdemonstrationとほぼ同様に行う
## 各batchに対してNormalizationとFindVariableFeatures
for (i in 1:length(your_cca.list)) {
    your_cca.list[[i]] <- NormalizeData(your_cca.list[[i]], verbose = FALSE)
    your_cca.list[[i]] <- FindVariableFeatures(your_cca.list[[i]], selection.method = "vst", nfeatures = 2500, verbose = FALSE)
}
## CCA施行
your_cca.anchors <- FindIntegrationAnchors(object.list = your_cca.list, dims = 1:30)
your_cca.integrated <- IntegrateData(anchorset = your_cca.anchors, dims = 1:30)

## 使用するassayをintegrated(=CCAしたdata）を使ってくださいと指定
DefaultAssay(your_cca.integrated) <- "integrated"

## Run the standard workflow for visualization and clustering
your_cca.integrated <- ScaleData(your_cca.integrated, verbose = FALSE)
your_cca.integrated <- RunPCA(your_cca.integrated, npcs = 30, verbose = FALSE)
your_cca.integrated <- RunUMAP(your_cca.integrated, reduction = "pca", dims = 1:30, verbose = FALSE)
your_cca.integrated <- FindNeighbors(your_cca.integrated, reduction = "pca", dims = 1:30)
your_cca.integrated <- FindClusters(your_cca.integrated, resolution = 0.5)

これでrefarence dataのCCAによるintegrationは完了です。

調べたいdata(= query data)を読み込み

# query data(mappingしたい方)のdataの入力
query.data <- merge(e.data, y= c(f.data, g.data), add.cell.ids = c("e.data","f.data","g.data"), project = "query_data")

# seuratobjectをlist化して、Normalization, FindVariableFeatureを行う
query.data <- SplitObject(query.data, split.by = "orig.ident")
for (i in 1:length(query.data)) {
    query.data[[i]] <- NormalizeData(query.data[[i]], verbose = FALSE)
    query.data[[i]] <- FindVariableFeatures(query.data[[i]], selection.method = "vst", nfeatures = 2500, verbose = FALSE)
}

# query側はlist化されてはいけないので、このままのdataをSeuratObjectとしてmergeさせる
merge.list <- c("f.data","g.data")
merge.query <- query.data[["e.data"]]
for (i in merge.list){
    merge.query <- merge(merge.query,y= query.data[[i]],merge.data=TRUE)
} ## merge.dataをTRUEにすることでNormalizationされたdataもそのまま引き継ぐ

reference data上にmerge.queryのdataをmappingする

your_cca.anchors <- FindTransferAnchors(reference = your_cca.integrated, query = merge.query, dims = 1:30, reference.reduction = "pca")

predictions <- TransferData(anchorset = your_cca.anchors, refdata = your_cca.integrated@active.ident, dims = 1:30)
merge.query <- AddMetaData(merge.query, metadata = predictions)

your_cca.integrated <- RunUMAP(your_cca.integrated, dims = 1:30, reduction = "pca", return.model = TRUE)

merge.query <- MapQuery(anchorset = your_cca.anchors, reference = your_cca.integrated, query = merge.query, refdata = list(celltype = "seurat_clusters"), reference.reduction = "pca", reduction.model = "umap")

このようにすればdata形成は完了です。あとはご自身の選択するplot方法で描写してください。

まとめ

最近はR Studioで解析を毎日やっているので、Blog の内容もそっちがメインとなってしまいますね。

今回のcodeはほぼSeuratのHPで確認できることと同じですが、seurat demoでは既にまとまっているデータなので、少し参考にならない部分もありました。私が少しデータの読み込み、形式の変換等で戸惑ったので、自分でデータを読み込ませる必要がある人に参考になれば幸いです。