이번 글에서는 전 글에 이어 BRIC의 scRNA-seq data 분석법 글을 따라해보도록 하겠습니다.
[22일차] 학부생, BRIC의 scRNA-seq data 분석법 글 따라해보기 01 :: Seurat 불러오기 & Quality Control (QC)
이번 글에서는 scRNA data를 기반으로 R 환경 Seurat 패키지를 통해 개별세포분석을 해보고자 합니다. 이를 위한 data로는 전에 수강을 완료한 KOBIC 교육센터의 '예제 데이터를 활용한 단일세포 전
tkmstudy.tistory.com
앞서 pbmc 예제 데이터를 R 환경에 가져와서 필터링하였는데, 이번엔 해당 데이터에 대해 Normalization(정규화)을 해보겠습니다.
pbmc <- NormalizeData(pbmc)
데이터 정규화 과정은 세포 간 비교가 가능하도록 발현 데이터를 정규화하는 것으로, 세포 내 유전자 발현량 차이를 조정하는 것인데 보통 'log-normalization 방식(로그 변환 + 스케일링)'을 사용한다고 합니다.
개별세포분석에서의 Data Normalization에 대해서는 아래 글에도 정리해봤던 기억이 납니다. 보기 쉽도록 몇 가지 구절을 인용해보겠습니다.
[15일차] 논문 공부 :: DMEM / F12, Higly Variable Genes (HVGs), Data Normalization
이번 글에서도 저번처럼 현재 읽고 있는 논문0)을 읽으면서 중요하다고 생각되는 몇 가지 개념들에 대해 가볍게 정리하고 가도록 하겠습니다. [12일차] 논문 공부 :: Chromogenic IHC와 Immunofluroscen
tkmstudy.tistory.com
"Data Normalization에 적용되는 첫번째 표준화 방법은 'size-factor normalization'인데요. 본 과정은 실험에 관찰하고자 하는 건 biological variation이지 technical variation이 아니기에, sequencing depth를 통일하는 표준화 과정이라고 합니다. 여기서 sequencing depth는 sequence의 copy number로 볼 수 있습니다. 즉, 샘플의 library size를 표준화함으로써 각 샘플의 reads의 total number 차이를 제거하기 위한 과정으로 볼 수 있다고 합니다. 이를 위해 cell-specific factor로 각각의 UMI count를 곱하여 모든 세포 샘플들이 각각 같은 UMI counts를 갖도록 한다고 합니다."
- 출처 : 위 링크 참조 -
"또 다른 표준화 방법으로 'log transformation'이 있는데, 이는 expression outilers에 의한 영향을 최소화하고 높은 발현량의 genes가 분석값을 독점하지 못하도록 하는 방법이라고 합니다. 즉, log tansformation을 활용하여 low/medium abundance genes만 관찰하게 하는 것입니다(보통 UMI의 수가 적은 세포(건강하지 못한 세포)들이 불균형하게 high abundance genes를 갖는다고 합니다)."
- 출처 : 위 링크 참조 -
다음으로 변동성이 높은 유전자를 찾은 FindVariable Features( )함수를 사용하여 분석에 중요한 유전자를 식별해보겠습니다(모든 유전자를 사용하면 계산량이 늘어나고 노이즈가 증가할 수 있기 때문입니다). 보통 2000개의 변동성 높은 유전자를 선택한다고 하네요.
pbmc <- FindVariableFeatures(pbmc)
이번엔 각 유전자의 발현값을 평균 0, 분산 1로 스케일링(표준화)하는 ScaleData( ) 함수를 사용해보겠습니다. 이는 PCA 및 클러스티링과 같은 분석을 하기 위한 단계로, 세포 간 발현량의 차이를 제거하여 유전자 간 비교를 용이하게 합니다. 참고로 함수들에 대한 설명은 챗GPT의 말을 빌린 것이니 오류가 있을 수 있으니 양해부탁드리겠습니다.
pbmc <- ScaleData(pbmc)
이제 RunPCA( ) 함수로 데이터의 차원을 줄이는 주성분 분석을 실행할텐데요, PCA는 변동성이 높은 유전자의 발현 패턴을 요약한 주요 축을 생성하고 이를 통해 이웃도 찾고 UMAP도 그릴 수 있습니다.
pbmc <- RunPCA(pbmc)
참고로 여기까지 pbmc 데이터의 head값은 다음과 같습니다.
이제 PCA 결과를 기반으로 세포 간 유사성 즉, 이웃 관계를 계산하기 위해 FindNeigbors( ) 함수를 사용해보겠습니다. 여기서 dims = 1:10은 PCA의 처음 10개의 주성분을 사용하여 계산한다는 의미로 이를 통해 이후 클러스터링을 할 예정입니다.
글쓴이의 말에 따르면 데이터에 따라 20차원, 50차원 정보까지 사용하는 경우도 있고 주로 30차원 정도를 쓰는 것 같다고 하네요. 일단 전 잘 모르니 글쓴이가 사용한 10차원으로 사용해봤습니다.
pbmc <- FindNeighbors(pbmc, dims = 1:10)
다음으로 FindClusters( ) 함수를 사용해 이웃 관계 그래프를 기반으로 세포를 클러스터로 나눠보겠습니다. 이때 Louvain 알고리즘을 사용하며, resolution 파라미터는 클러스터의 세분화를 지정합니다(낮은 값 = 더 큰 클러스터, 높은 값 = 더 작은 클러스터).
본 글에서는 resolution parameter로 0.5를 사용했기에 저 또한 0.5로 사용했고, 이 또한 절대적인 기준은 없으며 논문에 따라 0.005부터 1.2까지 사용하는 것을 봤다고 합니다.
pbmc <- FindClusters(pbmc, resolution = 0.5)
여기까지 돌리니 pbmc 데이터 구조가 다음과 같이 변화된 것을 확인할 수 있었습니다. 'RNA-snn-res.0.5'와 'seurat_clusters'이라는 이름을 가진 열이 추가된 듯 합니다.
마지막 열이 세포별 클러스터인 것 같아서, table ( ) 함수로 빈도수를 보니, 다음과 같이 10개의 클러스터가 등장했습니다.
재미삼아서 해상도를 0.3으로 낮춰서 pbmc1이라는 변수를 만들고 똑같이 빈도수를 확인해보았습니다.
pbmc1 <- FindClusters(pbmc, resolution = 0.3)
table(pbmc1$seurat_clusters)
그랬더니 이번엔 7개의 cluster가 등장했습니다.
그렇지만 다음단계로 넘어가는 건 글쓴이가 했던것처럼 해상도 0.5인 pbmc 데이터를 활용해서 해보겠습니다.
다음 단계는 UMAP 시각화를 하는 단계로, PCA 주성분 중 처음 10개를 사용해 차원을 축소하고, 세포 간 유사성을 기반으로 시각화를 해봤습니다. 시각화에는 DimPlot( ) 함수를 활용했습니다.
pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")
시각화를 하니 다음과 같은 UMAP이 등장했습니다.
UMAP 결과를 분석하기 전, 필터링하고 분석에 용이하게 정제한 pbmc 데이터를 saveRDS( )함수를 활용해 'pbmc_tutorial1.rds'라는 이름으로 저장해보겠습니다.
saveRDS(pbmc, "pbmc_tutorial1.rds")
그러면 다음과 같이 앞서 전 글에서 지정했던 디렉토리에 'pbmc_tutorial1.rds' 파일이 만들어지며, 다음에 이것을 클릭하면 다시 불러올 수 있습니다.
다음 글에서는 UMAP plot에서 지정할 수 있는 여러 옵션들에 대해 알아보면서, UMAP을 활용해 분석할 수 있는 추가적인 요인들에 대해서도 알아보도록 하겠습니다.
치킨 응모까지 이제 이틀 남았네요. 열심히 달려보겠습니다. 감사합니다!
참고자료
1) INCOEDU, KOBIC 교육센터, 예제 데이터를 활용한 단일세포 전사체 데이터 분석
2) 쿼카, [당신의 논문 동료] scRNA-seq data 분석법 – Plot 4 종류부터 DEG까지, BRIC
