RNASeq, 특히 싱글셀 RNASeq 같은 고차원의 high-throughput 데이터들은 noise가 많이 들어올 수 있기 때문에, PCA로 한번에 딱! 보려고 하면 막 겹치고 예쁘게 나오지가 않아서 논문에 내려고 하면 음... condition 별로 구분이 가긴가지만 확연하게 딱! 구분이 가지가 않아요.
이럴 때 쓰라고 있는 건 아니지만, PCA와 비슷한 분석 방법으로 PLS-DA라는 것이 있습니다. Partial Least Squares-Discriminant Analysis의 줄임말으로써, PCA는 데이터의 variance를 최대한 유지하면서 dimension reduction을 하려고 하는 반면에, PLS-DA는 사전에 그룹에 관한 정보를 미리 줌으로써 covariance를 최대한 유지하면서 분석하는 방법입니다. 자세한 계산방법이나 차이점 및 장단점은 다음에 다뤄보기로 하고 여기서는 이 알고리즘을 이용하여 플랏을 만들어 보겠습니다.
아주 간단합니다. 저희가 PCA에서 만들었던 AverageExpression 값을 그대로 써서 Mixomics 패키지를 이용해서 그려보겠습니다.
일단 mixomics 패키지를 깔고 불러옵니다.
그런 다음 아까 우리가 janitor 패키지를 이용해 0값을 제거해서 저장한 mySCT2 변수를 그대로 이용해서 PLSDA를 계산하고 그려줍니다.
그러면 아래와 같은 플랏이 나와요.
음... 아까 PCA랑 같이두고 비교해볼까요?
같은 데이터로 서로 다른 분석방법을 이용해 시각화를 했을 때에 왼쪽의 PCA랑 오른쪽의 PLS-DA랑 어떤게 condition간의 차이를 보기가 쉬울까요?
저 방법들을 다 쓰고 데이터베이스를 뒤지고, 그리고 또 최신 논문들을 찾아보아도 어려운 건 어려운 거에요 ㅠㅠ
아니 진짜 뭔가 확실하게 딱 이거다! 라는 것도 힘들고, 소프트웨어마다 차이도 있고 하니 어렵더라구요.
그래서, 요새 학계에서 유행하는AI에게 물어보았습니다.
일단 제일 유명한 chatGPT는 제외하고 제가 논문 쓸 때에 주로 쓰는 3가지를 써봤습니다.
chatGPT를 제외한 이유는 버전4로 바뀌면서 달라졌다고는 하지만, 없는 사실을 있는 듯하게 대답하는 경우가 많았기 때문입니다.
그래서 처음으로는 perplexity AI를 써봤습니다. UC Berkeley 출신의 팀이 만든 엔진으로, 논문 쓸 때에 특히 좋은데, 그 이유는 여러 학술 데이터베이스에서 레퍼런스를 가져옵니다.
처음에는 싱글셀 시퀀싱이 필요하다고 말하네요. 그래서 이거 싱글셀 시퀀싱 결과야 라고 말해줬습니다.
뭐 딱히 도움이 되진 않습니다. 그럴 거면 레퍼런스는 왜 가져온거니..
그래서 또 다른 셀 클러스터 마커들을 가져와봤습니다.
또 모른다고 하네요....
마지막으로 테스트 한 클러스터는:
NK 세포라고 말해주긴 합니다만. 아... 그건 이미 알고 있었어... 사실 NK세포의 어느 타입인지 알고 싶었던 거야.
Perplexity는 안쓰는 걸로..
다음은 Bing Chat 입니다. 마이크로소프트에서 만든 대화형 AI이며 검색에 도움을 주고 곧 마이크로소프트 오피스에도 포함될 예정이라고 하더군요. 평소에 검색할 때에 레퍼런스와 함께 결과를 가져오는 건 좋지만, 레퍼런스가 학술지가 아니라 블로그나 일반 뉴스사이트일 때가 많아서 그닥 좋아하진 않습니다. 하지만 검색을 할 때마다 마소 포인트를 주죠 ㅎ
앞서 Perplexity AI와 같은 유전자 마커들을 같은 순서로 질문해 봤습니다. 암튼 절대 모른대요. 알 수 없으니까 물어보지 마래요. 몇 개의 유전자들은 뭐 말해줄 수도 있는데 일단 자기는 모른대요.
Bing Chat도 패스하는 걸로...
마지막으로 구글에서 만든 Bard입니다. 신청자에 한해서 베타 테스트를 하게 해주는데 운좋게도 베타 테스터에 걸려서 요새 유용하게 쓰고 있습니다. 하지만, Bing Chat이나 Perplexity처럼 레퍼런스를 주진 않아요.
먼저 첫번째 데이터에 관해 물어봤습니다.
오오 쓸만하죠? B세포인 줄은 알고 있었는데 그 중 어떤거냐고 물어보니 플라즈마 세포라고 합니다. 그리고 각각의 유전자에 대해 왜 그런건지도 자세하게 설명을 해줍니다.
다음은 어떤 세포 종류인지도 모르는 데이터를 던져 주었습니다.
오 NK 세포래요. 사실 NK 세포인 줄은 알고 있었거든요. 그 중에 어떤 타입인지 궁금했었습니다. 그래서 다시 물어봤습니다.
NK세포 중에 CD56bright이래요. 음... 그럴듯한데? 근데 어떤 조직에서 나온건 줄 알면 더 정확하다고 말해서 어떤 세포인지 말해줬습니다.
그래서 말해줬더니 아까 자기가 말한 게 맞대요. 근데 사실 조금 곤란한게, 제가 다른 클러스터를 annotation 하면서 이거 말고 다른 클러스터가 아 이건 확실하게 NK CD56bright이다 라고 생각해 놓은 게 있었거든요. 그래서 다시 물어봤죠.
뭐 그럴 수도 있대요. 아 이거 헷갈리는데.. 하고 NK CD56bright라고 점찍어 놓은 클러스터의 마커들을 올려서 물어봤습니다.
그러니까 그건 MAIT 세포래요. 제가 틀렸고 자기가 맞다니까 뭐 어쩔 수 없죠. 암튼 다른 AI에 비해서 여러모로 자세히 알려주고 결론도 도출해줘서 일단은 구글 바드 최고!!
결론:갓구글 찬양합니다.
아, 마커들은 Seurat에서 클러스터별로 FindMakers 기능을 이용해서 찾아서 발현량 순으로 정렬하여 7-10개 정도 넣어봤습니다.
뭐 AI들이 백퍼센트 맞다고 할 순 없겠지만, 앞으로 한가지 더 테스트해볼 수 있는 방법으로 유용하게 쓰일 것 같습니다.
GSEA란 Gene Set Enrichment Analysis의 줄임말로 특정 유전자들의 집합이 있으면 이 유전자들이 어떠한 특성을 가지는 지 알아보는 분석 방법입니다. RNA 시퀀싱의 downstream analysis로 많이 쓰입니다. 간단하게 GSEA를 하는 방법은, 내가 원하는 Gene들을 p-value로든 Fold change로든 어떠한 점수를 기준으로 1위부터 꼴등까지 줄을 세웁니다. 그 후에 그 순위를 가지고 어떠어떠한 Gene Set에는 높은 순위의 유전자가 많이 들어있더라, 낮은 순위의 유전자가 많이 들어있더라 등을 확인합니다. 여기서 Gene set은 특정 pathway에 속하는 유전자들일 수도 있고, 공통된 기능을 가진 유전자 그룹일 수도 있고, 특정 질병에 노출되면 발현되는 유전자 리스트일 수도 있습니다. 사용자가 정하기 나름이죠. Broad Institute와 UC San Diego에서 이러한 유전자 세트, Gene set을 데이터베이스화 해서 정리해놓은 사이트가 있는데 그게 MsigDB입니다. 인간과 쥐의 유전자를 그룹별로 잘 정리해놓아서, GSEA에 유용하게 쓸 수 있습니다. https://www.gsea-msigdb.org/gsea/msigdb/
GSEA는 보통 그룹 간의 DEG (Differentially Expressed Gene), 즉 유의미하게 다르게 발현된 유전자들을 가지고 실행할 수가 있습니다. 이 때에 실행하는 GSEA는 내가 뽑은 DEG들이 어떠한 성격을 지니고 있는지 알아보는 것이겠죠. 예를 들어서, 두 환자 그룹에 한 그룹에는 비타민D를 1개월간 복용시키고 다른 그룹에는 플라시보를 복용시켜서 유전자 발현의 차이를 보는 경우에 두 그룹 간의 DEG를 구해 GSEA를 돌려 본다면 비타민D로 인하여 유의미하게 변한 유전자가 대체로 어떠한 그룹에 속하며 어떤 기능을 가지고 있는 지 알아볼 수 있습니다. 이 때, 이 DEG가 10개-20개 정도면 뭐 하나하나 살펴볼 수 있겠지만, 100개에서 500개 이렇게 나온다면 하나하나 보는 건 너무 많은 시간과 노력이 들겠죠. 그러면 GSEA가 이러한 시간을 아주 줄여줄 수 있습니다. 비슷하지만 다른 방법의 GSEA로는, DEG만 넣는 게 아니라 모든 유전자의 발현 정보를 넣어서 GSEA를 실행해볼 수도 있습니다. 이 때에는 입력 값이 아주 크겠지만, 위의 예제로 본다면 비타민D가 전반적으로 어떠한 영향을 미치고 어떤 pathway가 변화하는 지 종합적으로 알 수가 있습니다. 위의 두 방법 중에 어느 것이 옳고 어느 것이 틀리고 그러지는 않습니다. 다만, 내가 알아보려고 하는 질문에 따라서 전자 혹은 후자를 택하기도 합니다. 이 글에서는 Single cell RNA-Seq 분석을 할 때에 그 Downstream analysis로 Seurat 객체를 가지고 GSEA를 하는 방법을 간단히 기록하려고 합니다.
일단 먼저 Seurat을 불러오고 내가 이용하려는 데이터를 불러옵니다. 또한, 우리가 하려는 GSEA는 모든 유전자를 줄세워서 입력하는 위에 적은 후자의 방법을 쓰기 때문에, 유전자를 그룹 간에 비교한 값이 필요합니다. Seurat으로는 2만개 이상의 유전자를 모두 비교하는 데 엄청나게 오래 걸리기 때문에, presto 패키지에 있는 Wilcoxon rank sum test를 불러와서 사용하려고 합니다.
library(Seurat)
library(devtools)
install_github("immunogenomics/presto") #인스톨을 최초로 한 다음에는 이 줄은 생략합니다
library(presto)
#또 나중에 필요한 라이브러리 들을 여럿 불러와 보겠습니다.
library(msigdbr)
library(fgsea)
library(tidyr)
library(dplyr)
library(ggplot2)
library(tibble)
library(tidyverse)
library(data.table)
그런 다음, Wilcoxon rank sum test를 모든 유전자 대상으로 수행해줍니다. 제 Seurat 데이터는 myData3에 저장되어 있으며 저는 그 중 celltype.condition에서 특정 세포 (CD14+ Monocytes)의 질병 그룹(cALD)과 건강한 비교군 (Control) 간의 차이를 보고 싶습니다.
이제 앞서 말한 mSigDB에서 Gene Set들을 불러와야 합니다. 저는 가장 기본으로 유전자를 50개의 질병 및 생명 현상으로 나눈 Hallmark gene sets를 먼저 보려고 합니다. 그리고 이 불러온 Gene Set을 GSEA에 필요한 정보만 뽑아서 정리해 줍니다.
#msigdbr_show_species()
#위의 코드는 모든 가능한 생물 종을 표시해 줍니다.
#"H"로 Hallmark gene set을 불러 옵니다.
m_df<- msigdbr(species = "Homo sapiens", category = "H")
#불러온 Gene Set을 정리해 줍니다.
fgsea_sets<- m_df %>% split(x = .$gene_symbol, f = .$gs_name)
#불러온 Gene Set을 살펴봅니다.
head(m_df)
dplyr::count(myData.genes, group)
그런 다음에 제가 원하는 그룹을 뽑아서 logFC와 AUC score를 빼낸 뒤 정리해봅니다. 저는 CD14+ Mono_cALD 그룹에서 어떠한 pathway들이 더 혹은 덜 발현되는지 알아보고 싶습니다.
fgsea는 인풋으로 1) 미리 지정한 Gene Set 이 필요하구요 (m_df로 불러들여서 fgsea_sets로 포맷을 맞춤) 2) 유전자의 리스트와 각각 유전자의 순위 혹은 순위를 지정할 수 있는 점수가 필요합니다 (myData3에서 auc를 계산하여 myData.genes에 저장, 그런 다음 CD14Mono.genes에 순서대로 저장)
싱글셀 시퀀싱, scRNA-Seq 분석을 하다보면 가장 험난하고 어려운 단계가 각각 클러스터의 세포 종류를 지정하는 cell type annotation 인 것 같습니다. 저도 분석을 하면서 여러가지 방법을 시도해 봤는데, 그 기록을 남깁니다.
1. Seurat FindMarkers
Seurat에서 제공하는 가장 기본적인 메소드인 FindAllMarkers를 이용하면 각각 클러스터에서 어떠한 유전자가 많이 발현되는 지 알 수 있습니다. 그것을 이용하여 클러스터에서 마커들을 찾고 데이터베이스 혹은 논문에 비교해봅니다.
장점: 간단하다.
단점: 수동 비교 및 결정이라 시간이 오래 걸리고 헷갈릴 때가 많다
Seurat tutorial에서 제공하는 코드입니다. 이 코드를 실행시키면, 각각의 클러스터에서 다른 세포들보다 과발현된 유전자들을 정리해줍니다. 저는 주로 이 결과를 csv로 정리해서 엑셀로 열어본 다음 여러가지 데이터베이스에 넣어보기도 하고, 논문을 찾아보기도 합니다.
주로 찾아보는 데이터베이스:
a. PanglaoDB: 싱글셀 데이터가 방대하게 모아져있으며 유전자를 넣으면 그 유전자가 어느 세포에서 얼마나 발현되는지 알려줍니다.
레퍼런스는 celldex 레퍼런스를 쓰고 있으며, 본인의 데이터 종류에 따라서 여러가지의 레퍼런스를 불러올 수 있다는 것이 장점입니다. Mouse 싱글셀 시퀀싱 데이터를 분석 중인데 가장 간편한 방법이 SingleR이었습니다.
SingleR을 실행시킨 다음, 결과를 heatmap으로 그려볼 수 있습니다.
Heatmap 결과가 위와 같이 나오는 데 이걸 Annotation에 이용할 수 있습니다.
장점은 간편하고 실행 속도가 빠른 편이며 레퍼런스가 다양하다는 것이지만, 단점은 위와 보시다시피 정확한 세포 종류를 못 잡아내는 경우가 있으며, 서로 다른 클러스터가 수치 상으로는 같은 세포 종류가 되는 경우가 있습니다. 그리고 때때로 정말 이상한 결과가 나올 때도 있으니까 수동작업이 항상 필요합니다.
3. Azimuth
Azimuth는 NIH Human Biomolecular Atlas Project의 일환으로 개발된 cell type annotation 도구입니다.
여러가지 레퍼런스를 제공하고 있으며, 가장 큰 장점 중의 하나가 웹에서 구동이 가능합니다. 자신의 데이터를 업로드해서 계산해줍니다. 그만큼 Computing resource가 많이 들기도 하는데, 저는 주로 데이터 용량이 매우 큰 관계로 설치를 하여서 사용하고 있습니다.
장점은 위에서 언급한 것에 추가하여, 실행이 간편하며 특히 Seurat 패키지와 완벽하게 호환됩니다.
단점은 레퍼런스를 여러가지 제공하고 있긴 하지만 웹에서는 제공된 레퍼런스만 사용이 가능하며, 설치버전도 레퍼런스 만들기가 조금 복잡합니다. 또한 계산을 하면 UMAP이 나오는데 이게 좀 안맞을 때가 있습니다. 그래서 수동으로 하나하나 확인해주는 작업이 반드시 필요합니다.
Azimuth를 돌린 결과 UMAP입니다. 아, 그리고 때때로 이게 Seurat에서 작동하다보니 Azimuth를 돌리고 나면 내 원래 데이터를 건드릴 때도 있고, 나는 seurat_clusters에 나온대로 21개의 세포 종류로 annotation하고 싶은데 Azimuth는 resolution이 다를 때가 있습니다 (제가 잘 못하고 있는 거일 수도 있습니다.)
아, 그리고 설치버전은 좀 무거워서 데이터가 좀 큰 경우에 로컬머신에서 돌리기 버거울 수가 있습니다. 저는 주로 HPC 환경에서 돌리고 있습니다.
4. Reference Mapping
이건 딱히 적을 필요가 있나 싶은데 Seurat 패키지 자체에서도 Reference mapping을 제공하고 있습니다.
특히 자주 나오는 조직이나 샘플을 분석 중이라면 이 레퍼런스 맵핑이 유용하게 쓰입니다. 작동하는 방법은 Azimuth와 아주 비슷합니다. 그리고 결과도 Azimuth와 비슷한데, 가끔 Annotation이 다르게 나오는 경우가 있어서 참고용으로 유용하게 쓰고 있습니다.
아래는 제가 직접 쓰는 코드입니다:
Reference Mapping한 결과물입니다.
결론:
음.... 제가 싱글셀 분석을 하면서 가장 많은 시간을 쏟는 파트이자 가장 어려워하는 파트가 이 Cell type annotation인 것 같습니다. DEG나 clustering 처럼 딱 정답이 나오는 것이 아니고 레퍼런스에 따라 비교 방법에 따라 결과도 상이하게 나오는 경우가 있고 더군다나 마지막에는 결국 내가 결정해야하는 거구나 하고 결정해버리기 때문인 것 같습니다. 처음부터 끝까지 자동으로 해주는 Azimuth나 Reference mapping이 있긴 하지만, 그걸 맹목적으로 믿기에는 결과가 이상하게 나올 때가 많습니다 (위의 UMAP처럼요).
그래서 결국에는 어떤 방법을 쓰더라도 꼭 하나이상을 쓴 다음 그걸 수동으로 확인해서 검증하는 과정이 반드시 필요한 것 같습니다. 그리고 자신이 분석하고 있는 샘플의 특성을 잘 알고 각각 cell type의 특성이 어떤 것인지 잘 이해하고 있으면 더욱 도움이 되는 것 같습니다.
#scRNASeq 분석을 진행하던 도중에 UMAP을 보니 데이터에 #Doublet 이 제법 많이 있다는 사실을 알게 되었습니다. #UMAP 에서 보면 서로 다른 클러스터 사이에 브릿지같은 연결 다리가 보이죠? 그게 서로 다른 두 종류의 세포가 같은 방울 안에 묶여서 같이 시퀀싱 된 더블렛이라고 볼 수 있습니다.
이렇게 데이터 상에서 더블렛이 나타나는 경우에 다음 단계의 분석에 영향을 미칠 수도 있고, 클러스터링이 어려울 수가 있습니다. 그래서 #Seurat 에서는 기본적으로 feature number 혹은 gene number 등에 범위를 줘서 어느정도 doublet을 거르고 있습니다.
Seurat 기본 튜토리얼에서 너무 많거나 적은 RNA가 한 세포 안에 들어가 있는 걸 걸러주는 명령어, 출처: https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
하지만 데이터의 성격이나 품질에 따라서 이러한 스크리닝으로도 더블렛이 제거가 안되는 경우가 있습니다. 그럴 때에 따로 Doublet 을 찾아내는 프로그램을 사용하게 됩니다. 이러한 프로그램에는 여러가지가 있지만 저는 그 중에 #DoubletFinder 라는 프로그램을 선택해서 사용하였습니다. 이 소프트웨어를 쓴 이유는:
사용 설명이 비교적 자세하고 친절하게 되어있음 (그렇지 않은 프로그램들이 참 많죠 ㅠㅠ)
비교적 최근까지 업데이트가 이루어짐
Seurat 파이프라인과 연계가 간편함 (대부분의 다른 프로그램들은 SingleCellExperiment 스트럭쳐로 옮긴 다음에 실행을 해야해서, 데이터 구조를 변화를 시켜서 실행시킨 다음, 결과가 나온 후에 다시 Seurat 파이프라인으로 옮기는 과정이 들어가야 합니다.)
우리 학교 사람이 만들었어요.
이 소프트웨어는 아래의 링크에서 무료로 구하실 수 있습니다. R 패키지로 설치도 가능하구요.
DoubletFinder에 관한 배경이나 설명, 알고리즘은 위의 웹페이지에 자세하게 설명이 되어 있으니 넘어가고 실질적으로 사용하는 과정만 기록하겠습니다.
1. 일단 Seurat package를 불러들입니다. 만약 샘플끼리 이미 Integration을 했다면, 더블렛을 제거하고 다시 Integrate를 하는 것을 권장합니다. 아래의 R 코드는 Cellranger에서 count를 마치고난 결과를 읽어들여서 Seurat object로 만드는 것부터 시작합니다.
#다음 두줄은 나중에 합쳐서 분석을 위해 정보를 추가했습니다 myData$patient<-"p10" myData$visit<-"v3"
#익숙한 방법으로 mt 유전자들을 마크하고 걸러냅니다. 이때, Seurat 기본 튜토리얼에서는 5% 이상은 걸러내라고 적혀있는데, 실제 나중에 나온 실험 논문을 찾아보면 이 수치가 종마다, 조직마다 다르니까 찾아보시는 게 좋습니다 myData[["percent.mt"]]<-PercentageFeatureSet(myData, pattern = "^MT-") VlnPlot(myData,features=c("nFeature_RNA","nCount_RNA","percent.mt"),ncol=3) myData<-subset(myData,subset=nFeature_RNA > 200 & nFeature_RNA<4500 & percent.mt<15)
#이 시점에서 그래프가 나옵니다. pK를 계산해주는 그래프인데, 그래프만으로는 최대값 지점을 한눈에 알기가 어려우므로 다음 줄을 실행시킵니다 bcmvn_myData
pK 그래프, y값이 최고인 지점의 x값의 pK를 찾으면 됩니다.
pK 그래프 표. 맨 오른쪽 BCmetric이 가장 큰 줄의 pK를 찾아서 기억하세요.
4. 이전 단계에서 pK 최대값 지점을 확인한 후에 다음 코드를 실행시킵니다. 여기서 붉은색 숫자는 샘플 안에 Doublet이 얼마나 들어있는지를 예측하는 퍼센트입니다. 이 수치는 시퀀싱 테크닉마다 다르고, 샘플 내의 세포 갯수에 따라 또 다릅니다. 저의 경우에는 10X를 이용했는데 10X Genomics 홈페이지에 다음과 같이 퍼센트가 나와있습니다. 자신 샘플의 세포 갯수에 따라 조절해주세요.
#위의 코드를 실행하고 나면, myData의 구조를 다음 코드로 확인해봅니다 head(myData[[]])
#그 후에 아래의 pANN_0.25_0.005_763 을 자신의 데이터에 맞는 이름으로 고쳐주면 됩니다 myData <- doubletFinder_v3(myData, PCs = 1:30, pN = 0.25, pK =0.005, nExp = nExp_poi.adj, reuse.pANN = "pANN_0.25_0.005_763", sct = FALSE) head(myData[[]])
#그리고, DimPlot을 만들어 Doublet과 Singlet을 관찰해 봅시다 물론 DF.classifications 뒤에 오는 숫자는 알맞게 고쳐줘야합니다DimPlot(myData, group.by = "DF.classifications_0.25_0.005_763", raster = FALSE)
6. 이제 실제로 Doublet이 몇개인지 Singlet이 몇개인지 전체는 몇개인지 알아봅시다. 이런 수치들은 그때그때 적어서 기록으로 남겨두는 편입니다
