서로 다른 scRNA-seq data integration하기!_Public Data [R]

바이오인포메틱스

서로 다른 scRNA-seq data integration하기!_Public Data [R]_1

미이리 2022. 10. 20. 19:28

SMALL

# Visualize the Louvain clustering and the batches on the UMAP. 
# Remember, the clustering is stored in @meta.data in column seurat_clusters and the technology is
# stored in the column tech. Remember you can also use DimPlot
DimPlot(marge.data, reduction = "umap", group.by = "seurat_clusters")
DimPlot(marge.data, reduction = "umap", group.by = "each")
DimPlot(marge.data, reduction = "umap", group.by = "paper")
DimPlot(marge.data, reduction = "umap", group.by = "time")

single cell RNA seq 분석하는 것을 공부하면서 여러가지 이것도 가능한가? 저것도 가능한가? 하면서 궁금한 점이 정말 많았습니다. 그러다가 single cell RNA seq data를 이용하여 meta-analysis한 논문을 보게 되었는데요. 그래서 나도 저런 분석을 할 수 있게 되면, 나에게 정말 큰 도움이 될 것이라 생각해 이렇게 저렇게 분석을 해봤습니다. 어떤 방법이 맞는 방법인지는 해보면서 알아보려합니다.

1. Seurat tutorial에 나온 "Introduction to scRNA-seq integration"을 이용해서 분석해보기.

library(Seurat)
library(dplyr)
library(tidyverse)
library(tidyseurat)
library(cowplot)
library(patchwork)
library(DoubletFinder)

우선 필요한 package들을 시작을 해주고 실험 방법이 비슷한 두 논문에서 나온 public data를 다운받아 불러들였습니다.

# Import matrix data
Mock.data <- Read10X(data.dir = "C:/Rstudy/meta/Fiege/S1")
h24.data <- Read10X(data.dir = "C:/Rstudy/meta/Fiege/S2")
h48.data <- Read10X(data.dir = "C:/Rstudy/meta/Fiege/S5")

mock.R.data <- Read10X(data.dir = "C:/Rstudy/meta/Ravindra/mock")
h24.R.data <- Read10X(data.dir = "C:/Rstudy/meta/Ravindra/1dpi")
h48.R.data <- Read10X(data.dir = "C:/Rstudy/meta/Ravindra/2dpi")

각각 서로 다른 실험실에서 실험이 진행이 된 논문입니다. 하지만 실험에 쓰인 세포와 바이러스를 감염시킨 시간이 동일한 data를 이용해서 integration해보려 합니다. 우선 각각 SeuratObject를 만들어줍니다.

#createSeuratObject
Mock.F <- CreateSeuratObject(Mock.data,min.cells = 3, min.features = 200)
h24.F <- CreateSeuratObject(h24.data,min.cells = 3, min.features = 200)
h48.F <- CreateSeuratObject(h48.data,min.cells = 3, min.features = 200)

mock.R <- CreateSeuratObject(mock.R.data,min.cells = 3, min.features = 200)
h24.R <- CreateSeuratObject(h24.R.data,min.cells = 3, min.features = 200)
h48.R <- CreateSeuratObject(h48.R.data,min.cells = 3, min.features = 200)

그리고 하나의 Object로 만들어주기 전에 normalization을 진행합니다.

Mock.F <- subset(Mock.F, subset = nFeature_RNA > 1000 & nFeature_RNA < 12000 & percent.mt < 30)
h24.F <- subset(h24.F, subset = nFeature_RNA > 1000 & nFeature_RNA < 12000 & percent.mt < 30)
h24.F <- subset(h24.F, subset = nFeature_RNA > 1000 & nFeature_RNA < 12000 & percent.mt < 30)

mock.R <- subset(mock.R, subset = nFeature_RNA > 1000 & nFeature_RNA < 10000 & percent.mt < 30)
h24.R <- h24.R(Mock.F, subset = nFeature_RNA > 1000 & nFeature_RNA < 10000 & percent.mt < 30)
h48.R <- subset(h48.R, subset = nFeature_RNA > 1000 & nFeature_RNA < 10000 & percent.mt < 30)

Mock.F  <- NormalizeData(Mock.F , normalization.method = "LogNormalize", scale.factor = 10000)
h24.F <- NormalizeData(h24.F, normalization.method = "LogNormalize", scale.factor = 10000)
h48.F <- NormalizeData(h24.F, normalization.method = "LogNormalize", scale.factor = 10000)

mock.R <- NormalizeData(mock.R, normalization.method = "LogNormalize", scale.factor = 10000)
h24.R <- NormalizeData(h24.R, normalization.method = "LogNormalize", scale.factor = 10000)
h48.R <- NormalizeData(h48.R, normalization.method = "LogNormalize", scale.factor = 10000)

Mock.F <- FindVariableFeatures(Mock.F, selection.method = "vst", nfeatures = 2000)
h24.F <- FindVariableFeatures(h24.F, selection.method = "vst", nfeatures = 2000)
h48.F <- FindVariableFeatures(h24.F, selection.method = "vst", nfeatures = 2000)

mock.R <- FindVariableFeatures(mock.R, selection.method = "vst", nfeatures = 2000)
h24.R <- FindVariableFeatures(h24.R, selection.method = "vst", nfeatures = 2000)
h48.R <- FindVariableFeatures(h48.R, selection.method = "vst", nfeatures = 2000)

Mock.F <- ScaleData(Mock.F)
h24.F <- ScaleData(h24.F)
h48.F <- ScaleData(h48.F)

mock.R <- ScaleData(mock.R)
h24.R <- ScaleData(h24.R)
h48.R <- ScaleData(h48.R)

#paper
Mock.F[["paper"]] <- "Fiege"
h24.F[["paper"]] <- "Fiege"
h48.F[["paper"]] <- "Fiege"

mock.R[["paper"]] <- "Ravindra"
h24.R[["paper"]] <- "Ravindra"
h48.R[["paper"]] <- "Ravindra"

# time-point
Mock.F[["time"]] <- "mock"
h24.F[["time"]] <- "24h"
h48.F[["time"]] <- "48h"

mock.R[["time"]] <- "mock"
h24.R[["time"]] <- "24h"
h48.R[["time"]] <- "48h"

#each
Idents(Mock.F) <- "each"
Idents(h24.F) <- "each"
Idents(h48.F) <- "each"
Idents(mock.R) <- "each"
Idents(h24.R) <- "each"
Idents(h48.R) <- "each"

나중에 시간별, 저자별, 그리고 각각의 object로 나눌 수 있게 object에 또 다른 meta data를 만들어줍니다.

cluster를 하기 위해 merge를 해줍니다. 처음에는 batch correction 없이 어떻게 cluster되나 먼저 확인할 것 입니다.

add.cell.ids <- c("Mock.F", "h24.F", "h48.F", "mock.R", "h24.R", "h48.R")

marge.data <- merge(x = Mock.F, y = list(h24.F, h48.F, mock.R,h24.R,h48.R  ), add.cell.ids = add.cell.ids, merge.data = FALSE)

Idents(marge.data) <- "time"  # use identity based on sample identity
Idents(marge.data) <- "paper"
Idents(marge.data) <- "each"

# Look at how the number of genes per cell varies across the different technologies.
VlnPlot(marge.data, "nFeature_RNA", group.by = "time")
VlnPlot(marge.data, "nFeature_RNA", group.by = "paper")
VlnPlot(marge.data, "nFeature_RNA", group.by = "orig.ident")

그러면 이렇게 원하는 데로 나뉘어서 볼 수 있는 것을 알 수 있습니다.

하지만 이렇게 merge를 한 다음에 역시 normalize와 scaling을 해야 합니다. 하지만 저는 이번에 scale the variable genes만 할 것입니다.

# The merged data must be normalized and scaled (but you only need to scale the variable genes). 
# Let us also find the variable genes again this time using all the data.
marge.data <- NormalizeData(marge.data, normalization.method = "LogNormalize", scale.factor = 10000)
var.genes <- SelectIntegrationFeatures(SplitObject(marge.data, split.by = "each"), verbose = TRUE, selection.method = "vst")

VariableFeatures(marge.data) <- var.genes
marge.data <- ScaleData(marge.data, features = VariableFeatures(marge.data))

# Do PCA on data including only the variable genes.
marge.data <- RunPCA(marge.data, features = VariableFeatures(marge.data), npcs = 40, ndims.print = 1:5, nfeatures.print = 5)

# Color the PC biplot by the scRNA-seq technology. Hint: use DimPlot
# Which technologies look similar to one another?
DimPlot(marge.data, reduction = "pca", dims = c(1,2), group.by = "each")

SMALL

FindNeighbors와 FindClusters를 통해서 merge한 데이터에 각 cell이 가까운 것들끼리 모이게 clustering해줍니다.

# Cluster the cells using the first twenty principal components.
marge.data <- FindNeighbors(marge.data, reduction = "pca", dims = 1:20, k.param = 20)
marge.data  <- FindClusters(marge.data , resolution = 0.8, algorithm = 1, random.seed = 100)

# Create a UMAP visualization. 
marge.data <- RunUMAP(marge.data, dims = 1:20, reduction = "pca", n.neighbors = 15, min.dist = 0.5, spread = 1, metric = "euclidean", seed.use = 1)

완벽하게 batch effect가 나타나는 것을 볼 수 있다. 그러고 이제 batch correction을 진행할 것이다. 나는 seurat v4에 있는 방식을 사용할 것이다.

integration하는 것은 내일 이어서 정리하겠습니다. 그럼 좋은 하루 보내세요!!

LIST