시작하기 앞서¶

다음과 같은 몇가지 가정을 하고 시작하겠습니다.

• 가정1. 10x Genomics Single Cell Solution을 통해 얻은 서열 데이터를 사용
• 가정2. 리눅스(Ubuntu 22.04.3 LTS x86_64) OS를 사용하며 기본적인 shell 명령어를 앎
• 가정3. 패키지 매니저인 Anaconda에 대한 사용법을 알고 있음

터미널을 열고 원하는 폴더로 이동합니다. 저는 /home/fkt/Downloads에서 시작하도록 하겠습니다.

$pwd
/home/fkt/Downloads

가상환경 만들기¶

다음과 같이 사용가상환경을 만들어 사용합니다. 먼저 TutorialEnvironment.yml 파일을 다음 명령어로 다운받습니다.

해당 파일에는 RNA velocity 분석을 위한 velocyto 도구의 dependancy 정보가 들어 있습니다.

wget https://cf.10xgenomics.com/supp/cell-exp/neutrophils/TutorialEnvironment.yml

다음 명령어로 가상환경을 만들어 줍니다.

conda env create --file TutorialEnvironment.yml

만들어진 가상환경을 활성화 시킵니다.

conda activate tutorial

다운로드 받은 yml파일은 삭제해줍니다.

rm TutorialEnvironment.yml

Cell ranger 설치¶

10x Genomics Single Cell Solution을 통해 single cell RNA sequncing을 했다면 10x Genomics에서 Cell ranger라는 alignment 도구를 제공하고 있습니다. 속도는 좀 느려도 사용 용이성은 좋아서 저는 Cellranger를 추천합니다. 자세한 것은 공식사이트를 참고하시기 바랍니다.

시스템 요구사항¶

Cellranger는 리눅스 시스템에서 작동하며 다음의 최소 사양을 요구합니다:

• 8-core Intel or AMD processor (16 cores recommended).
• 64GB RAM (128GB recommended).
• 1TB free disk space.
• 64-bit CentOS/RedHat 7.0 or Ubuntu 14.04 - See the 10x Genomics OS Support page for details.

터미널에 다음 명령어를 입력해 cell ranger 7.2버전을 다운로드 합니다.

현재 기준으로 7.2 버전이 최신입니다.

wget -O cellranger-7.2.0.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-7.2.0.tar.gz?Expires=1695129121&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA&Signature=ZipqR8Pg4YvVDQ3MvAVuuiPwEOC5c39~Wj0WTAxfoJW6xtrxqIFgIDySFSFnsWcwDpmAovJGrHXU24Y9Cptt88OJSPdEupyFRXoGKJvVzRtDJChmuMbSpVCy-2N-QnMKwxNtd8Yt8Mdp2Vcq4wxx1hVC0Yx54c7U9o~RFXIVsIp48thKR6JnKhJmCAC5U4dFLa86~NcB4s5Ic4HATrQP2KyWexyYZWgmBEw13mlKYtlVRUil0zseoq0-CZyGmE8oB0iDBSUBAyIqo~XMVjv~lkMz4cRcyCbQKBRDr~U36FM2KnE3rhv-Rlp4KD-uXCReRnBsY6N6t-HxpS2YDpN3mQ__"

그리고 나서 다음 명령어로 압축을 풀어줍니다.

tar -xvf cellranger-7.2.0.tar.gz

cellranger-7.2.0 폴더로 이동하기

cd cellranger-7.2.0

정상적으로 되었다면 아래와 같은 폴더 구조를 볼 수 있습니다.

.
├── bin
├── builtwith.json
├── cellranger -> bin/cellranger
├── external
├── lib
├── LICENSE
├── mro
├── probe_sets -> external/tenx_feature_references/targeted_panels
├── sourceme.bash
├── sourceme.csh
└── THIRD-PARTY-LICENSES.cellranger.txt

5 directories, 6 files

이제 ./cellranger명령어를 입력하면 아래와 같이 cell ranger 실행이 가능합니다.

./cellranger
cellranger cellranger-7.2.0

Process 10x Genomics Gene Expression, Feature Barcode, and Immune Profiling data

Usage: cellranger <COMMAND>

Commands:
  count           Count gene expression and/or feature barcode reads from a single sample and GEM well
(중략)
  help            Print this message or the help of the given subcommand(s)

Options:
  -h, --help     Print help
  -V, --version  Print version

FASTQ 파일과 레퍼런스 데이터 준비¶

10x Genomics사에서는 학습을 위한 데이터를 제공하고 있습니다. 이미 많은 예제에서 Human 1k PBMC를 다루고 있으므로 저는 Mouse 1k neuron으로 진행해보겠습니다. Human 1k PBMC은 직접해보시기 바랍니다.

Human PBMC 1k 데이터 다운로드¶

터미널에 다음과 같이 입력합니다.

mkdir data && cd data
#  Human PBMC 1k = 5.17GB
wget https://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_1k_v3/pbmc_1k_v3_fastqs.tar
tar -xvf pbmc_1k_v3_fastqs.tar
# Human reference (GRCh38) = 11GB
wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz
tar -zxvf refdata-gex-GRCh38-2020-A.tar.gz

Mouse 1k neuron 데이터 다운로드¶

터미널에 다음과 같이 입력합니다.

# Mouse neron 1k data = 7GB
wget http://cf.10xgenomics.com/samples/cell-exp/3.0.0/neuron_1k_v3/neuron_1k_v3_fastqs.tar
tar -xvf neuron_1k_v3_fastqs.tar
# Mouse reference = 9.7GB
wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-mm10-2020-A.tar.gz
tar -zxvf refdata-gex-mm10-2020-A.tar.gz

다운로드된 파일 확인¶

정상적으로 진행되었다면 아래와 같은 폴더와 52개의 파일을 확인할 수 있습니다.

$ tree .
.
├── neuron_1k_v3_fastqs
│   ├── neuron_1k_v3_S1_L001_I1_001.fastq.gz
│   ├── neuron_1k_v3_S1_L001_R1_001.fastq.gz
│   ├── neuron_1k_v3_S1_L001_R2_001.fastq.gz
│   ├── neuron_1k_v3_S1_L002_I1_001.fastq.gz
│   ├── neuron_1k_v3_S1_L002_R1_001.fastq.gz
│   └── neuron_1k_v3_S1_L002_R2_001.fastq.gz
├── 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
├── refdata-gex-GRCh38-2020-A
│   ├── fasta
│   ├── genes
│   │   └── genes.gtf
│   ├── pickle
│   ├── reference.json
│   └── star
└── refdata-gex-mm10-2020-A
    ├── fasta
    ├── genes
    │   └── genes.gtf
    ├── pickle
    ├── reference.json
    └── star

12 directories, 52 files

Cell ranger 사용하기¶

데이터가 준비되었으니 Cell ranger를 사용해 count matrix를 만들어 보겠습니다. 먼저 상위 디렉토리인 cellranger-7.2.0으로 이동합니다.

cd ..

다음과 같이 cellranger count 명령어를 사용하면 count matrix데이터를 얻을 수 있습니다.

./cellranger count --id=run_count_mNeuron1k \
--fastqs=data/neuron_1k_v3_fastqs  \
--sample=neuron_1k_v3  \
--transcriptome=data/refdata-gex-mm10-2020-A

--id는 생성되는 폴더명이며 --fastqs는 다운로드한 fastq파일의 위치입니다. --sample은 fastq 파일명의 접두사(prefix)를 의미하고 --transcriptome은 레퍼런스데이터의 위치입니다.

시간이 좀 지나면 다음과 같은 출력을 볼 수 있습니다.

Running preflight checks (please wait)...
Checking sample info...
Checking FASTQ folder...
Checking reference...
(후략)

정상적으로 작업이 완료된다면 아래처럼 outs이라는 폴더가 생성됩니다. 여기에는 count matrix뿐만 아니라 여러가지 분석 결과도 포함되어 있어 복잡해보입니다. 그러나 저희가 필요 한 것은 filtered_feature_bc_matrix 폴더 안에 다들어 있습니다.

outs
├── 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

RNA velocity 분석을 위한 작업¶

RNA velocity이란? 스플라이싱된 RNA와 아직 스플라이싱되지 않은 RNA를 구분해 세포의 상태를 예측하는 기법으로 scRNA 데이터에서 클러스터간에 분화 과정에 대한 단서를 얻을 수 있습니다. 예를 들어 아래 그림에서 보면 cluster4에서 cluster2로의 방향성을 확인 할 수 있습니다.

안타깝게도 Cell ranger를 통해 얻은 count matrix로는 RNA velocity 분석을 할 수 없습니다. 추가적으로 velocyto이라는 도구를 사용해 unspliced , spliced 데이터가 포함된 loom파일을 만들어야 합니다. 아래 명령어를 사용하면 간단하게 얻을 수 있습니다.

velocyto run10x run_count_mNeuron1k/ data/refdata-gex-mm10-2020-A/genes/genes.gtf

실행이 완료되는 시간은 생각보다 오래 걸립니다.

2023-09-19 10:01:22,756 - DEBUG - Using logic: Default
2023-09-19 10:01:22,757 - DEBUG - Example of barcode: AAACGAATCAAAGCCT and cell_id: run_count_mNeuron1k:AAACGAATCAAAGCCT-1
(중략)
2023-09-19 10:27:23,985 - DEBUG - 1991314 reads were skipped because no apropiate cell or umi barcode was found
2023-09-19 10:27:23,985 - DEBUG - Counting done!
2023-09-19 10:27:23,988 - DEBUG - Collecting row attributes
2023-09-19 10:27:24,014 - DEBUG - Generating data table
2023-09-19 10:27:24,056 - DEBUG - Writing loom file
2023-09-19 10:27:25,286 - DEBUG - Terminated Succesfully!

완료가 되면 run_count_mNeuron1k폴더안에 velocyto라는 새로운 폴더가 생기고 run_count_mNeuron1k.loom 파일이 생성됩니다.

Count matrix 데이터 읽기¶

이제 힘든 데이터 준비 과정이 모두 끝났습니다. 데이터 분석을 시작해보죠. scRNA-seq 분야는 아직 초기이고 새로운 도구들이 계속 쏟아지고 있습니다. 나중은 어떻게 될지 모르겠지만 지금 시점에서는 R에 익숙하다면 Seurat, Python에 익숙하다면 Scanpy를 추천합니다.

참고로 속도적인 측면과 ML을 접목하기에는 Scanpy가 더 강점을 가지고 있습니다.

Scanpy 사용하기¶

다음 코드를 통해 count matrix를 scanpy로 불러 올 수 있습니다.

scVelo를 사용해 RNA velocity분석을 한다면 아래의 추가적인 코드가 필요합니다.

scanpy는 *.h5ad 포멧을 사용해 count matrix를 저장합니다. 아래 코드를 사용하면 단일 파일로 모든 것을 저장할 수 있습니다.

Seurat 사용하기¶

R에 좀 더 익숙하시다면 Seurat을 사용하세요. 저는 파이썬 문법을 더 좋아하지만 솔직히 분석 도구의 완성도 측면에서는 Seurat이 더 나은 것 같습니다. 아래 코드는 jupyter notebook에서 R코드를 사용하기 위해 필요한 것으로 rStudio를 사용하신다면 %%R이 포함된 셀의 코드만 사용하시면 됩니다.

필요한 Seurat 라이브러리를 불러옵니다.

아래의 코드를 사용하면 count matrix를 Seurat object로 만들 수 있습니다.

Seurat object를 파일로 저장하는 방법에도 여러가지가 있지만 R에서 많이 사용되는 *.rds를 사용하는 코드는 아래와 같습니다.

마치며¶

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