안녕하세요, 이번 글에서는 저번 글에 이어 단일세포 전사체 분석 기술의 프로세스를 정리해보도록 하겠습니다.
[44일차] 단일세포 전사체 분석 기술 개요 및 프로세스 정리 01 :: sequencing library 준비 및 구축, illum
안녕하세요, 이번엔 단일세포 전사체 분석 기술(single cell RNA-seq analysis)의 전반적인 개요와 고려사항들을 정리해보는 시간을 가져보고자 합니다. 공부 겸 정리하는 건지라 오류가 있을 수 있다는
tkmstudy.tistory.com
저번엔 데이터 분석을 위한 준비 절차였다면, 이번엔 본격적으로 데이터를 전처리하고 분석하는 절차로 보면 되겠습니다. 참고자료는 뒷부분에 적어두었으며, 딥리서치와 논문 'The Workflow for Computational Analysis of Single-cell RNA-sequencing Data'을 주로 참고하였습니다. 본 내용은 공부로 적은 내용인지라 오류가 있을 수 있다는 점 양해 부탁드립니다.
전 글에서 조직 샘플로부터 시퀀싱 데이터를 얻는 과정을 살펴보았습니다. 이제 10X Genomics의 Cell Ranger 파이프라인을 사용해 FASTQ 파일 양식으로 된 시퀀싱 데이터를 처리하여 '세포 별 유전자 발현 행렬(cell by genes matrics)를 얻을 차례입니다. 참고로 여기서 FASTQ 파일 양식은 DNA 서열 데이터를 저장하는 표준 양식 중 하나로, 시퀀싱 결과로 얻은 Read(시퀀싱된 DNA 단편)에 대한 서열과 품질 점수를 포함합니다. 아래와 같이 말이죠.
구체적으로, FASTQ 파일은 '@'를 포함한 서열 ID, 서열, +(설명), 품질 정보가 담긴 ASVII format의 '특수 문자' 이런 식으로 이어지는데요. 여기서 읽지못한 정보는 'N'으로 등장한다고 합니다. 보통 품질정보에 대문자가 많으면 퀄리티가 좋은 파일이라고 했던 거 같습니다.
Cell ranger 파이프라인이 FASTQ 데이터를 처리하기 위해선 분석에 맞는 참조 데이터(Reference)가 필요합니다. 우리는 de novo sequencing을 하려는 것이 아닌, resequencing을 하려고 하니까요. 우리는 과거 여러 사람들의 노력 덕에 인간 게놈 지도를 그릴 수 있었죠. 이제 우리는 사람에게서 얻은 특정 샘플 속 RNA의 염기서열을 보고(실제론 cDNA 서열을 보고), 그 염기서열이 어떤 유전자에 대응하는 염기서열인지 알 수 있습니다. 10X Genomics에선 미리 구축한 다양한 종(ex. 인간, 쥐 등)의 참조 데이터를 제공해줍니다. 따라서 우리는 이 참조 서열 데이터를 활용해 각각의 염기서열 '리드'가 어떤 유전자들에 대응되는지 확인합니다. 다시한번 정리하자면, Cell Ranger는 참조 서열을 기반으로 해당 전사체 샘플의 서열들이 어떤 유전자로부터 전사된 것인지 파악하고(Alignment), 각 유전자에 할당된 전사체들의 개수를 센 뒤, 세포별 유전자 발현 행렬(cell-by genes matrics)과 QC 리포트를 제공해줍니다.
| cell-by genes matrics (간단히 만든 예시) |
gene A | gene B |
| cell A | 30 | 10 |
| cell B | 0 | 30 |
그렇게 저희는 10X Genomics 시퀀싱 데이터 결과물로 세 가지 파일을 제공받을 수 있는데요, 각각 'barcodes.tsv', 'features.tsv', 'matrix.mtx' 파일입니다. 먼저, 바코드 파일은 앞서 비드에 붙어있던 Cell barcode 서열이 담긴 파일입니다. 이 바코드는 각 유전자가 어떤 세포에서 유래되었는지 알려주죠. 다음으로 feature 파일의 경우엔 각 유전자의 ID와 그에 맞는 유전자 이름을 제시해주는 파일인데요. 아래 이미지를 참고하시면 되겠습니다. 마지막 매트릭스 파일은 세포별 유전자들의 발현량 정보를 담은 파일입니다.
이 세 가지 파일들을 R 환경에 불러오면 CreateSeuratObject 함수로 Seurat 객체로 만들 수 있습니다. 여기서 Seurat은 대표적인 R 환경의 scRNA-seq 데이터 분석 툴이라고 보면 되겠습니다.
Tools for Single Cell Genomics
A toolkit for quality control, analysis, and exploration of single cell RNA sequencing data. Seurat aims to enable users to identify and interpret sources of heterogeneity from single cell transcriptomic measurements, and to integrate diverse types of sing
satijalab.org
# 예시코드, n1/엔 각각에 바코드 파일, 피쳐 파일, 메트릭프 파일 세 개의 파일이 존재
n1 <- Read10X(data.dir = "n1/")
n1 <- CreateSeuratObject(counts = n1, project = "normal", min.cells = 3, min.features = 200)
파이썬의 경우에는 대표적인 scRNA-seq 데이터 분석 툴로 scanpy가 있습니다. 저는 처음에 말했듯 Seurat을 기준으로 설명해보도록 하겠습니다. 당연한 말이지만, 좋은 데이터를 사용해야 올바른 분석 결과를 얻을 수 있습니다. 따라서 데이터 분석에 있어 데이터 전처리는 필수라고 볼 수 있겠습니다. 결국 저희가 seurat 객체에 대해 가장 먼저해주어야 할 일은 바로 'Quality Control(QC)'입니다.
1) Quality Control (QC)
QC 단계는 세포별 유전자 발현 행렬(raw count matris)에 대해 저품질 세포를 필터링함으로써 노이즈를 제거하는 단계입니다. 이 단계에서 저품질 세포는 물론, 드롭릿(droplet)에 세포가 들어가지 않았거나 두 개 이상의 세포가 들어간 경우도 필터링해줍니다. 드롭릿에 대한 기억을 잠시 되짚어봅시다. 앞서 말했듯, 본 시퀀싱 데이터는 10X Genomics의 droplet-based microfulidic 장치를 통해 시퀀싱 라이브러리를 준비해 시퀀싱을 돌림으로써 얻어낸 결과입니다. 여기서 각 droplet에는 하나의 세포와 하나의 비드만 들어가야 합니다(물론 역전사효소와 Master Mix도 들어가지만요). 그러나, 기술적 한계로 droplet에 세포가 들어가지 않거나 두 개 이상의 드롭릿이 들어가게 되는 경우가 존재합니다. 이러한 적절치 못한 드롭릿은 드롭릿의 양적 특성을 확인하여 QC 단계에서 제거해줍니다.
그렇다면 QC 단계에서 제거해야 할 '저품질 세포가 들어간 드롭릿'과 '적절치 못한 수의 세포가 들어간 드롭릿'은 어떻게 확인할까요? 첫번째 방법은 세포마다 감지된 고유의 유전자 개수를 확인하는 것입니다. 빈 드롭릿엔 세포가 없으니 드롭릿 내 총 유전자의 개수가 매우 낮겠죠(존재할 수는 있는데 뒤에서 설명하겠습니다). 이외에도 손상이 되거나 죽은 품질이 낮은 세포에서는 매우 적은 유전자가 검출되므로 QC 단계에서 세포마다 감지된 유전자 개수(nFeature_RNA)의 하한선을 적용합니다. 이때 하한선의 기준은 세포별 유전자 수를 시각화하고, 샘플의 특성 및 연구 목적을 고려해 연구자가 재량껏 설정해야 합니다. 그래도 보통 세포 당 200~300개 유전자를 넘는 세포로 필터링하는 것 같습니다. 주의해야 할 점은 "'유전자 수가 적은 세포'가 '죽거나 손상된 저품질 세포'가 아닌 단지 '생물학적으로 유전자 수가 적은 건강한 세포(건강한 암세포일 수 있음)'일 수 있다는 점"입니다. 따라서 연구자의 재량껏 하한선을 설정해야하는 것이죠.
앞서 제가 빈 드롭릿에 세포가 없어도 유전자가 검출될 수 있다고 했는데요. 이는 손상되거나 죽은 세포에서 나온 RNA(‘ambient RNA’)가 자유롭게 떠돌다가 빈 드롭릿으로 들어갈 수 있기 때문입니다. 이렇게 떠돌아다니는 RNA는 빈 드롭릿에서 유전자 수를 인위적으로 늘려, 실제로 유전자 수가 적은 세포를 가진 드롭릿과 혼동을 줄 수 있습니다. 따라서 연구자들은 빈드롭릿에서 ambient RNA가 일반적으로 나타나는 수준을 계산한 뒤, 빈 드롭릿과 유전자 수가 적은 세포를 포함하는 드롭릿을 구분할 수 있는 알고리즘을 개발했습니다. 그것은 바로 'EmptyDrops'입니다. EmptyDrops는 빈 드롭릿에서 보이는 RNA의 배경 수준(background level)을 추정하고, 그 배경 수준에서 크게 벗어나는 세포를 선별합니다 1).
사실 세포를 포함하고 있는 드롭릿에서도 ambient RNA가 포함될 수 있습니다. ambient RNA는 말 그대로 애매모호한 떠돌이 RNA이기 때문입니다. 그렇기에 QC 절차에서는 정상적인 드롭릿에서 ambient RNA만 필터링하기 위한 절차도 진행합니다. 대표적인 ambient RNA 제거 툴로는 'SoupX'가 있는데요. 이는 빈드롭릿이나 낮은 유전자 카운트를 보이는 드롭릿들을 이용해, ambient RNA의 발현 패턴('soup')을 추정합니다. 이때, 두 가지 가정을 전제로 합니다. 첫째, 빈 드롭릿에서 측정된 RNA들은 전적으로 외부 오염에 기인했다고 가정합니다. 둘째, 드롭릿마다 비드를 통해 검출된 총 RNA 수는 실제 세포에서 발현된 RNA와 ambient RNA 오염의 합으로 가정합니다. 이렇게 가정하고 SoupX는 각각의 세포 내에서 기대되지 않는 특정 유전자(즉, 특정 세포에서 발현되어야 할 필요가 없는 유전자들)의 발현을 확인하며, 각 세포에서 얼마나 많은 ambient RNA가 섞여 들어갔는지를 추정합니다. 이후 추정된 오염 비율을 이용해 각 세포의 관측된 유전자 카운트에서 ambient RNA의 기여를 빼주죠. 그렇게 실제 세포에서 유래한 순수한 전사체 데이터만을 복원하도록 합니다. 이것이 바로 SoupX가 ambient RNA에 의한 잘못된 신호를 제거하는 원리입니다. 챗GPT 내용을 참조했는데 자세하고 정확한 내용은 아래 링크의 설명과 논문을 참고하시길 바랍니다.
GitHub - constantAmateur/SoupX: R package to quantify and remove cell free mRNAs from droplet based scRNA-seq data
R package to quantify and remove cell free mRNAs from droplet based scRNA-seq data - constantAmateur/SoupX
github.com
지금까지 저희는 QC 단계 중 세포마다 감지된 유전자 개수가 적은 저품질 세포와, 빈드롭릿, 그리고 비지 않은 드롭릿에서의 ambient RNA에 의한 잘못된 신호를 제거하는 방법들을 알아보았습니다. 아직 끝이 아닙니다. 우리는 세포마다 감지된 유전자 개수가 비정상적으로 많은 세포도 제거해주어야 합니다. 이렇듯 유전자 개수가 많은 경우는 두 개 이상의 세포가 하나의 드롭릿에 포획되었을 수 있는데요. '두 개'나 '두 개 이상의 세포'가 캡처된 드롭릿을 각각 'Doublet'과 'Multiplet'이라고 표현합니다. 물론 유전자 개수가 많은 경우도 사실 생물학적으로 세포 자체가 부피가 커서 유전자를 많이 발현하는 세포일 수 있으니 이 또한 시각화와 도메인 지식을 기반으로 연구자 재량껏 주의하여 상한선을 두어 필터링을 진행해야 할 것입니다.
말은 연구자 재량껏이라고 했지만, 사실 과학에서는 누구나 동의할만한 '객관성'과 '연구 재현성'이 최대한 보장되어야 하죠. 따라서 최근에는 보다 신중한 Doublet 결정 및 필터링을 위한 툴들도 개발이 되었습니다. 예로, scDblfinder의 경우 두 개의 무작위로 선택된 드롭릿의 유전자 발현 데이터를 결합함으로써 '인공 더블릿(ariticial doublet)'을 생성합니다. 이후 인공 더블릿과의 유사도를 기준으로 각 드롭릿(세포)의 doublet score를 계산합니다. 이때, KNN(k-neareset neigborhood) 알고리즘을 활용한다고 하는데요. 각 드롭릿의 주변 이웃에 얼마나 많은 인공 더블릿이 포함되어 있는지를 보고, 각 드롭릿이 인공 더블릿과 얼마나 유사한지에 대한 점수를 매기는 것입니다. 이외에도 더블릿을 필터링을 하기 위한 여러 툴들이 개발되고 있는 상황입니다. 일반적으로 10x 플랫폼에서는 목표 셀 수의 ~5% 정도가 더블릿으로 추정되며, 필터링으로 이를 최대한 제거한다고 하네요.
사실 우리는 세포마다의 유전자 개수를 넘어, 각 유전자의 발현량까지 알고 있습니다. 즉, 세포마다의 총 RNA 수를 알고 있죠. 용어가 헷갈리는데 간단히 말해, RNA 시퀀싱 결과 특정 세포에서 '유전자A'와 '유전자B'에 해당하는 RNA가 각각 2개 3개 씩 캡쳐되었다면, 총 유전자 개수는 2개(A와 B), 총 RNA 수는 5개(2+3)로 계산합니다. 총 RNA 수는 총 UMI의 수라고도 볼 수 있는데요. 1편에서 정리했듯 시퀀싱 전에 각 RNA마다 고유의 UMI를 갖도록 했으니까요. 그렇게 우리는 총 UMI의 수를 통해 총 RNA 수를 알 수 있습니다. 이때 총 RNA 수도 유전자 개수에서처럼 너무 적은 RNA 수를 가진 세포는 죽거나 손상된 세포 혹은 빈 드롭릿으로 간주할 수 있고, 너무 많은 RNA 수를 가진 세포는 더블릿 혹은 멀티플렛으로 간주할 수 있습니다. 그러면서도 너무 적은 혹은 너무 많은 RNA 수를 가지는 세포가 단지 생물학적으로 다른 세포에 비해 부피가 작거나 큰 건강한 세포일 가능성도 있죠. 사실 유전자 개수가 특이하게 적거나 많은 세포는 총 RNA 수 또한 특이하게 적거나 많을 가능성이 높습니다. 그렇지만, 가끔 논문을 보면 총 RNA 수가 지나치게 많은 세포가 여럿 존재하는 경우 유전자 개수를 기준으로 필터링하는 것을 넘어 총 유전자 수로도 필터링을 진행하더군요. 또한, Seurat의 FeatureScatter 등을 이용해 두 지표(세포마다의 유전자 개수, 총 RNA 수)의 이상치를 함께 시각화하여 필터링 기준을 정할 수 있습니다(아래 오른쪽 이미지 참조).
한가지 더 중요한 필터링 기준이 있는데요. 세포에서의 미토콘드리아 유전자 비율을 확인하는 것입니다. 특정 세포에서 미토콘드리아 유전자 비율이 높으면 apoptosis 즉, 세포자멸사된 세포일 가능성이 높습니다. 세포자멸사된 세포에서는 미토콘드리아 유전자만 남고 핵에 있는 유전자는 분해되어 버리기 때문입니다. 그만큼 세포가 스트레스 받거나 죽어가면 미토콘드리아 유전자의 상대적 비율이 높아지는 경향이 있어 미토콘드리아 유전자 비율(percent.mt)을 저품질 세포를 판단하는 지표로 활용하곤 합니다. 그 기준 또한 연구자 재량껏 시료의 특성과 도메인 지식을 기반으로 설정하면 됩니다. 일반적인 조직 시료의 경우 5% 혹은 10%를 임계값으로 정해 그 이상의 미토콘드리아 유전자 비율을 보이는 세포를 제거합니다. 이떄 세포의 총 RNA 수(UMI 수) 별 미토콘드리아 유전자 비율도 Seuare의 FeatureScatter을 이용해 시각화하여 필터링 기준을 정할 수 있습니다. 아래 왼쪽 이미지를 보면 총 RNA수가 적은 세포가 미토콘드리아 유전자 비율이 높은 경향이 있는 걸 확인할 수 있는데요. 미토콘드리아 유전자 비율이 높은 세포는 핵에 있는 유전자가 분해되어버린 세포자멸사된 세포일 가능성이 높기 때문으로 볼 수 있겠습니다.
또한, 기본적으로 세포마다의 총 유전자 개수(nFeature_RNA), 총 RNA 수(nCount_RNA), 미토콘드리아 유전자 비율(percent.mt) 각각을 바이올린 플롯으로 나타냄으로써 세포들의 분포를 확인해 필터링 기준을 정할 수도 있습니다. 다음과 같이 말이죠.
결국 QC 단계에서는 샘플의 특성, 연구 목적, 그리고 QC 지표들의 분포를 확인하여 합리적인 임계값을 정해 고품질 세포만 남도록 필터링을 하는 것이 중요합니다.
# 예시코드, merged_cov가 seurat 객체
merged_cov[["percent.mt"]] <- PercentageFeatureSet(merged_cov, pattern = "^MT-")
p <- VlnPlot(object = merged_cov, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"))
plot1 <- FeatureScatter(merged_cov, feature1 = "nCount_RNA", feature2 = "percent.mt") + xlim(0, 40000)
plot2 <- FeatureScatter(merged_cov, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") + xlim(0, 40000)
options(repr.plot.width = 12, repr.plot.height = 6)
plot1 + plot2
# 필터링
pbmc <- subset(merged_cov, subset = nFeature_RNA > 300 & nFeature_RNA < 3000 & percent.mt < 15)
# 잘 필터링되었는지 확인
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
options(repr.plot.width = 12, repr.plot.height = 6)
plot1 + plot2
2. Normalization (정규화)
정규화(Normalization)는 간단히 말해 부피도 다르고 얻어낸 방식도 조금씩 다른 세포들에 대해 '세포 간 비교'가 가능하도록 만들어주는 단계입니다. 즉, 기술적 변이에 의한 세포마다의 총 RNA 수의 차이(ex. 기기에 따른 증폭 편향, RNA 수율 차이), 그리고 유전자의 비선형적 발현 패턴(어떤 유전자는 너무 많이 발현하고, 어떤 유전자는 너무 적게 발현)을 보정해주는 단계입니다. Seurat에서는 기본적으로 ‘전역 스케일 정규화(global scaling normalization)’ 방법으로 ‘로그 정규화(LogNormalize)’를 사용하거나, SCTransform 기반 정규화를 적용할 수도 있습니다.
먼저 로그 정규화(LogNormalize)는 모든 세포에 대해 총 UMI 카운트를 일정한 스케일로 맞춘 뒤 로그(log1p) 변환을 수행하는 전통적인 방법입니다. Seurat에서는 각 세포의 총 UMI를 10,000으로 맞춰주는 라이브러리 사이즈 정규화(library size normalization)를 진행합니다. 다음과 같이 말이죠.
이후엔 해당 값에 대해 log1p 변환 즉, 로그 변환을 진행합니다. 이를 통해 샘플마다 읽힌 길이가 달라도 상대적인 발현도를 비교할 있게 됩니다. 그럼 왜 정규화를 진행하지 않으면 세포 간 비교가 어려울까요? 예를 들어 설명해보겠습니다. RNA 시퀀싱 결과 특정 유전자의 RNA가 10개 관측되었다고 가정해봅시다. 이때 총 RNA 수가 20개인 세포에서 해당 유전자의 RNA 수가 10개를 차지하는 경우와 200개인 세포에서 10개를 차지하는 경우는 발현 정도는 다르게 평가해야 합니다. 20에서 10은 절반이지만 200에서 10은 20분의 1 정도니까요. 따라서 우리는 이러한 세포에 대해 다른 기준으로 판단할수 있도록 각 세포의 총 UMI 수를 동일하게 맞추어 줌으로써 차이를 보정하여 ‘세포 간 유전자 발현’을 비교할 수 있도록 합니다. 이에 더해 로그 변환을 적용하면 특정 유전자가 과도하게 발현되거나 거의 발현되지 않는 세포 별 유전자 발현의 비선형적 패턴을 완화하여, ‘세포 간 유전자 발현’을 좀 더 정확하게 비교할 수 있습니다. 이러한 로그변환을 통해 특정 세포에서 워낙 발현량이 적은 유전자의 RNA 수가 10개에서 30개로 변한 게 고발현 유전자의 RNA 수가 1000개에서 1050개로 변한 것보다 유의한 차이가 있다는 걸 우리는 이해할 수 있게 됩니다(RNA 수를 극단적으로 적게했다는 점 양해바랍니다).
# 예시코드, pbmc가 seurat 객체
pbmc <- NormalizeData(
object = pbmc,
scale.factor = 10000,
normalization.method = "LogNormalize"
)
또 다른 널리 사용되는 정규화 방법으로는 SCTransform 정규화가 있는데요. 이 방법은 정규화, 변동유전자 선정, 스케일링 및 회귀를 한번에 수행하는 글로벌 통계 모델링 접근법입니다. 구체적으로, 생물학적 차이에 기인한 변이가 아닌 기술적 차이(ex. 기기 차이, 실험자 차이)에 기인한 변이를 음이항회귀모델(negative binomical regression model)로 추정하여 제거하고, '잔차(Residual)'를 정규화된 표현값으로 사용합니다. 즉, "이러한 세포라면 이 정도의 RNA 수를 가지고 있을꺼야"라고 판단한 '예측값'을 '관측값'에서 빼줌으로써 기술적 변이를 계산할 수 있습니다. 너무 간단하게 설명했는데 보다 복잡하게 계산합니다. 그 결과 기술적 변이가 강하게 감소되어, 생물학적 신호가 더 선명해집니다. 실제로, SCTransfrom을 거친 데이터는 단순 로그정규화 데이터보다 미세한 차이까지 주성분에 반영될 수 있어, 다음 과정에서 진행할 PCA에 더 많은 주성분을 유의하게 활용할 수 있다고 합니다. 물론, SCTransform은 계산량이 많지만, glmGamPoi 등의 가속화 패키지를 활용해 속도를 개선할 수 있다고 합니다.
Analysis, visualization, and integration of Visium HD spatial datasets with Seurat
Seurat
satijalab.org
어떤 정규화 방법을 선택할지는 데이터의 특성과 목표에 따라 달라질 수 있습니다. 보통 소규모 데이터나 빠른 프로토타입 분석에는 'LogNormalize'로 충분할 수 있지만, 배치 효과 제거나 고정밀 클러스터링이 필요할 땐 SCTranfrom이 권장되는 추세라고 하죠. 중요한 건 분석 방법을 선택할 때, 앞에서 무엇을 했고 뒤에서 무엇을 할 것인지 생각해서 선택해야 한다는 점이 아닐까 싶습니다. 멋대로 했다가 전체 흐름을 놓치고 데이터를 이상한 형태로 변환시켜버리면 곤란합니다. 암튼 LogNormalize와 SCTransform 모두 사용 목적은 세포 간 비교 가능한 발현량 스케일을 얻도록 하는 것입니다. 또한, 정규화 후에는 데이터를 로그 스케일로 다루게 되므로 이후 PCA 등의 차원 축소 단계에서 분포 가정에 부합하게 될 수 있다고 하네요(참고로 Seurat v5에서는 SCTransform v2가 기본 정규화로 설정되어 있다고 합니다). 결국 scRNA-seq 데이터 분석에서 '정규화'는 QC와 마찬가지로 필수라고 볼 수 있겠습니다.
3. Feature Selection (특징 선택)
정규화된 데이터에서는 높은 변이를 보이는 유전자(Highly Variable Genes, HVGs)를 식별해야 합니다. 사실 SCTransform에서는 이미 진행을 한 절차입니다. 그래도 정규화에 SCTransform을 항상 사용하는 건 아니니 짚고 넘어가겠습니다.
단일세포 RNA-seq 데이터는 평균 발현이 높은 유전자일수록 변이(변동성)도 높게 관측되는 경향이 있다고 합니다. 그말은 즉슨 단순히 분산이 큰 유전자를 고르면 발현량이 높은 유전자들이 상위에 오를 수 있다는 것입니다. Seurat에서는 이러한 특징을 기반으로 평균 대비 분산(variance-to-mean ratio)을 고려한 통계적 방법으로 실제로 세포 간 발현 차이가 두드러진 유전자를 식별합니다. 평균을 고려하는 이유는 평균 발현이 높은 유전자들이 평균 발현이 낮은 유전자들보다 상대적으로 작은 변화임에도 변동이 크게 보일 수 있기 때문입니다. 앞서 말했듯 20개 중 10개와 200개 중 10개는 다르니까요. 보다 직접적인 비유로 설명하자면 20만원 있는 사람에게 10만원 뺏는거랑 200만원 있는 사람에게서 10만원 뺏는건 타격이 다르죠(물론 개인 차가 작용하겠지만 단순한 가정으로).
이러한 방식으로 우리는 세포와 세포 간의 차이를 가장 잘 설명해줄 수 있는 기준이 되는 고변동성 유전자들(HVGs)을 선택합니다. 그렇다면 왜 우리는 높은 변이를 보이는 유전자를 따로 뽑아내고자 하는걸까요? 그냥 전체 유전자를 대상으로 분석을 진행하면 안될까요? 음 그렇게 하기엔 세포당 약 2만 리드(RNA 수)를 읽은 만큼 데이터가 너무 크고 컴퓨터 자원은 한정적이죠. 심지어 적은 정보를 가진 유전자들이 핵심 유전자들을 분석하는 데 방해(노이즈)가 될 수 있습니다. 다시 말해, 고변이 유전자만 따로 뽑아서 분석을 진행하는 이유는 데이터에 포함된 모든 유전자를 사용하면 계산 비용이 너무 높아지고, 노이즈가 많이 포함될 수 있기 때문입니다. 따라서 세포 간 변이가 큰 유전자만 우선적으로 분석에 사용하는 것이 효율적입니다. 참고로, Seurat에서 HVG를 뽑아내기 위해 사용하는 FindVariableFeatures 함수는 보통 2000개 가량의 HVG를 선택하며, 충분한 정보를 담으면서 노이즈를 줄이도록 그 기준은 다르게 설정할 수 있습니다. 당연히 SCTransform을 이미 사용했다면 이 과정을 거쳐 3000개 정도의 HVG가 뽑혔으니 FindVariableFeatures 함수를 실행할 필요가 없습니다.
종합적으로, 본 단계는 데이터 차원을 크게 줄이면서도 핵심 정보를 유지시켜주는 단계로, 전체 유전자보다 이런 고변동성 유전자에 초점을 맞출 때 단일세포 데이터의 의미 있는 군집 구조를 더 잘 포착할 수 있습니다. 허나, 특정 핵심 유전자가 변동성이 낮아 HVG에서 제외되었다면 후속 해석에서 간과될 수 있으므로, 분석자는 필요한 경우 해당 유전자를 추가 포함하거나 결과를 해석할 때 주의해야 한다고 하네요.
FindVariableFeatures 함수로 HVG를 선별한 뒤엔 보통 Scaledata 함수를 사용해서 Seurat 객체 내 유전자 발현 데이터를 정규화하는데요. 정규화는 아까 세포마다 UMI 카운드를 10000으로 맞춰주고 로그 변환하면서 해준거 같은데 왜 또 등장했을까요? 물론 앞서 Normalize data 함수로 다음과 같은 방법의 정규화를 해주긴 했지만, Scaledata 함수의 경우엔 Normalizedata 함수로는 보정되지 않는 유전자별 평균, 표준편차 차이 등을 맞춘다고 합니다. 보통 평균 0, 분산이 1이 되도록 데이터를 스케일링한다고 하네요. 물론 SCTransform에서는 스케일링까지 해주기에 본 절차는 필요가 없곘죠.
# 예시코드, 디폴드값 HVG 2000개 선택
pbmc <- FindVariableFeatures(pbmc)
pbmc <- ScaleData(pbmc)
참고로, 정규화 절차에서는 기술적 변이에 의한 차이 말고도, 연구 목적에 따라 세포 주기에 따른 세포의 유전자 발현 차이를 추가로 보정해주기도 합니다. scRNA-seq 분석을 통해 보고자 하는 세포 유형 간 유전자 발현 차이가 세포 주기 진행에 따른 세포 유형 간 유전자 발현 차이에 가려질 수 있기 때문입니다. 즉, 분석 시 세포 주기의 차이를 보고자 하는 것이 아닐 때는 세포 주기 효과(cell cycle effects)를 교란 변수(confounding factors)로 간주하고 이를 보정해주기도 합니다. 이를 위해 Seurat에서는 세포 주기 효과를 보정할 수 있는 기능을 제공합니다.
Analysis, visualization, and integration of Visium HD spatial datasets with Seurat
Seurat
satijalab.org
4. Dimensionality Reduction (차원 축소)
본 단계는 수천 개에 이르는 유전자(피처) 차원을 축소하여 주요한 변이 구조를 저차원 공간에 표현하는 단계인데요. 보통 1차적으로 주성분 분석을 수행한 후, 2차적으로 UMAP 또는 t-SNE와 같은 비선형 방법으로 임베딤됩니다.
먼저 주성분 분석(Principal Component Analysis, PCA)은 데이터가 주로 어느 방향으로 흩어져 있는지를 파악하는 기법으로, 분산이 가장 큰 방향을 찾아 주성분을 구해냅니다. 분산이 가장 큰 방향으로 직선을 그으면 그 직선이 다른 직선에 비해 해당 데이터들의 좌표와 가장 비슷할테니까요. 즉, 그 직선만으로 전체 데이터를 설명한다고 해도 비교적 전체 데이터를 잘 설명할 수 있습니다. 물론 주성분 분석은 하나의 주성분만 구하는 것은 아닙니다.
첫 번째 주성분을 구한 뒤에는 분산이 감소하는 순서로 또 다른 주성분들을 찾습니다. 이때 각각의 주성분은 이전 주성분들과 직교하는 방향으로 찾는데요. 이는 각 주성분이 이전 성분으로부터 독립적인 정보를 추출하도록 하기 위함입니다. 이러한 방식을 통해 주요 성분들을 파악하는 주성분 분석은 전체 샘플 간 차이를 가장 잘 설명하는 축(주성분)을 찾아 데이터 차원을 줄이는 데 활용됩니다. 이를 통해 고차원 데이터를 소수의 주성분으로 요약하여, 데이터 구조를 단순화할 수 있습니다.
PCA를 진행하는 대상은 앞서 구한 고변이 유전자들(HVGs)입니다. 이 HVG 발현 행렬에 PCA를 적용하면 서로 연관된 유전자들의 공변동 패턴을 요약한 '주성분(PC)'을 얻을 수 있습니다. PCA를 통해 주성분을 찾아 주성분 몇 개 만으로 차원을 축소를 하고 나면 노이즈를 억제하고 중요한 신호를 부각시킬 수 있습니다. 다시 말해, Seurat에서는 기본적으로 선택된 HVG를 대상으로 PCA를 수행하며, 그 결과로 각 세포에서의 유전자들의 PC(성분) 좌표 및 각 PC별 기여 유전자 목록을 확인할 수 있고, 이후 몇 개의 주성분만 남도록 차원을 축소하여 전체 HVG 데이터에서 중요한 정보만 부각할 수 있습니다. 그렇다면 몇 개의 성분을 남길지는 무엇을 기준으로 결정할까요? 보통 Elbow plot과 같은 여러 기준을 통해 유의한 PC의 개수를 결정한 뒤, 후속 분석에 사용합니다(통상 10~20개의 PC를 사용한다고 하지만 데이터가 클 땐 50 PC까지 쓰기도 하는 만큼 꼭 그렇진 않은 것 같습니다). Elbowplot의 경우 그래프가 성분 번호가 증가할 수록(x축이 증가할수록) 해당 성분의 표준편차가 점차 낮아지고 완만해지는 구간이 등장하는데요. 표준편차가 낮아지면 데이터 변이를 크게 설명하지 못하므로(분산이 작으니까) 완만해지기 이전의 성분들까지만 사용합니다.
# 예시코드, pbmc가 seurat 객체
pbmc <- RunPCA(pbmc)
ElbowPlot(pbmc)
PCA로 주성분을 구한 뒤에는 UMAP 또는 t-SNE를 이용한 비선형적 방법으로 차원 축소를 진행합니다. 사실 그 전에 통상적으로 진행하는 중요한 단계가 있어서 잠깐 짚고 넘어가겠습니다. 바로 'Batch Correction'입니다.
+. Batch Correction (배치 보정)
단일세포 RNA 시퀀싱(scRNA-seq) 실험에서는 실험 일자, 실험실, 사용 키트, 장비 등의 차이로 '배치 효과(batch effect)'가 발생할 수 있습니다. 예로, 같은 조직에서 유래한 시료라도 두 실험실에서 서로 다른 세포 분리 방법을 사용하면, 한 배치의 세포들은 스트레스 관련 유전자를 더 많이 발현하는 차이가 나타날 수 있다고 하죠. 결국 '배치 효과'란 여러 배치(샘플) 간에 기술적인 차이로 인해 유전자 발현 수준이 체계적으로 달라지는 현상을 말합니다. 이러한 차이는 의도된 생물학적 차이가 아님에도 불구하고 데이터에 반영되어, 실제 생물학적 신호를 가리거나 '가짜 군집 구조'를 만들어 기술적 변이로 인한 차이를 생물학적 변화로 인한 차이로 오인하게 될 가능성이 있습니다. 따라서 배치 효과를 교정하지 않으면 동일한 세포 유형이 배치에 따라 따로 따로 클러스터링 되어 버려 다운스트림 분석 결과를 왜곡할 위험이 있습니다. 특히, 단일세포 수준의 분석은 미세한 발현 차이를 포착하는 데 초점을 맞추기에 작은 기술적 변이도 크게 두드러져 배치 효과가 특히 문제가 될 수 있는데요. 따라서 배치 효과를 교정함으로써 실험 간 기술적 변이를 제거하고 공통의 생물학적 구조를 드러내는 것이 중요합니다. 즉, 배치 효과를 제대로 제거해야 동일한 세포 유형들이 배치에 무관하게 모여 클러스터링 되고, 배치 간에 가려져 있던 희귀 세포 아형도 발견할 수 있습니다. 결국 배치 보정(batch correction)은 단일세포 RNA-seq 분석의 신뢰성을 높이고 다양한 샘플을 합쳐 일관된 분석을 하기 위해 필수적인 전처리 단계죠.
GitHub - immunogenomics/harmony: Fast, sensitive and accurate integration of single-cell data with Harmony
Fast, sensitive and accurate integration of single-cell data with Harmony - immunogenomics/harmony
github.com
배치 보정에 자주 사용되는 방법으로 배치 효과 제거 알고리즘인 'Harmony'가 있는데요. 이는 Seurat과 연계하여 사용할 수 있는 독립적인 툴입니다. Harmony는 저차원 임베딩 공간(PCA 공간)에서 배치 효과를 제거하는 접근을 사용합니다. 본 프로세스를 간단하게 설명해보겠습니다. 먼저 여러 배치의 데이터를 모두 포함한 통합된 PCA를 계산하고(앞선 절차에 따라), 이 PCA 공간에서 각 세포를 클러스터링합니다. 이때 클러스터들은 배치가 섞인 상태로 구성되도록 유도되며, Harmony는 각 클러스터 내에서 '배치 다양성(batch diversity)'를 최대화하는 방향으로 k-평균 군집화를 수행합니다.
그런 다음 위의 이미지에서 보듯 클러스터의 중심을 기준으로 세포들의 좌표를 배치 간 차이가 줄어들도록 이동시키고, 이 과정을 반복적으로 수행해 배치 효과를 점진적으로 제거합니다. 이후 배치가 보정된 Haromy 결과를 기반으로 UMAP이나 클러스터링을 진행하면 됩니다.
Harmony 알고리즘은 다른 배치 보정 툴에 비해 특히 여러 배치를 통합하는데 계산 효율성과 확장성이 뛰어나다고 하는데요. 배치 효과를 제거하면서도 세포 유형 간 고유한 변이는 보존하는데 뛰어난 성능을 보이는 것은 물론, 배치가 여러 계층(ex. 환자, 실험일, 기술 플랫폼 등)으로 이루어진 복잡한 실험 디자인에서도 각 배치를 설명하는 메타 데이터를 여러 개 입력하여 동시에 교정할 수 있는 유용성이 있다고 하죠. 다만 Harmony의 결과는 배치 효과가 제거된 임베딩(저차원 공강) 형태로 제공되며 조정된 count matrix를 직접 산출하지는 않는다는 점에서 한계가 있습니다. 이는 클러스터링이나 시각화에는 문제가 없지만, 통합 후의 데이터로 차등 발현 분석(DEG analysis) 등을 수행할 땐 임베딩 결과가 아닌 count matrics를 활용해야 하기에 주의가 필요하다고 합니다. 물론, seurat 객체에 원본 데이터가 보존되어 있기에 차등발현분석 시엔 배치보정을 하지 않은 원본 카운트 데이터를 사용할 수 있습니다(참고로 seurat 객체의 데이터 구조를 쭉 보면 여러 layer들이 존재하고, 각각의 값들이 분석 시에 잘못된 방향으로 변형되지 않게 해둠을 알 수 있습니다). 이외에도 최근엔 딥러닝 기술의 발달로 scVI, scANVI 등의 자동인코더(VAE) 기반 방법을 통해 배치 효과를 최선의 방법으로 줄이고자 하는 노력이 이어지고 있습니다.
# 예시코드
# orig.ident는 배치 명이 담긴 열
library(harmony)
pbmc <- RunHarmony(pbmc, "orig.ident")
이제 주성분 분석 결과를 바탕으로 배치 보정까지 완료했습니다. 다음으로 해당 결과를 기반으로 비선형 차원 축소 알고리즘인 UMAP(Uniform Manifold Approximation and Projection)이나 t-SNE(t-distributed Stochastic Neighbor Embedding)를 적용하여 2차원 시각화를 진행합니다. 이 알고리즘들의 목표는 고차원에서 서로 가까운 거리에 있는 점(세포)들은 저차원에서도 가깝게, 멀리있는 점들은 멀게 배치하여 데이터의 내재적 구조를 보존하면서 저차원 공간으로 차원 축소를 진행하는 것입니다. 이때 두 방법 모두 지역적인 유사도를 최대한 보존하긴 하지만, 전역적인 거리는 희생하는 경우가 많으므로 플롯 상의 군집 간 거리 해석에 주의가 필요하다고 합니다.
최근에는 UMAP이 t-SNE보다 많이 사용되고 있는데요. UMAP이 글로벌(전역적인) 구조를 약간 더 잘 유지하면서도 속도가 빠르다는 장점이 보고되었기 때문입니다. 즉, t-SNE는 국소적 구조는 잘보존하지만 전역적 구조는 왜곡할 수 있기에 전역적 구조도 어느정도 보존할 수 있는 UMAP을 사용하는 것입니다. 물론 UMAP도 한계가 있는데요. 다음 글을 참고하시길 바랍니다.
[스크랩] 생물학자들이 무분별하게 UMAP plot을 쓰는 것에 대한 경고
https://simplystatistics.org/posts/2024-12-23-biologists-stop-including-umap-plots-in-your-papers/ Simply Statistics: Biologists, stop putting UMAP plots in your papersUMAP is a powerful tool for exploratory data analysis, but without a clear understandin
bio-kcs.tistory.com
정리하자면, 차원 축소 단계는 고차원 단일세포 데이터를 이해하기 쉽게 줄여주는 과정입니다. 다시 한번 절차를 정리해보겠습니다. 우선 PCA를 통해 HVG에서 주요 분산 성분(주성분)을 추출하고 노이즈를 제거합니다. 이후 주성분을 대상으로 배치 보정을 하고, 보정된 결과값을 기반으로 UMAP 혹은 t-SNE로 시각적 직관이 좋은 2D 좌표에 투영합니다. 이렇게 얻어진 저차원 표현은 이후 '클러스터링'이나 '계통 추론(trajectory inference)' 등의 분석 결과를 보기 좋게 시각화하는 데 활용될 수 있습니다.
8. Clustering (군집화)
클러스터링은 차원 축소된 표현을 바탕으로 '세포 군집(cluster)'를 식별하는 단계입니다. 사실 앞서 차원 축소 단계에서 UMAP(혹은 t-SNE)을 먼저 설명하겠지만, 코드로선 클러스터링을 먼저한 뒤 UMAP 임베딩 공간에 투여합니다. 아래 코드와 같이 말이죠. 사실 생각해보면 클러스터링 과정은 고차원 공간에서 서로 거리가 가까운 혹은 거리가 먼 세포들을 찾는 과정이기에, UMAP을 통한 저차원 임베딩보다 먼저 진행해야 하는 것은 당연합니다. 그렇다면 세포 간 거리는 무엇을 기준으로 판단할까요?
# 예시코드, pbmc는 seurat 객체, 20 성분 선택, harmony 결과 기반해서 클러스터링 진행
pbmc <- FindNeighbors(pbmc, dims = 1:20, reduction = "harmony")
pbmc <- FindClusters(pbmc, resolution = 0.5)
set.seed(1234)
# UMAP으로 차원 축소
pbmc <- RunUMAP(pbmc, dims = 1:20)
Seurat에서는 위의 코드에서 보듯 먼저 FindNeighbors 함수를 사용해 세포들 간의 유전자 발현의 유사성을 계산하여 최근접 이웃 그래프(KNN graph)를 구성합니다. 즉, 유전자 발현 패턴이 유사할 수록 거리가 가까운 세포, 유사하지 않을 수록 거리가 먼 세포라고 보면 되겠습니다. 이렇게 만들어진 셀-셀 그래프에서 Louvain 알고리즘이나 Leiden 알고리즘과 같은 모듈러리티 최적화 기법을 적용하여 커뮤니티(클러스터)들의 강도를 측정하는 커뮤니티 탐지를 수행합니다. 커뮤니티 탐지란 각 커뮤니티 내 존재하는 연결이 우연히 형성되는 연결보다 얼마나 더 밀집되어 있는지 보는 것입니다.
클러스터링 알고리즘들은 모듈리티값을 최대화하는 방향으로 네트워크를 분할하며 커뮤니티 구분 즉, 클러스터링을 진행하는데요. 여기서 높은 모듈리티 값은 각 커뮤니티(클러스터) 내부의 노드(세포)들 간의 연결은 매우 조밀하고, 커뮤니티 간 연결은 상대적으로 적게 되는 값을 의미합니다. 이러한 클러스터링은 Seurat의 FindClusters 함수가 구현하며, 기본값으로는 Louvain 알고리즘을 사용한다고 하네요(Leiden은 옵션으로 선택). 참고로, Leiden은 Louvain의 단점(비연결성 세포 처리 등)을 개선하여 나온 알고리즘으로, 더 견고한 군집화를 제공한다고 합니다. 그래도 보통 기본 Louvain으로도 충분하지만, 대규모 데이터나 군집 경계가 애매한 경우엔 Leiden이 유리할 수 있다고 합니다.
Louvain 혹은 Leiden과 같은 그래프 클러스터링에서 중요한 매개변수는 resolution(해상도) 값입니다. 이 값은 군집의 세분화 정도를 조절하며, 높게 설정할수록 더 많은 소규모 집단이 나오고 낮게 설정하면 적은 수의 거대 군집이 나오게 됩니다. 일반적으로 세포 수가 3천개 정도일 때 resolution을 0.4~1.2 사이로 조정하면 적절한 결과를 얻는다고 알려져 있으며, 데이터 규모가 커질수록 최적 resolution 값도 증가하는 경향이 있다고 합니다. 그렇지만 이 또한 연구자 재량껏 선택해야 할 하이퍼파라미터입니다. 즉, 연구자는 해상도를 변화시켜 가면서 군집의 수와 형태를 비교하고, 생물학적으로 의미있는 분류가 되는 지점을 선택해야 합니다. 이때 해상도가 너무 높아 너무 세분화되면 세포 유형이 쓸데없이 여러 군집으로 갈라질 수 있고, 너무 낮으면 이질적인 세포들이 뭉뜽그려질 가능성이 있습니다. 이러한 비교 기준 외에도 연구자가 연구 문제에서 보고자 하는 세포 유형(혹은 아형)이 등장하는 수준으로 해상도 값을 조정할 수도 있습니다(ex. memory CD8+ T cell). 뒤에서 잠깐 설명하겠지만, 서브 클러스터링을 차후에 진행하기 위해 처음엔 낮은 해상도로 설정해 뭉뚱 그려서 세포 유형을 파악하고, 특정 세포 유형을 따로 분리(subset)하여 본 클러스터의 하위 클러스터들을 파악할 수도 있습니다. 예로, classical monocyte만 따로 분리하여 해당 세포가 분비하는 사이토카인(유전자 발현을 통해 확인)에 따라 아형을 구분할 수 있습니다.
클러스터링을 진행하면 각 세포에는 군집 ID가 부여됩니다(seurat_cluster에서 확인 가능). 이를 UMAP 등 임베딩 플롯에 색깔로 표시하면 분류된 군집들을 시각화할 수 있습니다. 아래 이미지는 Seurat 홈페이지에 있는 이미지로, PBMC 데이터의 UMAP 시각화 결과, 9개의 클러스터를 색상으로 표현한 것입니다.
9. Cell Type Annotation (세포 유형 주석)
클러스터링을 한 뒤엔 각 클러스터가 어떤 세포 유형인지 파악하는데요. 이러한 절차를 'Cell Type Annotation'이라고 합니다. 즉, Cell Type Annotation은 클러스터링을 통해 구분된 각 세포 군집이 어떤 생물학적 세포 유형에 해당하는지를 설명하고 명명하는 단계입니다. 이 단계는 굉장이 주관적일 수 있고 표준화되어 있지 않아 까다롭습니다. 단순히 이름을 지정해주는 것이 아닌, 단일세포 분석에서 결과를 도메인 지식(생물학 지식)과 적절하게 연결해야 하기 때문입니다. 이러한 주석 방법은 크게 전문가에 의한 '마커 기반 수동 annotation'과 '데이터베이스/참조 기반 자동 annotation'으로 나눌 수 있습니다.
첫째, 마커 기반 수동 annotation은 숙련된 연구자가 문헌에 보고된 '세포 유형별 마커 유전자'들을 활용하여 각 군집의 발현 프로파일을 해석하는 방법입니다. 예로 샘플이 면역세포를 대상으로 하는 PBMC 샘플이라면 CD3E, MS4A1 같은 잘 알려진 마커 유전자의 발현을 확인해 T세포, B세포 등을 결정합니다. 추가로, T 세포 중에선 CD4 T 세포인지, CD8 T 세포인지, 그중에서도 naive인지 effector인지 memory인지도 마커 유전자를 통해 확인합니다. 보통 Seurat에서는 FindAllMarkers 또는 FindMarkers 함수를 사용하여 각 클러스터에서 다른 클러스터에 비해 높게 발현하는 유전자를 뽑은 후, 알려진 마커 목록과 대조하면서 세포 유형을 지정합니다(scRNA-seq을 진행한 논문을 보면 supplementary table에 활용한 마커 유전자를 적어두곤 합니다). 이 과정에서 featureplot 함수나 vlnplot 함수로 특정 유전자의 클러스터별 발현을 시각화하여 판단하게 됩니다. 허나 수동 주석은 전문가의 생물학적 지식을 동원하기에 신뢰도가 높을 수 있지만, 주관적이고 시간이 많이 소요되는 단점이 있습니다. 특히 전문가가 아닌 저와 같은 학생은 수동 주석을 하면서 제가 하고 있는게 맞는지 계속 의문을 들게 합니다. 암튼 보통 새로운 단일세포 데이터셋에서 30개의 군집을 일일이 수동 주석하려면 20~40시간이 걸릴 수 있다고 보고되었으며 6), 해석자에 따라 결과가 달라질 수도 있습니다. 무엇보다 잘 알려지지 않은 세포 아형의 경우엔 마커 정보가 부족해 식별이 어려울 수 있습니다.
# 예시 코드
# 1) FeaturePlot으로 확인
gene_ID <- c('CD3D', 'CD3E', 'CD3G', 'CD4', 'CD8A', 'CD8B', 'NCAM1', 'KLRB1', 'NKG7', 'CD19', 'MS4A1', 'CD38', 'CD14', 'CD68', 'FCGR3A', 'CD1C', 'LILRA4', 'CD34', 'HBB', 'PPBP')
for (gene in gene_ID) {
plot <- FeaturePlot(pbmc, feature = gene)
file_name <- paste0(gene, "_FeaturePlot.png")
png(filename = file_name, width = 800, height = 800)
print(plot)
dev.off()
}
2) VlnPlot으로 확인
p <- VlnPlot(
object = pbmc,
features = marker.genes,
group.by = "seurat_clusters",
assay = "RNA",
stack = TRUE,
flip = TRUE)
p
이러한 수동 작업의 한계를 보완하기 위해 세포 유형을 자동 추론하는 도구들이 개발되었습니다. 큰 틀에서 두 가지 접근이 있는데요. '기존 마커 유전자 세트에 대한 스코어링을 하는 접근'과 '참조 데이터세트와 비교하는 접근'이 있습니다. 첫번째 접근법은 미리 구축된 세포 유형별 마커 리스트에 현 클러스터의 발현 데이터를 대조하여 가장 일치하는 세포 유형을 찾는 방식입니다. 간단하죠. 그저 마커 유전자 하나하나를 우리가 직접 클러스터별로 비교할 필요 없이 컴퓨터가 각 세포 유형에서 마커 유전자의 평균 발현량을 확인하여 점수를 매겨줍니다. 두번째 접근법은 SingleR과 같은 툴로 대표되는 접근법인데요. 참조로 지정된 대규모 싱글셀 데이터(ex. 면역세포 아틀라스 등 이미 세포 타입이 라벨링된 데이터)의 표현 패턴과 현 샘플 세포들을 비교합니다. SingleR의 경우에는 라벨린(annotation)된 참조 세트와 현 데이터의 각 세포 클러스터의 유전자 발현을 비교해 현 데이터의 각 세포 클러스터에 가장 유사한 타입 라벨을 부여한다고 합니다. 이외에도 머신 러닝 분류 접근도 있는데요. 대표적인 예시가 바로 CellTypist입니다. CellTypist는 방대한 참조로 훈련된 분류 모델을 적용해 신규 세포들을 예측한다고 합니다. 올해 초 인턴 과정에서 scRNA-seq 데이터 분석을 할 때 CellTypist를 굉장히 요긴하게 사용했던 기억이 납니다.
[29일차] CellTypist :: scRNA-seq, cell-type annotation 자동화 툴
이번 글에선 최근에 개별세포분석 리뷰논문을 읽으면서 관련 툴들을 노션에 정리하다 찾게 된 흥미로운 툴인 CellTypist에 대해 소개해보겠습니다. 대학원 준비생의 노트 ( 24. 12 ~ ) | Notion'과학
tkmstudy.tistory.com
이처럼 자동 주석 도구들은 한번 설정으로 전체 데이터셋을 빠르게 주석할 수 있어 효율적이고 사람마다 결과가 달라지지 않아 재현성이 높다고 합니다. 실제로 많은 파이프라인에서 클러스터링을 진행한 뒤 SingleR, CellTypist, Azimuth 등의 도구로 빠르게 세포 타입을 할당하고 있다고 하죠. 그렇지만 자동 주석 도구들도 항상 완벽하지는 않으므로 최종 검증과 보정은 결국 연구자의 몫이라고 볼 수 있겠습니다. CellTypist를 돌렸을 땐 각 클러스터의 세포유형 라벨과 클러스터별 세포 유형 예측 점수들을 파악할 순 있었지만, 어떻게 그런 결과를 산출했는지 직관적으로 이해하긴 어렵더군요. 물론 수동 주석은 더 헷갈렸지만요. 그래도 자동 주석의 경우 참조 데이터에 존재하지 않는 새로운 세포 상태나 유형이 있다면 잘못된 라벨이 붙을 수 있고, 유사한 타입 간엔 혼동이 생길 수 있으니 주의해야겠습니다. 따라서 자동 라벨링 후에도 주요 마커 발현을 한번 더 점검하고, 필요한 경우 수동으로 수정하는 등 수동주석과 자동주석을 함께 진행하는게 신속하게 분류하면서도 주석의 신뢰도를 높이는데 좋을 듯합니다. 여기서 annotation을 주석이라고 적긴 했지만, 그냥 거의 annotation이라고 부르는 것 같습니다.
10. Downstream Analysis (후속 분석)
cell type annotation까지 완료했다면, 이제 가장 중요한 절차가 기다리고 있습니다. 우리는 이제 연구 가설을 입증하기 위한 '다운스트림 분석(Downstram Analysis)' 즉, 다양한 심화 해석을 진행해야 합니다. 다운스트림 분석에는 여러 가지 방법이 있는데요. 어떤 분석을 진행할 지는 연구 목적과 실험 디자인에 따라 달라질 수 있습니다. 저는 그 중 딥리서치가 소개한 대표적인 몇 가지 방법만 설명하겠습니다.
10-1) 차이 발현 분석 (Differential Gene Expression Analysis, DGE analysis)
DGE 분석(혹은 DEG 분석)은 서로 다른 조건이나 세포 유형(클러스터) 간에 어떤 유전자들이 유의하게 발현 차이를 보이는지 통계적으로 검정하는 방법입니다. 예로, 세포 유형 A와 B를 비교하여 클러스터 특이적 마커 유전자를 찾거나, 특정 질병의 중증도별(mild, severe)로 동일 세포 유형 간에 유전자 발현 차이를 비교할 수도 있습니다.
Seurat의 FindMarkers 함수는 기본적으로 윌 콕슨 순위합(Wilocoxon rank-sum) 검정을 수행하여 클러스터 간 혹은 조건 간 유의하게 발현 차이가 나타나는 유전자를 파악합니다. 참고로 여기서 윌 콕슨 순위합 검정은 두 그룹 간의 분포 차이를 평가하는 비모수적 방법으로, 두 그룹에서의 세포마다 유전자 발현량 값을 모은 후 각세포에서의 각 유전자들의 순위를 매깁니다. 그리고 각 유전자들의 순위합(rank sum)을 각각 그룹마다 계산합니다. 예시를 들어보겠습니다(임의로 드는거라 맞지 않을 수 있습니다). 각 그룹마다 세포 두개씩만 있다고 가정해봅시다. 이때 첫번째 그룹 A에서의 첫번째 세포에선 유전자 'po'의 발현량 순위는 2위, 두번째 세포에선 'po' 유전자의 순위는 1위라면 첫번째 그룹에서 유전자 'po'의 순위합은 3일 것입니다. 두번째 그룹 B에선 한 세포에선 순위가 100위, 다른 세포에선 순위가 105위라면 순위합은 205겠죠. 그렇다면 두 그룹에서의 'po' 유전자의 순위합의 차이가 202로 큽니다(크다는 기준도 상대적으로 비교를 해야겠지만 크다고 가정해봅시다). 202가 유의한 순위합 차이를 보인다면, 두 그룹에서의 해당 유전자의 분포는 유의하게 차이가 난다고 볼 수 있겠습니다. 다시 말해, 한 그룹에서 다른 그룹보다 유의미하게 높은 혹은 낮은 순위를 보이는 유전자는 두 그룹 간 발현 차이가 있다고 볼 수 있습니다.
FindMarker 함수에서 윌콕슨 순위합 검정을 통해 계산한 차이 발현 결과는 각 유전자의 logfold change와 adjusted p-value로 제공되는데요. logfold change는 나머지 클러스터 대비 특정 클러스터, 혹은 조건 A의 클러스터 대비 조건 B의 동일 클러스터의 평균 유전자 발형량을 로그 변환을 하여 계산합니다. 아래가 그 예시인데요. 여기서 'Mean Expression'은 조건 별 동일한 특정 클러스터(세포 유형)에서의 특정 유전자의 발현량을 나타낸 것입니다.
이때 보통 log2를 기반으로 하기에, 로그 변환된 fold change 값이 1 이상이라는 것은 실제 발현량이 비교 대상에 비해 최소 2배 이상 차이남을 의미할 수 있겠습니다. 반대로 -1 이하라는 것은 비교 대상에 비해 실제 최소 2배 이상 평균 발현량이 적다는 의미겠죠. adjusted p-value는 일반적인 p-value의 정확도를 높이기 위해 FDR(False Discovery Rate)을 보정해준 값으로, 보통 0.05 미만이면 통계적으로 유의하게 차이가 난다고 판단할 수 있습니다. p-value는 이전 딥러닝 글에서도 정리했지만, 동일 모집단에서 왔다고 가정할때(귀무가설이 맞다고 가정할 때) 두 집단의 평균이 다를 확률입니다. 결국 p-value가 크면 두 집단이 동일 모집단에서 왔다고 추정할 수 있기에 통계적으로 유의한 차이가 나지 않는다고 볼 수 있고, p-value가 작으면 통계적으로 유의하게 차이가 난다고 볼 수 있습니다. 보통 0.05를 판단 기준으로 하지만 그 또한 임의의 설정 값이기에 맞지 않을 수 있으며, 이러한 한계를 보완한 게 FDR을 고려한 adjusted p-value라고 할 수 있겠습니다.
[39일차] 프로그래밍 & 수학 책 정리 01 : 딥러닝을 위한 수학
안녕하세요, 이번엔 최근에 완독한 '딥러닝을 위한 수학' 책을 정리해보려고 합니다. 본 책은 수학책이면서 수학을 프로그래밍으로 구현하며 여러 심층학습 개념들을 소개해주는 게 좋았는
tkmstudy.tistory.com
결국, FindMarker 함수 결과에서 특정 유전자의 logfold change가 1 이상에 adjusted p value가 0.05미만이라는 것은 비교 대상에 비해 유의하게 높게 발현하는 유전자라는 뜻이고, logfold change가 -1 이하에 adjusted p value가 0.05미만이라는 것은 비교 대상에 비해 유의하게 낮게 발현하는 유전자라고 볼 수 있겠습니다.
# mild와 normal의 유전자 차등 발현을 비교하는 코드
up_DEGs <- FindMarkers(object = nc,
ident.1 = "mild",
ident.2 = "normal",
group.by = "orig.ident",
only.pos = TRUE,
min.pct = 0.1, # min.pct는 해당 클러스터 내 세포들이 해당 유의한 수준으로 발현하는 비율
logfc.threshold = 0.25,
test.use = "wilcox")
up_genes <- rownames(up_DEGs)[
up_DEGs$avg_log2FC >= 1 & up_DEGs$p_val_adj < 0.05]
차후엔 유의하게 상향 조절되는 유전자들(up-regulated genes)을 대상으로, 해당 유전자들이 어떤 생물학적 경로에 관여하는지 Gene Set Enrichment Analysis(ex. GO & KEGG analysis)를 진행할 수 있습니다. 참고로 Enrichr이라는 통합 웹 툴을 활용하면 보다 편리하고 신속하게 각 유전자들이 어떤 경로에 관여하는지 여러 가지 기준들로 확인할 수 있습니다.
Enrichr
Recently, Enrichr was upgraded to support inputting a background. This can be done via the UI and the API. Another new feature that was recently added is the ability to export up to 100 annotated gene sets after performing a metadata search using the Term
maayanlab.cloud
사실 저는 해당 유전자들이 유의하게 관여하는 경로가 등장해도 이게 어떤 의미를 갖는지 종합적으로 판단하긴 아직 어려운 듯 합니다. 물론 여타 논문들을 보면 조건별 혹은 클러스터별로 유의하게 상향 조절되는 유전자를 대상으로 해당 유전자가 관여하는 경로나 기능을 파악하는 다운스트림 분석을 많이 하는 것 같습니다. 암튼 이렇게 유의하게 상향 조절 되는 유전자 혹은 하향 조절되는 유전자는 다음과 같이 VolcanoPlot으로 나타낼 수 있습니다.
10-2) 세포 궤적 추론 (Trajectory Inference)
Trajectory Inference는 세포가 분화, 발달, 세포 주기 진행 등을 통해 연속적인 상태 변화를 보이는 경우, 각 세포를 시계열상의 한 지점으로 배치하기 위해 pseudotime 개념을 도입하는 다운스트림 분석 방법입니다. 여기서 pseudotime이란 실제 시간이 아닌 세포가 잠재적 생물학적 변화를 따라 이동하는 과정을 ‘유사 시간축’으로 나타낸 값으로, 특정 생물학적 과정에서 세포가 어느 지점에 위치하는지를 상대적으로 추정하기 위해 계산된 지표입니다. 즉, 세포들의 발현 양상을 통해 발달 혹은 분화의 진행 정도를 가늠할 수 있게 해 주는 값이며, Trajectory Inference를 통해 이를 예측하고 시각화할 수 있습니다.
Trajectory Inference를 위한 대표적인 도구로는 Monocle 3나 Slingshot 등이 있으며, 차원 축소 공간에서 연속적인 경로(나무 가지 모양의 궤적)을 찾아내는 알고리즘을 사용합니다. 최근 들었던 KOBIC 강의에서 "Monocle의 핵심원리는 단일세포 RNA-seq 데이터를 저차원 공간으로 축소하면서, 세포들이 변화하는 주요 축을 찾고, 그 축들을 트리(그래프)로 연결하여 세포 간 가상 시간 관계를 부여하는 것"이라고 어떤 교수님께서 설명해주었던 기억이 나네요. 이러한 궤도 분석을 통해 우리는 클러스터링만으로는 놓칠 수 있는 동적 과정의 단계별 차이를 모델링할 수 있습니다. 예로, 조혈모 데이터에서 조상세포 -> 전구세포 -> 말단세포로 이어지는 분화 경로를 단일세포 수준에서 재구성하고, 가지(branch)마다 어떤 유전자 어떻게 변하는지 분석할 수 있습니다. 또한, 세포 재프로그래밍(reprogramming)이 특정 조건(ex. 조직 손상 or 야마나카인자 투입)에서 제대로 이루어지는지 확인할 때도 본 분석을 이용할 수도 있습니다. 이러한 단일세포 궤적 분석은 발달 생물학이나 세포 상태 전환 연구에 도움이 될 수 있지만, 결과가 알고리즘 및 파라미터에 민감하고 표준화되어있지 않으므로 충분한 검증이 필요하겠습니다.
최근에는 분화의 방향은 물론 속도까지 예측하 수 있는 RNA velocity 계산 툴들도 등장했는데요. RNA velocity 추론을 위해 분기점에 키가 되는 특정 유전자에 대한 unspliced RNA와 spliced RNA의 양을 시간대별로 정량화합니다. 만약 해당 유전자가 활성화되면 unspliced RNA가 spliced RNA보다 먼저 나타날 것이고, 다음으로 spliced RNA가 늘어나겠죠. 해당 유전자가 단백질로 발현되려면, 해당 유전자로부터 전사되어 만들어진 pre-mRNA에서 가변 스플라이싱 과정을 통해 인트론이 떨어져나가고 엑손의 조합으로 mRNA가 만들어지는 가공 과정을 거쳐야 하니까요. 결국 이를 시간대별로 각각 정량화하면 스플라이싱 속도를 추정할 수 있고, 이를 통해 분화 혹은 발달 속도를 추정할 수 있습니다. 이러한 속도는 플롯에서 화살표의 길이로 표현하며, 길이가 짧을수록 속도가 빠른 것이라고 합니다. 다음 링크의 이미지를 참고하시길 바랍니다.
Pumping the brakes on RNA velocity by understanding and interpreting RNA velocity estimates - Genome Biology
Background RNA velocity analysis of single cells offers the potential to predict temporal dynamics from gene expression. In many systems, RNA velocity has been observed to produce a vector field that qualitatively reflects known features of the system. How
genomebiology.biomedcentral.com
그 외에도 앞서 말했듯 클러스터별 마커나 DGE 결과 리스트를 가지고 특정 클러스터에서 어떤 기능적 경로가 특이적으로 활성화되어 있는지 분석할 수도 있고, 특정 기능과 관련된 유전자 세트(ex. GSEA) 혹은 유전자 시그니쳐(ex. fetal like 시그니쳐)가 어떤 세포 유형에서 풍부하게(enrich) 나타나는지도 확인할 수 있습니다. 또한, 세포 사이 의사소통 분석(Cell-Cell interaction), 다중오믹스 통합(ATAC-seq, spatial transcriptomics 등과 연계) 등 여러 확장이 가능합니다.
종합적으로, 다운스트림 분석 단계에서는 발현 데이터로부터 생물학적 의미를 해석하는 다양한 방법들이 활용됩니다. 정리하자면 단일세포 RNA-seq 데이터 분석은 세포 간 유전자 발현 패턴의 차이를 확인하고, 연속적 변화를 추적하며, 기능적 맥락을 해석할 수 있는 도구로써 단순한 클러스터링을 넘어 세포 기능과 상호관계에 대한 종합적인 그림을 그릴 수 있게 해줍니다. 이를 통해 우리는 왜 특정 질병을 앓게 되는지, 보다 구체적으로 특정 질병의 중증도가 어떤 세포의 기능들로부터 영향을 받는지, 그리고 개인에게 특정 질병에 취약하게 만드는 요인은 무엇이며, 왜 그 사람은 약물 저항성이 나타나는지 이해할 수 있습니다. 이는 질병 진단 및 예방, 그리고 예후 개선을 위한 전략에 활용될 수 있겠죠. 중요한 것은 scRNA-seq과 같은 기술은 빠른 속도로 발전하고 있으며, 앞으로 또 다른 기술의 발전이 과학의 발전을 가져올 수 있다는 점입니다. 표준화된 scRNA-seq에 대한 교재가 없는데에도 그 이유가 있겠죠(그렇지만 우리에겐 앞으로 정확하고 더 정교해질 것으로 기대되는 딥리서치가 있습니다).
이것으로 단일세포 전사체 분석 기술(scRNA-seq) 개요 및 프로세스 정리를 마무리 짓겠습니다. 많은 부분 딥리서치 내용을 참고했으며, 옳지 못한 부분이 있을 수 있다는 점 양해부탁드립니다. 감사합니다!
참고자료
1) Sung-Hun WOO, Byung Chul JUNG. The Workflow for Computational Analysis of Single-cell RNA-sequencing Data. Korean J Clin Lab Sci 2024;56:10-20
2) 챗GPT, 딥리서치
3) 10X Genomics, Cell Ranger's Gene Expression Algorithm, URL : https://www.10xgenomics.com/support/software/cell-ranger/latest/algorithms-overview/cr-gex-algorithm
Cell Ranger's Gene Expression Algorithm - Official 10x Genomics Support
A set of analysis pipelines that perform sample demultiplexing, barcode processing, single cell 3' and 5' gene counting, V(D)J transcript sequence assembly and annotation, and Feature Barcode analysis from single cell data.
www.10xgenomics.com
4) Seurat, URL : https://satijalab.org/seurat/articles/pbmc3k_tutorial
Analysis, visualization, and integration of Visium HD spatial datasets with Seurat
Seurat
satijalab.org
5) Harmony Tutorial, URL : https://portals.broadinstitute.org/harmony/articles/quickstart.html
Quick start to Harmony
portals.broadinstitute.org
6) Singleron, Decoding the Biological Meaning of Your Data: The Power of Accurate Automated Cell Type Annotation, 2023
The Power of automated cell type annotation - Singleron
Discover the differences between manual and automated cell type annotation, a crucial step in analyzing single cell RNA sequencing data.
singleron.bio
