scRNA-seq Preprocessing

전반적인 pre-processing 과정에 대해 정리.

(1) Create Seurat object
(2) Doublet Filter
(3) Quality Control Processing
(4) Data Integration
(5) Cell-Type Annotation

Step 1. Create Seurat object

우선은 가지고 있는 scRNA-seq data 데이터를 읽고 불러온 이후 seurat object 형태로 변환을 해줌. 아래 코드는 간단하게 survived, non-survived 환자군 데이터 한 개씩 불러와서 코드 진행하였음.

library(Seurat)
library(tidyverse)
library(scDblFinder)
library(harmony)
library(furrr)
library(future)
library(here)

here::i_am("main.R")

options(future.globals.maxSize = 4 * 1000 * 1024^2)
options(future.rng.onMisuse = "ignore")
plan("multicore", workers = 12)

set.seed(42)

survived_patients <- c(5)
survived_timepoints <- c(1)
survived_samples <- cross2(survived_patients, survived_timepoints) %>% map_chr(~ paste(.[1], .[2], sep = "_"))

non_survived_patients <- c(4)
non_survived_timepoints <- c(1)
non_survived_samples <- cross2(non_survived_patients, non_survived_timepoints) %>% map_chr(~ paste(.[1], .[2], sep = "_"))

create_seurat_object <- function(sample) {
      result <- Read10X_h5(here(paste0("data/processed/survived/", sample, "_output_filtered.h5")), use.names = TRUE, unique.features = TRUE)
      result_seurat_obj <- CreateSeuratObject(result, project = "Survived_sepsis")
      return (result_seurat_obj)
}

create_nonsurvived_seurat_object <- function(sample) {
      result <- Read10X_h5(here(paste0("/home/jhkim/2023_scRNA_sepsis/data/processed/non_survived/", sample, "_output_filtered.h5")), use.names = TRUE, unique.features = TRUE)
      result_seurat_obj <- CreateSeuratObject(result, project = "Non_survived_sepsis")
      return (result_seurat_obj)
}


s_single_objects <- survived_samples %>% future_map(create_seurat_object) %>% setNames(survived_samples)
ns_single_objects <- non_survived_samples %>% future_map(create_nonsurvived_seurat_object) %>% setNames(non_survived_samples)

survived <- s_single_objects$'5_1'
non_survived <- ns_single_objects$'4_1'

survived$status <- "Survived"
non_survived$status <- "Non-Survived"

Step 2. Doublet filter

이제 이렇게 seurat object가 만들어졌으면, 이후 과정으로 doublet filtering을 진행. Doublet filtering 과정은 scRNA-seq 데이터에서 doublets 또는 multiplets를 식별하고 제거하는 과정으로 sequencing 단계에서 하나의 droplet에 두 개 이상의 세포가 포함되어 발생한 경우를 말하는 것. 이제 분석을 진행하면 Doublet 같은 경우는 세포 클러스터링 과정에서 실제로 존재하지 않는 가짜 세포 유형을 만들어낸다거나 잘못된 상태가 존재하는 것 처럼 왜곡된 결과를 불러올 수 있으니 꼭 제거를 하고 진행.

이제 아래 코드를 확인해보면 우선 Seurat object를 SingleCellExperiement object로 변환을 해주고 scDblFinder package를 활용해 doublet 탐지를 진행하고, 해당 결과를 이제 metadata에 추가하여, singlet에 해당되는 데이터만 subset을 진행함.

check_doublet <- function(sample) {
    sample_sce <- as.SingleCellExperiment(sample)
    sample_sce <- scDblFinder(sample_sce)
    sample_class <- sample_sce$scDblFinder.class
    sample <- AddMetaData(object = sample, metadata = sample_class, col.name = "doublet_class")
    return (sample)
}

s_sample <- check_doublet(s_single_objects)
ns_sample <- check_doublet(ns_single_objects)

s_sample <- subset(s_sample, subset = doublet_class == "singlet")
ns_sample <- subset(ns_sample, subset = doublet_class == "singlet")

Step 3. Quality Control Processing

Doublet filter 이후 이제 한번더 데이터 품질 개선하고 노이즈 줄이기 위한 전처리를 진행.

• 우선 미토콘드리아 DNA 비율을 계산함. PercentageFeatureSet 함수를 통해서 MT-로 시작하는 유전자에 해당하는 feature들의 발현 비율을 계산함. 즉 mtDNA (미토 DNA) 비율을 계산하는 것인데, 높은 mtDNA 비율은 세포가 손상되었거나 죽어가고 있음을 나타내기 때문에 대체적으로 threshold 5-10 이상인 세포들에 대해서는 필터링을 진행.
• 또한 세포에서 발현되어 있는 유전자들의 수를 기반으로도 필터링을 진행. nFeature RNA는 각 세포에서 검출된 고유한 유전자의 수를 의미하는 것으로, 너무 적은 유전자가 검출된 세포는 데이터 품질이 낮을 가능성이 높고, 또 너무 많은 유전자가 검출된 세포는 이중 세포일 가능성이 있기도 하기에 VlnPlot을 통해서 outlier에 해당되는 부분들에 대해서도 filtering을 진행.

이 과정을 이제 샘플별로 진행을 하는지 아니면 샘플들 데이터를 묶어서 이제 풀링된 상태에서 VlnPlot을 통해 값을 확인해보고 threshold를 정하는데, 사람들 깃허브를 확인해보면 대체적으로 그냥 풀링된 상태로 qc를 진행함. 근데 사실 각각의 샘플별로 상태가 전부다 다를 것인데 정석대로 한다면 개별적으로 qc를 진행하는 것이 아닌가 싶긴 하지만 .. 하나씩 하는게 귀찮아서 그런건지 개별적으로 진행하는 사람은 거의 못봄.

s_sample[["percent.mt"]] <- PercentageFeatureSet(s_sample, pattern = "^MT-")
ns_sample[["percent.mt"]] <- PercentageFeatureSet(ns_sample, pattern = "^MT-")

VlnPlot(s_sample, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
VlnPlot(ns_sample, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

s_sample <- subset(s_sample, subset = nFeature_RNA > 200 & nFeature_RNA < 4000 & percent.mt < 10)
ns_sample <- subset(ns_sample, subset = nFeature_RNA > 200 & nFeature_RNA < 4000 & percent.mt < 10)

Step 4. Data Integration

이제 기본적인 qc 절차들을 진행했으면 각각의 샘플들에 대해서 합쳐주는 단계.
(우선은 seurat으로 integration 하는 과정, 근데 실제로 해보면 R은 너무 오래걸려서 나는 사실 여기까지 과정을 scanpy로 진행한 이후 RDS파일로 다시 변환해서 downstream 분석을 하는 것을 더 선호함)

그냥 단순하게 데이터 merge 하면 되지 않을가 생각할 수도 있지만, 결국 integration의 주 목적은 batch effects를 보정을 하기 위한 것. Batch effects란 실험 조건이나 실험자 등 실험의 다양한 세팅에서 발생하는 변동성을 의미하는데, 이러한 변동성 때문에 scRNA-seq 분석을 진행할 때 실제 생물학적 차이가 아니라 실험적 차이에서 비롯되 분석 결과에 오류를 초래할 수 있음. 그렇기 때문에 서로 다른 샘플들에 대해서 통합할 때 batch correction을 꼭 해줘야되고, 이거 안하면 데이터 신뢰성 너무 떨어지고 .. 나도 논문들 볼 때 이런 batch correction을 보정하지 않은 사람들 결과보면 신뢰성이 너무 떨어지니 꼭꼭 진행하기.

Batch correction을 위한 다양한 방법들이 있는데, "A benchmark of batch-effect correction methods for single-cell RNA sequencing data" 논문 확인해보면 여러 benchmark 기법들에 대한 비교를 하고 있음. 이전에 자문을 받았을 때 harmony 기법이 짱이라고 말씀하셨던 기억이 있고, 실제로 다른 기법들이랑 비교해보았을 때 harmony가 다른 BBKNN 같은 기법들보다 훨씬 보정이 잘 됬던 경험이 있음.

그래서 아래 코드를 보면
(1) 우선은 각 샘플 데이터를 merge 해준 다음 orig.ident로 각 데이터별로 나누어줌. 이렇게 샘플별로 데이터를 나누어준 이후 process 과정을 진행한 다음 다시 합쳐주면 배치 효과를 좀 줄일 수 있는 장점이 있음. (내 코드는 샘플 2개만 존재해서 split.by = "status"로 정의했지만 실제 분석을 하면 샘플 n수가 훨씬 많을 것이니 해당 부분에 sample 같은 metadata 컬럼을 삽입하면됨)
(2) 이후 이제 각 obj 들에 대해서 NormalizeData와 FindVariableFeatures 를 찾는 과정을 거치고,
(3) SelectIntegrationFeatures, FindIntegrationAnchors, IntegrateData 함수는 여러 데이터 세트를 통합하고 배치 효과를 보정하는 과정을 수행, 또 추가적으로 나의 경우는 "Harmony"를 통해서 추가적인 배치 보정을 적용해서, 각 배치의 데이터를 더 잘 조정한 이후 통합하였음.
(4) 이후 integrated된 데이터에 대해서 데이터 scaling -> PCA -> UMAP -> FindNeighbors -> Clustering 단계 진행 (자세한 사항은 Seurat 홈페이지 참고)

merge <- merge(s_sample, ns_sample)

obj_list <- SplitObject(merge, split.by = "status")

obj_list <- lapply(X = obj_list, FUN = function(x) {
  x <- NormalizeData(x)
  x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

features <- SelectIntegrationFeatures(object.list = obj_list)
immune_anchors <- FindIntegrationAnchors(object.list = obj_list, anchor.features = features)
immune_combined <- IntegrateData(anchorset = immune_anchors)

DefaultAssay(immune_combined) <- "integrated"
immune_combined <- ScaleData(immune_combined, verbose = FALSE)
immune_combined <- RunPCA(immune_combined)
immune_combined <- RunHarmony(immune_combined, group.by.vars = "orig.ident",reduction=‘pca’,reduction.save = “harmony”)
immune_combined <- RunUMAP(immune_combined, reduction = "harmony", dims = 1:10)
immune_combined <- FindNeighbors(immune_combined, reduction = "harmony", dims = 1:10)
immune_combined <- FindClusters(immune_combined, resolution = 1.0)

그래서 이제 최종 결과물을 확인해보면 아래와 같은 결과 확인 가능. Batch correction이 잘 보정된 것을 확인 가능하다.

Step 5. Cell-Type Annotation

이제 이렇게 데이터가 전부 합쳐졌다면 가장 귀찮고 손이 많이 가는 cell-type annotation를 진행해야됨.
사실 실제 manuscript를 위한 결과물이라면 각 cell-type들에 대한 cell marker들을 전부 찾아서 해당 유전자의 발현값을 확인해보고 해당 cluster에 메뉴얼로 어떤 세포에 속하는지에 대해 판별하는 과정을 거쳐야됨.. 근데 이게 실제로 해보면 정말로 오래 걸려서 이번에는 Azimuth package를 사용해서 cell-type annotation을 진행하였음.

Azimuth는 이제 cell-type annotation을 자동으로 진행해주는 패키지로, 이미 사전에 정의된 reference dataset과 비교해서, 각 세포에 가장 적절한 세포 유형을 할당해주는 패캐지임. 이제 우리 데이터셋의 세포 유형별로 유전자 발현 프로파일을 reference set과 비교를하고 두 데이터 간의 similarity를 측정 (유클리디안 distance, cosine similarity, correlation) 을 하게 됨.

그래서 단 한 줄로 cell-type annotation 과정을 끝낼 수가 있는 것이고, 실제로 돌려보면 predicted.celltype.score가 나와서 세포 유형에 대한 예측 값을 정량적으로 제공해줌. 값들을 확인해보면 0-1사이의 값들로 구성이되어 있고 이에 대한 threshold를 정해서 filtering을 진행해줘야됨. 근데 사실 이런 threshold를 내가 자의적으로 정하면 왜 이렇게 했어? 라는 질문들에 꼬리가 물려버릴 수 있으니 Nature Medicine에 published 된 method를 그대로 따라함. Nature Medicine "Single-cell multi-omics analysis of the immune response in COVID-19" method 부분을 확인해보면 이 논문도 동일하게 azimuth 를 사용하고 predicted.celltype.score가 0.5미만인 세포들에 대해서는 아예 필터링을 하고 분석을 진행하였음.

library(Azimuth)

df <- RunAzimuth(immune_combined, reference = "pbmcref", verbose = TRUE)
subset_df <- subset(df, subset = predicted.celltype.l2.score > 0.5)

그래서 최종 annotation 된 결과를 확인해보면 사실 어느정도 결과가 나쁘지는 않지만, 메뉴얼로 한 것보다는 결과가 잘나오지는 않는 것 같음. 특히 T cell 같은 경우는 명확하게 나누어주지를 못하는 것 같지만 모 시간을 아주 많이 줄여줬으니 threshold값을 좀더 올리거나 하면 좀더 명확히 구분되지 않을가 싶음.

그래서 여기까지가 Seurat 을 이용한 전체적인 pipeline에 대한 내용이였고 추후 글에는 scanpy로도 진행하는 방법과 여러 논문을 위한 figure 그리는 방법들에 대해서도 포스팅할 예정.
-끝-