Seurat으로 scRNA-seq 데이터 다루기

사전 정보¶

scRNA seq과 10xGenomics¶

scRNA-seq는 single-cell RNA sequencing의 줄임말로, 하나의 세포에서 mRNA를 측정하는 방법입니다. 이 기술은 기존 bulk RNA-seq 방법과는 달리 하나의 세포에서 RNA를 추출하여 분석합니다. 이를 통해, 개별 세포의 유전자 발현 패턴, 전사체 감지, 변형과 발현의 상호작용 등을 이해할 수 있습니다.

10xGenomics는 scRNA-seq 분석에서 매우 인기있는 플랫폼으로 droplet-based 방법을 사용합니다. droplet-based 방법은 cell barcoding 및 unique molecular identifier(UMI)를 사용하여 RNA-seq 라이브러리를 생성하는 공정으로 사실상 현재 scRNA seq분야에 표준으로 사용됩니다.

사전 작업¶

10xGenomics 기술을 사용해 scRNA-seq을 진행하면 결과로 fastq 파일을 얻게 됩니다. fastq 파일을 Seurat 패키지에 바로 사용할 수는 없고, 10xGenomics에서 제공하는 Cell Ranger를 사용해 맵핑을 진행해야 합니다.

시퀀싱 데이터 준비¶

실험을 통해 다음과 같은 fastq 파일을 가지고 있다고 간주합니다. 여기에서는 학습 목적으로 아주 작은 데이터셋으로 진행합니다만 실제 데이터는 훨씬 큽니다.

.(pbmc_1k_v3_fastqs)
├── pbmc_1k_v3_S1_L001_I1_001.fastq.gz
├── pbmc_1k_v3_S1_L001_R1_001.fastq.gz
├── pbmc_1k_v3_S1_L001_R2_001.fastq.gz
├── pbmc_1k_v3_S1_L002_I1_001.fastq.gz
├── pbmc_1k_v3_S1_L002_R1_001.fastq.gz
└── pbmc_1k_v3_S1_L002_R2_001.fastq.gz

cell ranger 설치¶

여기에서는 설치 방법은 생략하고 공식 홈페이지 링크를 참조하시기 바랍니다.

cell ranger 실행¶

cellranger count --id=run_count_1kpbmcs \
   --fastqs=./pbmc_1k_v3_fastqs \
   --sample=pbmc_1k_v3 \
   --transcriptome=./refdata-gex-GRCh38-2020-A
   --nosecondary

위의 명령어를 통해 cell ranger를 실행할 수 있습니다. --id 는 생성되는 결과의 폴더명이며, --fastqs는 fastq 파일이 있는 폴더의 위치, --sample은 metadata에 들어가는 샘플 정보, --transcriptome는 참조 전사체의 위치 입니다.

저의 경우는 10x Genomics 사이트에서 다음의 명령어로 다운로드 받았습니다. 참고로 참조 전사체는 실험에 사용된 시료에 따라 다르게 사용해야 합니다.

wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz
tar -zxvf refdata-gex-GRCh38-2020-A.tar.gz

cell ranger 결과¶

cell ranger가 문제없이 작동했다면 out 폴더에 다음과 같은 파일과 폴더가 생겨납니다.

.(out)
├── analysis
├── cloupe.cloupe
├── filtered_feature_bc_matrix
├── filtered_feature_bc_matrix.h5
├── metrics_summary.csv
├── molecule_info.h5
├── possorted_genome_bam.bam
├── possorted_genome_bam.bam.bai
├── raw_feature_bc_matrix
├── raw_feature_bc_matrix.h5
└── web_summary.html

이중에 Seurat 패키지가 필요로 하는 것은 filtered_feature_bc_matrix 폴더 입니다. 폴더안에는 다음과 같은 파일이 들어있습니다.

.(filtered_feature_bc_matrix)
├── barcodes.tsv.gz
├── features.tsv.gz
└── matrix.mtx.gz

이걸로 모든 사전 준비가 완료되었습니다.

사용할 패키기 불러오기¶

Seurat¶

Seurat은 R 프로그래밍 언어로 작성된 scRNA-seq 데이터 분석을 위한 유명한 패키지 중 하나입니다. Seurat은 높은 차원의 scRNA-seq 데이터에서 유전자 발현 패턴을 탐색하고 이를 이용하여 세포 및 클러스터의 식별과 분석, 특성 제시, 시각화 등 다양한 분석 작업을 수행할 수 있습니다. Seurat은 다양한 데이터 전처리 및 정규화 기능과 함께 차원 축소, 클러스터링, 시각화, 서브셋 생성, 유전자 발현 분석, 세포간 상호작용 분석 등 다양한 분석 도구를 제공합니다. Seurat은 현재까지 업데이트와 기능 추가가 활발하게 이루어지고 있으며, scRNA-seq 분석에 필수적인 유틸리티 패키지 중 하나입니다.

scDblFinder¶

scDblFinder는 단일 세포 RNA 시퀀싱 데이터에서 더블렛(두 개의 세포가 동시에 캡처되어 하나의 세포로 보이는 것) 현상을 탐지하고 제거하기 위한 R 패키지입니다. 이 패키지는 UMI(Unique Molecular Identifier)를 기반으로하여 두 개 이상의 세포에서 동시에 탐지된 UMIs를 찾아서 더블렛으로 추정하고, 더블렛으로 추정된 셀을 제거합니다. 이를 통해 scRNA-seq 데이터의 정확도와 해석력을 높일 수 있습니다. 또한, scDblFinder는 Seurat 및 SingleCellExperiment 형식의 데이터를 지원하며, 다양한 분석 옵션을 제공하여 사용자가 데이터에 맞게 조정할 수 있습니다.

tidyverse¶

tidyverse는 데이터 분석에 필요한 필수 R 패키지들의 모음으로 데이터 처리, 시각화, 모델링, 프로그래밍 등의 다양한 작업을 수행하는 데 사용됩니다. 주요 패키지로는 ggplot2, dplyr, tidyr, readr, purrr, tibble, stringr, forcats 등이 있습니다.

패키지 버전 확인¶

사용한 Seurat 패키지의 버전¶

사용한 scDblFinder 패키지의 버전¶

Seurat 개체 만들기¶

먼저 데이터를 읽어오는 것으로 시작합니다. Read10X() 함수를 사용하면 간단하게 UMI count 행렬을 불러 올수 있습니다. UMI count 행렬은 앞서 cell ranger를 통해 만들어 진 것으로 scRNA-seq 실험에서 얻어진 유전자 발현 데이터를 행렬 형태로 나타낸 것으로, 행은 각각의 세포를 나타내고, 열은 유전자를 나타냅니다.

그런 다음 CreateSeuratObject()함수를 사용해 Seurat 개체를 생성합니다. 그러면 UMI count 행렬은 seurat_obj[["RNA"]]@counts에 저장되게 됩니다.

count 행렬은 어떻게 생겼을까?¶

Seurat 개체의 메타데이터는 어디에 저장될까?¶

seurat_obj@meta.data 혹은 seurat_obj[[]]을 통해 메타데이터를 확인할 수 있습니다.

Pre-processing¶

서열 데이터 품질 관리¶

scRNA seq 데이터의 분석의 신뢰성을 얻기 위해서 데이터의 품질 관리는 필수이빈다. Seurat을 사용하면 품질관리(QC) 지표를 쉽게 탐색하고 사용자 정의 기준에 따라 세포를 필터링할 수 있습니다. 일반적으로 사용되는 QC 기준은 다음 세가지 입니다.

• 각 세포에서 검출된 고유 유전자의 수.
  • 품질이 낮은 세포는 종종 유전자가 매우 적습니다.
  • 이중 또는 다중의 세포가 들어간 drop-let에는 비정상적으로 유전자 수가 높습니다.
• 각 세포에서 검출된 총 서열의 수(고유 유전자의 수와 밀접한 상관관계가 있음)
• 미토콘드리아 게놈에 매핑되는 서열의 비율
  • 품질이 낮거나 죽어가는 세포에는 미토콘드리아 유전자가 많이 발견됩니다.

QC 지표 시각화하기¶

바이올린 플랏¶

Scatter 플랏¶

QC 및 추가 분석을 위한 세포 선택하기¶

위의 결과를 토대로 nFeature_RNA가 200개에서 6000개 사이이고 percent.mt가 20 이하인 세포들만 선택합니다.

아직 모든 QC 과정이 끝난 것은 아닙니다. PCA를 진행하고나서 scDblFinder 패키지를 사용해 추가적인 더블렛 데이터를 제거하겠습니다.

데이터 정규화하기¶

QC를 통해 일부 데이터를 제거한 다음 단계는 데이터를 정규화하는 것입니다. 여기서는 각 세포의 발현 값을 전체 발현으로 나누고 스케일 계수(기본적으로 10,000)를 곱한 다음 로그 변환하는 LogNormalize방법을 사용합니다. 이렇게 정규화된 값은 seurat_obj[["RNA"]]@data에 저장됩니다.

Elbow 플랏 그리기¶

Seurat에서 clustering을 수행하기 전에는 몇 개의 차원(dimension)을 사용할지 결정해야 합니다. 차원의 수는 PCA와 같은 차원 축소 기법을 사용하여 줄여진 차원의 수를 의미합니다. 그리고 이 차원의 수는 클러스터링 알고리즘에 사용됩니다.

그러나 차원의 수가 너무 적거나 많으면 적절한 클러스터링이 어려울 수 있습니다. 차원이 적을 경우 정보 손실이 크게 발생하고, 차원이 많을 경우에는 불필요한 차원의 포함으로 인해 과적합(overfitting)이 발생할 가능성이 있습니다.

따라서 적절한 차원의 수를 선택하기 위해 elbow plot을 사용합니다. elbow plot은 차원의 수를 x축으로, 해당 차원의 데이터를 잘 설명하는 정도(예: variance)를 y축으로 나타냅니다. 이 때, 차원의 수를 늘리면 y축 값은 점점 증가하게 됩니다. 그러나 어느 지점 이후로는 y값이 더 이상 크게 증가하지 않고 평평해지는 지점이 나타나는데, 이 지점이 elbow point입니다. 이 지점 이후로는 차원을 늘려도 데이터를 잘 설명하지 못하므로, elbow point를 기준으로 적절한 차원의 수를 선택합니다. 이를 통해 데이터의 차원을 축소할 때, 적절한 차원의 수를 선택하여 과적합을 방지하고 필요한 정보만을 추출할 수 있습니다.

Clustering¶

elbow 플랏에서 elbow point가 대략적으로 10임을 알 수 있습니다. 따라서 clustering의 dims 변수를 1:10 으로 지정합니다. 그런다음 UMAP을 그려보겠습니다.

UMAP 그리기¶

UMAP은 Uniform Manifold Approximation and Projection의 약자로, scRNA-seq 데이터를 시각화하기 위한 비선형 차원 축소 방법 중 하나입니다. t-SNE와 유사한 기능을 가지고 있지만, 대규모 데이터셋에서 더욱 빠르고 정확한 임베딩을 제공합니다.

UMAP은 데이터의 국부적인 구조를 보존하는데 초점을 둡니다. 즉, 비슷한 특성을 가진 데이터들이 서로 가깝게 묶이고, 서로 다른 특성을 가진 데이터들은 더 멀리 배치되도록 임베딩을 생성합니다. 이를 통해, scRNA-seq 데이터의 복잡한 구조를 파악하고 시각화할 수 있습니다.

위의 UMAP 플랏을 통해 총 10개의 cluster로 나누어 졌음을 알 수 있습니다. 이제 추가적인 QC를 진행해보겠습니다.

scDblFinder를 사용해 추가 QC 하기¶

scDblFinder는 클러스터 정보를 기반으로 인공적으로 생성된 더블렛을 찾아냅니다. 아래 코드를 통해 더블렛을 찾고 UMAP 플랏에 표시해보겠습니다.

scDblFinder 결과에서 Threshold 값이 0.519 라는 것과 총 30(2.6%)개의 더블렛이 계산되었습니다. 다시 subset()함수를 사용해 scDblFinder 값이 0.519 이하인 세포만 고르는 QC 과정을 진행하겠습니다.

더블렛 제거하기¶

Cell type identification¶

각각의 클러스터들이 어떤 세포인지 알아내는 작업을 합니다.Seurat의 FindAllMarkers()함수를 사용하면 나머지 모든 세포와 비교해 클러스터에 대한 유전자 마커를 찾을 수 있습니다.

Marker gene 찾기¶

FindAllMarkers() 결과는 데이터프레임입니다. avg_log2FC는 다른 클러스터와 비교해 발현량이 얼마나 차이나는지를 의미합니다. 해당 열을 가지고 각 클러스터당 상위 5개의 유전자 마커를 추려서 heatmap 을 그려봅니다.

Marker gene heatmap¶

위의 heatmap은 멋져보이기는 하지만 그렇게 큰 정보를 제공해주지는 않습니다. 일반적으로 클러스터가 어떤 세포인지 알아내는 작업이 scRNA seq 분석방법에서 가장 중요한 부분이며 여러가지 접근법이 있지만 가장 덜 자동화된 부분이기도 합니다.

가장 일반적인 방법은 FindAllMarkers() 함수를 이용해 얻은 각 클러스터의 유전자 마커와 참조 데이터 세트를 비교 하는 것입니다. 참조 데이터 세트란, 이미 잘 정의된 세포 유형들의 scRNA-seq 데이터로 현재 분석 중인 데이터의 클러스터들을 참조 데이터 세트와 비교하여 유사한 패턴을 가진 세포 유형을 찾아내는 것입니다. 이 작업은 시간이 아주 많이 필요하며 작업자에 따라 다른 결과가 나옵니다. 그래서 이번에는 ChatGPT를 사용하는 방법으로 해보겠습니다.

ChatGPT 사용하기¶

ChatGPT에 대한 자세한 설명은 생략합니다.

프롬프터에 클러스터의 유전자 마커를 쉽게 입력하기 위해 다음의 코드를 사용합니다.

위 코드의 출력을 가지고 ChatGPT에 입력했더니 다음과 같이 말했습니다.

주어진 유전자 목록만으로는 각각의 클러스터의 정확한 세포 유형을 알 수 없습니다. 유전자는 여러 종류의 세포에서 발현될 수 있기 때문입니다. 유전자의 발현 패턴을 분석하거나 추가적인 실험 결과와 함께 고려해야만 각 클러스터의 세포 유형을 추론할 수 있습니다.

그러나 일반적으로 각 클러스터에 대한 예상 세포 유형은 다음과 같습니다:

$cluster0: 이 클러스터는 PBMC에서 neutrophils을 나타내는 것으로 추정됩니다.
$cluster1: 이 클러스터는 PBMC에서 CD4 T cells을 나타내는 것으로 추정됩니다.
$cluster2: 이 클러스터는 PBMC에서 B cells을 나타내는 것으로 추정됩니다.
$cluster3: 이 클러스터는 PBMC에서 B cells을 나타내는 것으로 추정됩니다.
$cluster4: 이 클러스터는 PBMC에서 monocytes을 나타내는 것으로 추정됩니다.
$cluster5: 이 클러스터는 PBMC에서 NK cells을 나타내는 것으로 추정됩니다.
$cluster6: 이 클러스터는 PBMC에서 B cells을 나타내는 것으로 추정됩니다.
$cluster7: 이 클러스터는 PBMC에서 CD8 T cells을 나타내는 것으로 추정됩니다.
$cluster8: 이 클러스터는 PBMC에서 NK cells을 나타내는 것으로 추정됩니다.
$cluster9: 이 클러스터는 PBMC에서 platelets을 나타내는 것으로 추정됩니다.

위의 결과를 통해 대부분의 클러스터 타입을 알 수 있지만, B cell이 3개 NK cell이 2개로 중복됩니다. 좀 더 세분화하기 위해 B cell의 유전자 마커를 물어봅니다.

B cell은 다양한 sub type으로 분류될 수 있습니다. 일부 sub type들과 그에 따른 대표 마커들은 다음과 같습니다.

- Naive B cell
    CD19, CD20, CD27-, IgD+, IgM+
- Memory B cell
    CD19, CD20, CD27+, IgD-, IgM+/-
- Plasma cell
    CD19, CD20-, CD38++, CD138+, MUM1+, IRF4+, XBP1+
- Germinal center B cell
    CD19, CD20, CD77+, CD10+, BCL6+, AID+
- Regulatory B cell (Breg)
    CD19, CD20, CD24++, CD38++, CD25+, IL-10+
- B1 B cell
    CD19, CD20, CD27+, IgM+, CD43+
- Marginal zone B cell
    CD19, CD20, CD27-, CD43-, IgM+, IgD+, CD23-, CD21+

각 sub type마다 대표적인 마커들이 있지만, 이것들만으로 모든 sub type을 완벽하게 분류하는 것은 어렵습니다. 따라서, sub type 분류는 보통 여러 마커를 조합하여 결정됩니다.

이번에는 NK cell 세분화를 위한 유전자 마커를 물어봅니다.

NK 세포는 크게 두 가지 subset으로 나눌 수 있습니다: CD56bright NK 세포와 CD56dim NK 세포입니다. CD56bright NK 세포는 주로 면역 조절에 관여하며, cytokine을 분비하는 역할을 합니다. 반면 CD56dim NK 세포는 주로 세포 독성에 관여하며, 인체 내에서 암세포와 바이러스에 대한 첫 번째 방어선 역할을 합니다.

CD56bright NK 세포를 구별하는 유전자 마커로는 CD56(bright) (NCAM1)과 CD16 (FCGR3A)이 있습니다. 반면 CD56dim NK 세포를 구별하는 유전자 마커로는 CD56(dim) (NCAM1)과 KIR (Killer-cell immunoglobulin-like receptors)가 있습니다. KIR은 NK 세포의 분화 상태, 활성 상태, 대상 세포와의 상호 작용 등을 조절합니다.

위의 결과와 pbmc1k_marker.csv 파일을 비교하면서 다음과 같이 클러스터의 이름을 지정합니다.

ChatGPT를 사용한 방법은 정확하지 않고 힌트정도만 알 수 있습니다. 실제 실험결과라면 SingleR 패키지 혹은 문헌 검색을 통해 직접 찾아보시기 바랍니다.

그리고 UMAP 플랏을 그려서 결과를 확인합니다.

마치며¶

클러스터에 포함된 세포 수를 확인해보고 RDS 파일을 저장하면서 마무리 하겠습니다.

클러스터당 세포 수 확인¶

클러스터당 세포 비율 확인¶

RDS 파일로 저장하기¶

RDS 파일을 저장해두면 추후 분석에 위의 과정을 반복할 필요가 없습니다.

지금 까지 살펴본 내용을 요약해보면 다음과 같습니다.

• Pre-processing: 데이터 전처리를 통해 불필요한 변수 제거, 정규화 등을 수행합니다.
• Dimensionality reduction: 차원 축소 기법을 사용해 데이터의 주요 구조를 파악합니다.
• Clustering: 유사한 특성을 가진 데이터들을 그룹화합니다.
• Cell type identification: 각 클러스터에 대해 유전자 발현 패턴 등을 비교하여 cell type을 추론합니다.

이후 진행되는 scRNA-seq Downstream analysis는 여기서 얻은 결과를 기반으로 합니다. 따라서 여기서의 결과가 부정확하거나 신뢰성이 떨어지면 추가 분석에서 얻은 결과 또한 의미가 없습니다. 그러므로 분석에서 사용된 데이터의 품질, 분석 방법의 적절성, 도구의 성능 등을 철저히 검토하고 확실한 기준에 따라 분석을 수행하세요. 또한, 추후에 데이터나 분석 방법이 변경되는 경우 이전의 결과와의 비교를 통해 신뢰성을 유지할 수 있도록 관리하는 것도 잊지 말아야 합니다.