[전문기술] 단일세포 다중오믹스 분석

[전문기술] 단일세포 다중오믹스 분석

 

​

최근 단일세포 분석기술이 발달함에 따라 다양한 종류의 분자(mRNA, 단백질, DNA 접근성, 후성 유전학적 변형 등)정보가 얻어지고 있다. 개별 단일세포 분자 정보에서 더 나아가 동일 세포에서 두 가지 이상의 분자정보(mRNA+단백질발현량, DNA 접근성+mRNA 발현량)를 분석할 수 있는 단일세포 다중오믹스 실험 분석법들이 개발되고 있다. 이들 분석기법들은 단일세포의 다양한 분자 정보를 제공하지만, 측정 민감도에 있어 여전히 한계가 있다. 따라서 개별 세포에서 측정할 수 있는 분자 수가 한정적이며 많은 분자에 대해서 결측치들이 존재한다. 이러한 다량의 결측치(sparsity)에 강건한 데이터 분석법들이 필요하다.

​

​

본 보고서에서는 다양한 단일세포 다중오믹스 실험기법들과 데이터 분석 방법론을 소개한다. 또한 다중오믹스 데이터 분석은 각 분자 유형 데이터에 대한 결론의 신뢰도를 높이기 위해 동종의 여러 데이터 셋의 통합(integration)을 요구한다. 이러한 데이터 통합을 위해 사용되는 데이터 분석 방법론과 통합된 데이터를 대상으로 세포들이 이행하며 분화해 가는 것을 이해하기 위한 궤적분석(trajectory analysis) 방법론을 소개한다.

​

​

2. 단일세포 다중오믹스 분석

​

가. 단일세포 수준의 분자 분리

​

단일세포 다중오믹스 분석은 단일세포 단일오믹스 분석과 마찬가지로 개별 세포 수준으로 세포들을 분리하는 것으로 시작한다. 개별 세포 분리는 기계적인 방법 혹은 효소를 처리해 생존 세포들을 서로 떼어 놓은 후, 단일세포를 캡처하는 방식으로 진행된다. 단일세포 캡처 방법에는 1)레이저 포착 현미해부[1], 로봇식 현미조작[2] 등의 수십~수백 개의 세포를 캡처하는 저처리량 방법들과, 2) 형광 활성 세포 분리 분석 후 판상 분리, 미세 유체 플랫폼 방식[3] 등의 수만 개의 세포를 캡처하는 고처리량 방법들이 있다. 단일세포 다중오믹스 분석은 이렇게 캡처된 단일세포에서 여러 종류의 분자들을 분리, 추출해내야 한다. Genomic DNA(gDNA)의 경우 핵에 위치하며, 대부분의 mRNA는 세포질에 위치한다. 따라서 원형질막 특이적 용해 및 원심분리를 통해 핵을 분리하여 gDNA를 분리하고, 세포질 부분으로부터 mRNA를 분리하는 방식을 사용하기도 한다[4]. 그러나 이는 핵에 위치한 mRNA를 측정할 수 없다는 한계가 존재한다. 이를 극복하기 위해 oligo-dT가 코팅된 자성 입자를 이용해 mRNA를 캡처하여 mRNA를 gDNA로부터 분리하는 방식이 이용되기도 한다 [5]. 그러나 이 방법 역시 gDNA-mRNA 분리 과정에서 샘플 손실이 일어난다는 문제가 생긴다. 따라서 이후 gDNA-mRNA 분리가 필요 없는 방법이 개발되었는데[6], 이 방법에서는 분리가 되지 않은 상태에서 poly-dT primer를 이용해 mRNA를 역전사하여 cDNA로 만들고, gDNA와 cDNA를 동시에 증폭시키게 된다. 이후 이 product를 둘로 나누어 한쪽에서는 gDNA만 PCR을 이용해 증폭시키고, 다른 한쪽에서는 cDNA만 in vitro transcription으로 증폭시키는 방식으로 진행하게 된다.

​

​

나. 단일세포 유전체(gDNA)-전사체(mRNA) 분석

​

단일세포 다중오믹스 분석에서 gDNA와 mRNA의 분리 이후에는 각 분자 유형별로 단일세포 단일오믹스에서의 시퀀싱 방법들이 응용되어 왔다. 분리된 gDNA에 대해서는 단일세포전장유전체시퀀싱(singlecellWholeGenome Sequencing) 방법인 MDA[7],MALBAC [8],PicoPLEX(Rubicon Genomics PicoPLEX Kit) 등이 적용 되었으며, 분리된 mRNA에 대해서는 단일세포전사체시퀀싱(scRNA-seq) 방법인 Quartz-seq[9], Smart-seq [10], CEL-seq[11] 등이 적용되고 있다. 방법마다 서로 다른 전략을 취하게 되는데, 예를 들어 Quartz-seq의 경우 전사체의 3’ 말단을 측정하는 반면 Smart-seq은 전체 길이의 전사체를 측정하게 된다. 현재까지 다수의 유전체-전사체 다중오믹스 분석 방법들이 개발되었는데, scTrio-seq[4], G&T-seq[5], DR-seq[6], SIDR[12], TARGET-seq [13] 등이 있다(그림 1). scTrio-seq은 동일한 하나의 세포로부터 원형질막 특이적 용해 및 원심분리를 통해 핵과 세포질을 분리하고, 이로부터 분리한 gDNA와 mRNA 각각에 대해 scWGS protocol(MDA, PicoPLEX 등)과 Smart-seq2를 적용한다. G&T-seq은 oligo-dT가 코팅된 자성 입자를 이용해 gDNA로부터 mRNA를 분리하고, 분리한 gDNA와 mRNA 각각에 대해 scWGS protocol과 Smart-seq2를 적용한다. DR-seq은 gDNA와 mRNA를 분리하지 않은 상태에서 동시다발적 MALBAC-like quasilinear preamplification을 수행하고, gDNA와 cDNA 각각에 대해 MALBAC과 CEL-seq을 적용한다. 그러나 DR-seq은 다양한 scWGS protocol을 사용하는데 제약이 있으며, 전체 길이의 전사체를 측정할 수 없다는 한계가 있다. 이외에도 저삼투압 용해를 통해 핵을 보존하면서 세포액 RNA만을 분리하는 SIDR, 단백질 분해 효소를 통해 gDNA와 mRNA 분리의 효율을 높인 TARGET-seq 등의 방법이 있다.

​

 

​

다. 단일세포 전사체-후성유전체(DNA 메틸화, 히스톤 변형) 다중오믹스 분석

​

DNA 메틸화, 히스톤 변형, 염색질 접근성은 유전자 발현에 영향을 주며, 단일세포 수준에서 이들을 측정하려는 방법들이 개발되어왔다. DNA 메틸화의 경우 단일세포 bisulfite 시퀀싱(scBS-seq) 방법들을 이용하게 되는데, scRRBS[14], scWGBS[15], snmC-seq[16], sci-MET[17] 등이 있다. 단일세포에서 DNA 메틸화체와 전사체를 동시에 측정할 수 있도록 개발된 최초의 다중오믹스 방법은 scM&T-seq이며[18], 이 방법에서는 G&T-seq 과정에서처럼 gDNA와 RNA를 분리하고 증폭한 후, 증폭된 gDNA에 대해 scBS-seq 방법을 적용하여 DNA 메틸화체 데이터를 생산하게 된다. 이후 개발된 scMT-seq은 micropipetting으로 단일세포 용해물로부터 핵을 분리하고, scRRBS와 Smart-seq2를 이용해 DNA 메틸화체와 전사체 데이터를 각각 생산하는 방식으로 진행된다[19]. 또한, scTrio-seq은 유전체, DNA 메틸화체, 전사체를 동시에 측정하는 방법인데[4], scRRBS를 이용해 DNA 메틸화체 데이터를 생산하게 된다. 그러나, 이러한 방법들은 여전히 bisulfite 처리에 의해 DNA 분해가 유도될 수 있다는 한계점이 존재하며 [20], 추후 이를 극복하기 위한 방법들의 개발이 필요하다고 할 수 있다. 단일세포에서 히스톤 변형을 측정하는 방법으로는 Drop-ChIP 시퀀싱이 있다[21]. 이 방법에서는 먼저 미세유체 장비를 이용해 단일세포, 용해제, micrococcal(endo-exo) nuclease를 하나의 droplet에 가두어 뉴클레오좀을 생성한 후, 이를 세포 특이적 바코드를 함유하고 있는 droplet과 일대일로 병합하여 바코딩된 염색질 fragment들을 생산한다. 이후 이 fragment들을 취합한 후 ChIP-seq을 진행하여 히스톤 변형 위치를 동정하게 된다. 그러나 아직은 이 방법을 이용해 세포 당 측정되는 히스톤 변형 위치가 적으며, 이 방법을 적용한 단일세포 다중오믹스 분석 방법도 거의 개발되어 있지 않다. 염색질 접근성 프로파일링은 전사조절인자의 결합 위치(promoter, enhancer등) 예측 분석을 위해 사용된다. 단일세포에서의 염색질 접근성 분석 방법으로는 scDNase-seq[22], sci-ATAC-seq[23], scATAC-seq[24], NOMe-seq[25], scMNase-seq[26] 등이 개발되어 있다. 이 방법들을 이용한 단일세포 다중오믹스 분석 방법도 개발되었으며, sci-CAR[27], SNARE-seq[28], scNMT-seq[29] 등이 있다. 단일세포에서 gDNA 분리 후, 이들 방법을 적용하여 염색질 접근성을 측정함과 동시에 분리된 mRNA에 대한 전사체 분석법을 적용하여 염색질 접근성-전사체 다중오믹스 분석이 활발히 진행되고 있다.

​

​

라. 단일세포 전사체-단백체 다중오믹스 분석

​

단백질의 생물학적 중요성에도 불구하고, 단백질은 DNA나 mRNA(cDNA)와 같이 증폭시킬 수 없다는 특성 때문에 단일세포 단백체 분석을 통해 측정 가능한 단백질 수는 여전히 제한적이다. 단일세포에서 전사체와 단백체를 동시에 측정하는 방법들에는 PEA/STA[30], PLAYR[31], CITE- seq[32], REAP-seq[33], RAID[34] 등이 있다(그림 2). PEA/STA 방법은 DNA oligonucleotide가 부착된 항체를 이용해 단백질을 DNA oligonucleotide와 연결시키는 과정과 역전사를 통해 mRNA를 cDNA로 전환하는 과정을 진행한 후, DNA oligonucleotide와 cDNA를 각각 분리하고 증폭시켜 이를 qPCR이나 시퀀싱으로 항체가 타겟팅하는 단백질의 양을 측정하는 방법이다. PLAYR 방법에서는 단백질에 대해 동위원소를 포함하는 항체로 표지하게 된다. 또한 mRNA는 PLAYR probe를 결합시키고 이후 PLAYR probe pair들이 RNA 특이적 insert-backbone oligonucleotide의 결합 부위를 제공함으로써, rolling circle amplification을 통해 동위원소 표지된 probe와 결합되도록 한다. 이후 동위원소 표지된 mRNA와 단백질의 발현을 mass cytometry로 측정하게 된다. 최근에 개발된 CITE-seq과 REAP-seq의 경우 oligonucleotide가 표지된 항체를 이용하여 세포 표면 단백질과 mRNA를 동시에 측정하는 방법들이다. CITE-seq 방법에서는 먼저 세포 표면에 존재하는 표적 단백질을 PCR handle과 항체 식별 바코드, poly-A tail을 가지고 있는 DNA oligo- nucleotide가 결합된 표적 특이적 항체와 결합시킨다. 이후 세포들을 oligo- dT primer를 포함하는 bead가 들어있는 droplet 안에 가두고 droplet 안에서 세포 용해가 일어나게 함으로써, bead가 mRNA와 항체에 결합되어 있는 DNA 바코드를 캡처하게 하게 한다. 이후 이들을 역전사와 PCR을 통해 증폭시켜 library를 만든 후 시퀀싱을 수행하게 된다. REAP-seq도 이와 비슷한 방법이며, bead에 결합되는 바코드의 구성만 약간 다르다. 세포 표면 단백질만 측정하는 CITE-seq이나 REAP-seq과 달리 RAID 방법의 경우 세포 내부에 존재하는 단백질이나 인산화 단백질도 측정이 가능하다. 이 방법에서는 세포막 투과화 이후 단일세포들에 존재하는 세포 내 단백질들을 RNA 바코드가 결합된 항체로 면역염색하게 된다. 이후 세포들을 CEL-seq2-compatible primer를 포함하는 well plate에 sorting한 후 reverse crosslinking을 통해 RNA를 추출하고 역전사시켜 cDNA로 전환하게 된다. 이러한 방법들은 수천 개의 세포들에 대해 수십 개의 단백질들을 측정할 수 있다.

​

 

 

3. 단일세포 다중오믹스 데이터 분석

​

가. 단일세포 단일오믹스 데이터 분석 방법

​

단일세포 단일오믹스 분석은 단일세포 수준의 유전체 변형(돌연변이, 복제변이 등), DNA 메틸화, 염색질 접근성, mRNA 또는 단백질 발현 등의 정보를 제공해주며, 데이터 종류 및 분석 목적에 맞게 다양한 분석 방법들이 개발되어왔다. 예를 들어, scRNA-seq 데이터의 경우 세포군, 조절 네트워크, 세포 trajectory 등을 동정하기 위한 분석 방법들이 개발되어 있다. 먼저 세포군 동정은 세포들을 유전자 발현의 유사도에 따라 군집화하고 각 세포군에서 유의하게 높은 발현을 보이는 마커 유전자들을 선정하는 과정을 의미한다. 세포 군집화에 앞서 보통 데이터 차원 축소를 수행하게 되는데, 이때는 주성분 분석이 보통 사용되며 최근 들어서는 결측값을 고려한 신경망 기반 모델이 사용되기도 한다. 이후의 세포 군집화는 세포 간 거리 측정 기준(유클리드 거리,상관관계 등)과 군집화 알고리즘 종류(k-means, 계층적, 그래프 기반 군집화 등), 유전자와 세포의 동시다발적 군집화가 가능한지 여부 등에 따라 방법들이 분류될 수 있다. 자주 사용되는 분석 툴로는 Seurat[35], pcaReduce[36], SC3[37], BackSPIN[38], SNN-cliq[39] 등이 있다. 다음으로, 마커 유전자 간의 동시 발현 패턴에 기반하여 조절 관계를 예측하는 방법들이 사용되기도 한다. scRNA-seq 데이터의 경우 zero inflation이 존재하고 많은 세포 수로 인해 데이터의 차원이 높기 때문에 조절 네트워크를 추론하는 데에 어려움이 있다. 이를 해결하기 위해 네트워크 추론 방법들은 zero-inflated 분포 모델링, 데이터상의 조절 관계의 이산화 (binary value), 세포 및 유전자 발현 유사도 기반 군집화 등을 수행하게 된다. 네트워크 추론을 위해 자주 사용되는 방법으로는 SCNS toolkit [40], inferenceSnapshot[41], SCODE[42], SCENIC[43] 등이 있다. 마지막으로, 세포 분화 trajectory와 같은 세포 trajectory 분석을 위한 방법들도 사용되고 있다. 초기의 trajectory 추론 방법들은 trajectory의 위상(선형, 두 갈래형, 순환형 등)을 고정한 상태에서 세포들을 알맞게 정렬하는 것에 초점이 맞춰져 있었지만, 최근에 개발된 방법들은 trajectory의 위상 추론과 세포의 정렬을 동시에 진행하기도 한다. 자주 사용되는 방법으로는 Monocle[44], DPT[45], Wishbone[46], Waddington-OT[47] 등이 있다. scWGS 데이터의 경우 분석의 주된 목적은 단일세포 수준에서 유전자 복제 변이나 단일 염기서열 변이를 동정하는 것이다. 이를 위해 bulk WGS 데이터에서 사용되는 기존 분석 방법들을 단일세포 오믹스와 같은 낮은 유전체 coverage 데이터를 다룰 수 있게 변형하여 사용되는 경우가 많다. scWGS에서의 복제 변이 분석의 경우 Ginkgo[48], baseqCNV [49], SCNV[50], SCCNV[51], SCOPE[52] 등이 사용되며, 대부분 circularbinary segmentation 방식의 segmentation 과정을 차용하고 있다. 단일 염기서열 변이 분석의 경우 SCcaller[53], baseqSNV[49], MonoVar [54], SCAN-SNV[55] 등이 사용되며, 이들은 낮은 유전체 coverage와 높은 PCR 증폭 에러로 인해 발생하는 높은 allele coverage bias를 고려한 분석을 수행한다. 단일세포 후성유전체 데이터 분석의 경우 주로 열린 염색질이나 DNA 메틸화 자리 동정을 목적으로 한다. 단일세포 후성 유전체 데이터의 경우 bulk 데이터에 비해 유전체 coverage가 낮기 때문에 열린 염색질이나 DNA 메틸화 자리에 해당하는 peak을 정확히 동정하는 데에 어려움이 있다. 이를 극복하기 위한 방법으로 100여 개의 단일세포에 해당하는 데이터를 병합하여 기존 bulk 데이터에서의 방법을 이용해 peak을 먼저 동정하고, 이후 각 peak가 존재하는 단일세포를 찾아내는 전략을 사용하기도 한다. 일례로, scABC[56]가 이러한 전략으로 scATAC-seq 데이터에서 열린 염색질을 동정한다. 그러나 이러한 세포 데이터의 병합으로도 peak 동정이 어려울 수 있는데, 이러한 경우 인접한 영역에서의 signal을 병합하거나 서로 유사한 regulatory element를 포함하는 영역에서의 signal을 병합하는 전략을 사용하기도 한다. 일례로, chromVAR[57]가 scATAC-seq 데이터에서 열린 염색질을 동정할 때 이러한 전략을 사용한다. 또한 cisTopic[58], SCALE[59]와 같은 방법은 세포 데이터 병합과 영역 데이터 병합의 두 전략을 혼용하기도 한다.

​

​

나. 단일세포 다중오믹스 데이터 통합 분석 방법

​

단일세포 다중오믹스 데이터 통합 분석을 위해서는 앞서 기술한 단일세포 단일오믹스 데이터 분석 방법들을 확장하거나 조합한 방법들이 사용된다. 이러한 방법들은 크게 1) 단일세포 단일오믹스 데이터간의 상관관계 분석, 2) 한 종류의 오믹스 데이터 분석 후 다른 종류의 오믹스 데이터를 통합하는 분석, 3) 모든 종류의 오믹스 데이터를 한 번에 통합하여 진행하는 분석으로 구분 지을 수 있다(그림 3).

​

​

 

​

1) 단일세포 단일오믹스 데이터 간 상관관계 분석 방법

​

단일세포 수준에서 mRNA 발현량과 유전자 복제 변이 혹은 DNA 메틸화 정도 간의 상관관계를 보는 여러 연구에 적용되었다. 예를 들어, Angermueller et al.[18]에서는 scM&T-seq을 쥐 배아줄기세포에 적용하여 개별 유전자들에 대해 위치별로(promoter,distalregulatoryelement,genebody) DNA 메틸화 정도와 mRNA 발현량 간의 상관관계를 계산하였다. 이를 통해 DNA 메틸화와 mRNA 발현량 간의 음의 상관관계가 non-CpG island promoter에서 우세하게 관측되며, distal regulatory element에서는 양의 상관관계와 음의 상관관계가 모두 관측됨을 알아내었다. 상관관계 분석은 mRNA와 단백질 발현량 간에도 이루어지는데, 예를 들어 Peterson et al.[33]에서는 PBMC에 REAP-seq을 적용하여 면역 세포 마커들의 mRNA와 단백질 발현량 상관관계를 계산하였다. 이를 통해 mRNA와 단백질 발현량 간의 상관관계가 전반적으로 약하며, 낮은 mRNA 발현량을 갖는 마커들에 대해 단백질 정량이 mRNA 정량에 비해 더 민감함을 확인하였다.

​

​

2) 한 종류의 오믹스 데이터 분석 후 다른 종류의 오믹스 데이터 통합 분석

​

scRNA-seq이 기준이 되고 여기에 다른 종류의 오믹스 데이터를 통합하여 분석하는 경우가 많다. 왜냐하면, scWGS, scBS-seq, scATAC-seq 데이터는 낮은 유전체 coverage를 갖고, CITE-seq 데이터 역시 단백체 coverage가 낮은 반면에, scRNA-seq은 상대적으로 높은 전사체 coverage를 가지기 때문이다. 예를 들어, Cao et al.[27]에서는 sci-CAR를 쥐 신장 세포에 적용한 뒤 scRNA-seq 데이터를 먼저 이용해 10,727개의 세포들을 14개의 세포군으로 분류하였다. 이후 14개의 각 세포군에 특이적인 열린 염색질 위치들을 동정하였고, 이를 통해 세포군 특이적인 마커 유전자들의 발현에 기여할 수 있는 cis-regulatory element(전사조절인자 결합 위치)들을 동정하였다.

​

​

3) 모든 종류의 오믹스 테이터를 한번에 통합하여 분석

​

서로 다른 단일세포 오믹스 데이터가 비교 가능한 수준의 coverage를 가지는 경우 적용할 수 있으며, 그렇지 않은 경우에는 통합 분석이 더 높은 coverage를 갖는 데이터에 편향된 결과를 낼 수 있다. 분석을 위해 여러 matrix factorization 기반 방법들이 개발되어왔으며, LIGER[60], MOFA[61] 등이 있다. LIGER는 integrative nonnegative matrix factorization을 이용해 통합 분석을 수행한다. Welch et al.[60]에서는 LIGER를 쥐 전두피질 뉴런으로부터의 단일세포 전사체 및 단일세포 DNA 메틸화체 데이터에 적용하였다. 세포군을 정의할 때 mRNA 발현과 DNA 메틸화 정보가 모두 존재하는 유전자들과 둘 중 하나의 정보만 존재하는 유전자들을 사용하였는데, 이 중 두 정보가 모두 존재하는 유전자들을 이용해 DNA 메틸화를 통한 mRNA 발현 조절 효과를 예측하는 분석을 수행하였다. MOFA는 multiway matrix decomposition 방법을 사용하며, 이를 통해 하나의 factor loading matrix(세포들이 각 세포군에 어느 정도로 기여하는지)와 데이터 유형별 하나의 weight matrix(molecular feature들이 각 세포군에 어느 정도로 기여하는지)를 생성하게 된다. Argelaguet et al.[61]에서는 MOFA를 쥐 배아줄기세포 scM&T-seq 데이터 분석에 적용하였으며, 이를 통해 유전자들 중 mRNA 발현 또는 DNA 메틸화 정도가 각 세포군에 기여하는 정도가 큰 유전자들을 선별하였고 이로부터 각 세포군에서 중요한 기능을 하는 mRNA-DNA 메틸화 조절 관계를 예측하였다.

​

​

다. 단일세포 전사체 멀티데이터 셋 통합 분석

​

앞서 언급된 대로 단일세포 다중오믹스 데이터 통합 분석에서 단일세포 전사체 데이터는 가장 depth가 크고, 변이 및 후성유전학적 변형(DNA 메틸화, 히스톤 변형, 염색질 접근성 등)의 기능적 결과를 포함하므로, 세포군을 정의하고 그 특성을 파악하기 위해 가장 많이 도입되는 단일세포 오믹스 분석이다. 따라서 다중오믹스 데이터의 통합 이전에, 단일세포 전사체로부터의 결론의 신뢰도를 높이기 위해 같은 조건에서 생산된 여러 단일세포 전사체 데이터 셋을 통합 분석하는 것이 빈번하게 수행되고 있다. 주어진 데이터에 대해 통합 대상이 되는 데이터 셋은 Single Cell Portal(https://singlecell.broadinstitute.org/single_cell), Gene Expression Ominibus[62]와 Human Cell Atlas[63] 같은 공개 데이터 저장소에 축적된 단일세포 전사체 데이터들이 이용되고 있다. 여러 전사체 데이터 셋을 통합하면, 여러 데이터 셋에서 일관되게 관찰되는 고신뢰성 세포특성(cellular features)의 동정과 개별 데이터 셋에서 작은 수의 세포들이 모여 통계적 신뢰도가 있는 사이즈의 세포군 동정이 가능하다. 하지만 데이터 생성 단계에서의 차이로 인해 데이터 셋 간 또는 데이터 셋 내에는 관심의 대상인 생물학적인 변동 외에도 원치 않는 기술적 변동이 불가피하게 수반되는데, 이를 ‘배치 효과(batch effect)’라고 부른다. 배치 효과는 생물학적 요인에 의한 변동과 혼합되어 둘 사이의 구분을 어렵게 함으로써 세포 특성을 규명하는 것을 방해할 수 있다. 따라서 단일세포 전사체 데이터 셋을 통합 분석할 때 배치 효과를 교정할 필요가 있다.

​

​

1) 멀티데이터 셋에서의 샘플 배치의 정의

​

배치 효과는 주로 샘플 특성(기부자, 조직의 종류, 종, 질병 상태 등), 실험 프로토콜 및 시퀀싱 플랫폼의 차이에 의해 발생한다(그림 4A). 데이터 셋의 통합 분석을 위해 주된 배치 효과를 유발하는 요인을 주관적으로 결정한 후 이것에 기반해 특성이 비슷한 샘플들의 집합을 배치(batch)로 정의하는 경우가 많다. 예를 들어, 배치를 개별 기부자에서 유래한 샘플 집합으로 정의하거나[64-67] 서로 다른 프로토콜 또는 시퀀싱 플랫폼을 이용해 생산된 개별 데이터 셋으로 정의한 경우[68,69], 개별 샘플로 정의한 경우가 있다[70,71]. 두 가지 이상의 주요 요인이 있는 경우 각 요인에 대해 배치를 정의해 순차적으로 교정해 줄 수도 있다[68, 69]. 배치를 정의한 후에는 여러 데이터 셋의 발현 행렬을 병합하고 정규화한 뒤, 각 배치 내에서 가장 분산이 큰 유전자들(highly variable genes, HVG)을 선별하여 여러 배치에서 동시에 HVG로 선별된 유전자들을 데이터 통합에 활용할 수 있다(그림 4B). 일반적으로 HVG만 사용할 때가 전체 유전자를 사용할 때보다 원치 않는 변동을 제거하거나 세포군 간의 생물학적 차이를 규명하는 데에 효과적이라고 보고된 바 있다[72]. 위의 과정을 통해 얻은 병합되고 정규화된 HVG에 대한 발현 행렬을 배치 효과 교정과 그 이후의 세포 군집화 및 세포 종류 규명에 활용할 수 있다. 대표적인 배치 효과 교정 방법들을 배치 효과를 모델링하는 방식에 따라 분류하면 1) 선형 분해를 통한 배치 효과 교정 방법과 2) 축소된 차원에서 유사도에 기반한 교정 방법, 3) 변이형 오토인코더(variational autoencoder)를 사용한 교정 방법으로 나눌 수 있다.

​

 

​

➀ 선형 분해를 통한 배치 효과 보정 방법

​

일반화 선형 모델(generalized linear model)을 사용해 배치 효과를 모델링하는 방식은 bulk 전사체 데이터셋을 통합할 때 배치 효과를 제거하기 위해 처음 도입되었다. 예를 들어, limma(ref)의 ‘removeBatchEffect’ 함수[73]와 ComBat[74]이 이런 선형 분해 모델을 사용하였다. 선형 분해 방법은 대체로 알고리즘이 가장 단순하기 때문에 분석 속도면에서 장점을 갖는다[75]. limma의 ‘removeBatchEffect’ 함수는 배치 효과가 발현 평균에 미치는 영향을 모델링하기 위해 배치 회귀 계수 행렬과 배치 설계 행렬(design matrix)의 곱으로 표현되는 추가 배치 항(additive batch term)을 사용한다(그림 5A, Batch term matrix). 배치 설계 행렬은 샘플의 배치 정보를 원핫 인코딩 방식으로 표현한 행렬로, 같은 배치에 속하는 샘플들은 설계 행렬에서의 값이 같기 때문에 두 행렬의 곱으로 계산되는 배치 효과 또한 같다. ComBat은 분산에 대한 배치 효과를 모델링 하는 곱셈형 배치항(multiplicative batch term)을 추가적으로 포함한다(그림 5A,괄호안). 회귀 계수들을 추정하는 방식에 있어서 limma는 오차의 제곱 합을 최소화하는 선형 회귀 방식을 사용하는데 반해 ComBat은 샘플 사이즈가 적은 상황에서 과도한 교정을 예방하기 위해 empircal Bayes 방식을 사용한다. 하지만 두 방식은 기본적으로 같은 배치에 속한 샘플들에 대해서는 동일한 배치 효과를 적용해 보정해 준다는 공통점을 갖는다. 위의 접근 방식은 단일세포 전사체 데이터를 이용한 초기 연구에서 데이터 셋을 통합할 때 배치 효과를 제거하기 위해 사용되었다[65, 76-78]. 그러나 또 다른 선형 분해 모델인 ZINB-WaVE[79]는 단일세포 전사체 데이터의 고유한 특징인 zero inflation, overdispersion 및 bulk RNA-seq 데이터와는 다른 count 분포를 고려하여 단일세포 전사체(scRNA-seq) 데이터의 분석에 특화된 방법으로 개발되었다. ZINB-WaVE는 병합된 발현 카운트 행렬을 zero-inflated된 음이항 분포(ZINB)를 따르는 확률 변수로 정의한다. 음이항 확률 질량 함수의 평균(μ)에 대한 log μ 행렬과 결측 확률(π)에 대한 logit π 행렬 각각을 샘플 수준의 공변량 행렬과 유전자 수준의 공변량 행렬, 알려지지 않은 샘플 수준의 공변량 행렬로 분해하여, 샘플 수준의 공변량 행렬이 배치 효과를 모델링 하게 함으로써 배치 효과가 완화된 축소된 차원의 행렬을 얻을 수 있다(그림 5B).

​

 

 

​

② 축소된 차원에서 유사도에 기반한 배치 효과 보정 방법

​

위의 방법들은 같은 배치 내의 모든 세포들이 동일한 변동 요인에 대해 평균과 분산이 동등하게 변한다고 가정한다. 그러나 세포는 세포 상태(cell state) 등에 따라 다른 민감도를 가지므로 배치 간 및 샘플 간에 다른 변동 양상을 보일 수 있다. 배치 교정에서 이러한 세포 수준 공변량 문제를 해결하기 위해 다양한 방법들이 개발되었으며, 이러한 방법에는 CCA[80], MNN[81], fastMNN[81], Seurat의 ‘IntegrateData’함수[35], Scanorama[82], BBKNN[83], Conos[84], Harmony[85], DESC[86], LIGER[60], scMerge[87] 및 SAUCIE[88]가 있다. 이 방식들은 축소된 차원으로 세포들을 투영하는 것에서 시작하며, 이후 개별 세포 수준(그림 6A, 왼쪽)이나 세포 군집 수준(그림 6A, 오른쪽)에서 배치 간에 비슷한 발현 프로파일을 갖는 세포들을 찾아, 이러한 유사한 세포들이 축소된 공간에서 같은 분포를 따르도록 교정한다.

​

 

​

a. 개별 세포 수준 유사성 탐색

개별 세포 수준에서 배치 간에 유사한 세포들을 찾는 방법은 주로 최근접 이웃(nearest neighbor) 기반 방법을 사용하며, 이에는 MNN, Seurat의 ‘FindIntegrationAnchors’ 함수 및 Scanorama, BBKNN, Conos 등의 알고리즘이 포함된다. 두 배치(배치 1과 배치 2)를 정렬해야 하는 경우 MNN은 배치 1의 임의의 세포 i(그림 6B, 어두운 파란색 점)에 대해 배치 2에서 최근접 이웃 k개를 찾고, 마찬가지로 배치 2의 임의의 세포 j (그림 6B, 빨간점)에 대해 배치 1에서 최근접 이웃 k개를 찾는 작업을 수행한다. 이때 세포 i와 세포 j가 서로의 최근접 이웃 k개에 속하면 두 세포를 한 쌍으로 고정(anchoring)한다. 이렇게 고정된 세포 쌍은 같은 세포 상태에 있는 세포로써 이들 사이의 발현 값 차이는 배치 효과에 의한 것이라고 가정한다. 이 가정에 기반해 비참조(non-reference) 배치의 세포 분포를 참조(reference) 배치의 분포에 맞춰주게 되며, 이를 위해 각 비참조 배치 세포에 대해 인근의 고정된 세포 쌍들의 발현 차이를 이용하여 그 가중 합으로 발현 값을 보정하는 작업을 수행한다(그림 6B, 오른쪽). MNN은 데이터의 원래 차원에서 이런 탐색을 수행하지만, fastMNN과 Seurat의 ‘FindIntegrationAnchors’ 함수는 각각 PCA와 CCA로 축소된 공간에서 수행한다. MNN에 기반한 방식은 종종 계산 부하(속도 및 최대 메모리 사용)가 매우 커서 수십만 개 이상의 세포를 포함하는 데이터셋에는 적용이 제한될 수 있다[86]. 이런 문제를 해결하고자 Scanorama는 최대 메모리 사용량을 줄이도록 고안되었으며, 이로써 scalability를 개선하였다[72, 75, 86]. BBKNN과 Conos는 명시적으로 배치 보정을 수행하지 않고 보정된 발현 행렬 대신 가중 그래프를 제공한다. 이들은 기본적으로 각 세포에 대해 배치 내부에서와 배치 간에 거리가 가까운 세포들을 찾아 연결하고, 축소된 차원에서의 거리가 가까울수록 연결의 가중치를 높게 부여하여 그래프를 만든다. Conos의 경우 배치 내부의 연결에 대한 가중치를 조정하여 배치 내부의 연결을 약화시킨다. BBKNN과 Conos를 통해 얻은 가중 그래프는 세포 군집화와 같은 이후 분석에 사용될 수 있지만, 기능적 유전자 프로그램 식별[89]이나 trajectory 분석[90]과 같은 보정된 발현 행렬을 요구하는 일부 하향 분석에는 적용이 제한될 수 있다[72].

​

b. 세포 군집 수준 유사성 탐색

위의 방법과 달리, 같은 세포 군집에 속하는 비슷한 RNA 발현 프로파일을 갖는 세포들이 배치 요인에 의해 유사한 변동을 겪을 것이라는 가정하에 세포 군집 수준에서 배치 효과를 교정해 줄 수 있다. 이 개념은 Harmony, DESC, LIGER 및 scMerge에서 사용되었다. Harmony는 PCA로 축소된 차원에서 1) 변형된 소프트 k-means 클러스터링(그림 7A, maximum diversity clustering)과 2) 세포 군집 특이적 선형 혼합 모델을 이용한 배치 효과 제거 작업(그림 7B, linear mixture model correction)을 반복한다. 이때 클러스터링 알고리즘은 배치와 세포 군집 할당 간의 독립성을 가정하여 다양한 배치에서 유래된 세포들을 포함하는 세포 군집을 선호하도록 하며, 세포들을 세포 군집에 확률적으로 할당한다. 클러스터링을 통해 세포 군집이 정의되면, 군집별로 배치 종속적인 부분을 선형 모델로 모델링하여 제거한다. 이렇게 보정된 투영 값에 대해 다시 클러스터링과 배치 효과 모델링 및 제거 작업을 수행하며, 세포들의 투영 값이 수렴할 때까지 이 과정을 반복한다. DESC 역시 세포 군집 수준에서 배치 효과를 교정하지만, 일반적으로 세포 종류 간의 변동이 배치 효과보다 크다는 점에 착안하여 배치 정보를 직접적으로 이용하지 않고 세포 군집 내의 순도를 높이는 방식으로 배치 효과를 교정할 수 있도록 설계되었다. 이를 위해 먼저 stacked 오토 인코더를 사용하여 세포들의 비선형 투영을 얻고(그림 7B, 왼쪽), Louvain 클러스터링을 수행하여 클러스터의 수와 중심의 초기 값을 얻는다(그림 7B, 중간). 이후 각 세포 군집 내의 세포들에 대해 중심 값에서 가까운 세포일수록 높은 값을 갖는 metric을 이용하여 보조 확률 분포를 만들고, 이 보조 확률분포와의 KL divergence를 최소화하도록 인코더의 가중치를 조정하는 과정에서 군집 내의 순도가 높아진다(그림 7B, 오른쪽). LIGER 역시 세포 군집 수준에서 배치 보정을 수행한다. LIGER는 차원 축소 방법으로 NMF의 변형인 iNMF를 사용하는데, NMF는 기본적으로 차원 축소와 동시에 클러스터링의 의미를 가지기 때문에 각 세포를 특정 iNMF 요인에 할당할 수 있다. LIGER는 이 할당된 정보를 활용하여 각 세포에서 가장 가까운 k개의 세포가 어떤 요인에 할당되었는지 그 빈도를 참고함으로써 할당된 요인의 정확도를 높이는 전략을 사용한다. 이를 바탕으로 세포를 더 정확하게 군집화한다. 이후 세포 군집 사이의 분포를 quantile normalization[91]을 통해 맞춰준다. LIGER는 생물학적 신호의 보존보다는 배치 효과 교정에 중점을 두는 경향이 있어 교차 종 데이터 셋 통합에 더 유용한 것으로 나타났다[72]. 마지막으로, scMerge는 먼저 각 배치의 세포에 대해 k-means 클러스터링을 수행하고, 상호 간에 가장 가까운 군집 쌍(mutual nearest cluster, MNC)을 찾아 각 군집을 node로 하고 MNC를 edge로 하는 그래프를 만든다. 이 그래프를 서브 그래프로 나누고 각 서브 그래프에 속하는 모든 군집에 대해 군집의 중심에서 가까운 세포들을 핵심 세포로 정의한 뒤, 핵심 세포 사이의 변동을 배치 효과로 간주하여 이 세포들을 pseudoreplicate으로 정의한다. 이후 pseudoreplicate들과, 샘플 사이에 발현 변동이 적은 유전자들(stably expressed genes)을 활용하여 전체 세포에 대한 배치 효과를 추정하는 방식을 사용한다[92].

​

​

 

​

③ 변이형 오토인코더를 사용한 배치 효과 교정 방법

​

유사성 기반 방법은 종종 수십만 개의 세포를 통합할 때 유사성 탐색과 배치 교정을 위한 계산 부하가 매우 크다는 문제가 있다. scVI[93]는 데이터의 비선형적 특성을 효과적으로 모델링하고 계산 부하 문제를 해결하기 위해 개발된 변이형 오토인코더에 기반한 방식이다. scVI는 각 세포의 유전자에 대한 발현 카운트(xm)가 ZINB 분포를 따른다고 가정하며, 'variational posterior'와 'generative model'의 두 부분으로 구성된 모델을 사용한다(그림 8). 첫 번째 부분은 각 세포의 전사체 xm에 대해 정규 분포를 따르는 저차원 잠재 변수(zm)로의 변환을 인코딩한다. 이때 zm은 캡처 효율 및 시퀀싱 깊이로 인한 변동을 나타내는 세포별 크기 요인(lm)과, 배치 요인를 담은 변수(sm)의 존재로 인해 배치 효과가 보정된 투영으로써 세포 군집화와 시각화에 사용될 수 있다. 이후 모델의 두 번째 부분은 비선형 변환을 통해 잠재 변수들을 디코딩하여 각 세포의 각 유전자의 분포 파라미터들을 사후 추정한다. scANVI[94]는 scVI를 확장한 모델로, 세포 유형에 대한 추가적인 정보를 활용하여 잠재 변수 zm의 분포를 추론할 때 세포 유형을 반영하도록 함으로써 배치 보정의 정확도를 높이고자 한 방법이다. 변이형 오토인코더를 사용한 다른 배치 교정 방법으로는 scGen[95]과 trVAE[96]가 있다.

​

 

​

4. 단일세포 오믹스 데이터 기반 궤적 추론

​

궤적 추론(trajectory inference)은 단일세포 오믹스 데이터에 담긴 세포주기, 세포분화, 또는 세포 활성화와 같은 세포의 동적 과정에 대한 정보를 효율적으로 모델링하고 해석할 수 있게 하는 계산 방법론이다[97]. 각 세포는 동적 과정의 여러 단계를 나타내며, 이는 단일 세포들의 분자 정보(mRNA 발현량, 후성유전학적 변형, 염색질 접근성 등) 패턴의 이질성으로 나타나는데, 궤적 추론 방법론들은 세포 간 분자 정보 패턴의 유사성에 기반하여 세포를 정렬하고 다양한 형태의 궤적으로 모델링한 다음, 모델링된 궤적 위에서 동적 과정의 가상의 타임라인을 반영하도록 계산되는 의사 시간(pseudotime)을 각 세포에 부여한다. 이때 사용되는 분자 정보는 주로 mRNA 발현량을 사용하며, 이러한 mRNA 발현량 패턴 변화를 일으키는 전사인자(transcription factor)를 동정하기 위해 동일 세포에서 생산된 단일세포 염색질 접근성 데이터가 통합되기도 한다. 따라서 단일세포 전사체 기반 궤적 추론 방법들을 설명하였다. 궤적 추론은 일반적으로 데이터 정규화, 변동성 높은 변수의 선택, 차원 축소, 궤적 모델링, 의사 시간 추정, 차등 발현 분석의 단계로 진행된다(그림 9A). 데이터 정규화와 변수 선택, 차원 축소의 단계들은 상기된 방법들을 사용하며, 이 단계들은 효율적인 궤적 추론을 위한 데이터 가공 단계이다. 데이터 가공이 끝나면 궤적 추론 방법론들은 먼저 데이터가 대변하는 동적 과정의 궤적을 모델링하는데, 이때 이 궤적은 여러 형태(선형, 갈라지는 형, 트리형, 순환형 등)를 가질 수 있다(그림 9B). 다음으로 이들은 모델링된 궤적을 기반으로 각 세포에 그 동적 과정 내에서의 가상의 순서에 따라 의사 시간을 부여한다(그림 9C). 의사 시간은 모델링된 궤적 내의 여러 경로에 따라 각각 부여될 수 있다. 일반적으로 세포의 동적 과정의 시작점이 되는 세포에 0, 마지막 시점의 세포에 1이 부여된다. 부여된 의사 시간은 특정 경로에서 의사 시간에 따라 특이적이고 유사한 발현 패턴을 갖는 변수(e.g. 유전자)의 군집을 찾을 수 있게 한다(그림 1D). 예를 들어, 유전자 클러스터 3과 4는 각각 경로 1과 2가 분기된 후 발현량이 높아지므로 분기 후 각각의 세포의 상태를 반영하는 유전자 군집임을 알 수 있다. 이러한 유전자 군집들의 생체 경로 농축 분석(pathway enrichment analysis)을 통해 세포의 동적 과정의 특정 시점에서 작용한 유전자들과 관련된 생체경로를 알 수 있다(그림 9E).

 

현재까지 수많은 궤적 추론 방법론들이 고안되었고, 이들은 차원 축소 방법, 세포 군집화 방법, 궤적 모델링 알고리즘, 요구하는 사전 지식의 종류 등에 있어서 다양한 개별적 특성을 가진다. 개별적 특성 중에서 가장 두드러지는 것은 이들이 추론할 수 있는 궤적 위상(topology)의 종류와, 그 궤적 위상의 고정 여부이다. 초창기 제안된 Wanderlust[98]나 waterfall[99]과 같은 방법론들은 선형의 궤적만을 고정하여 모델링 할 수 있도록 고안되었다. 그러나 뒤이어 보고된 Slingshot[100], TSCAN [101]과 같은 방법론들은 다양한 가지 구조를 갖는 트리 형태의 궤적을 데이터에 따라 자유롭게 모델링할 수 있었고, 이후 Monocle3[102], PAGA[103]와 같이 연결되지 않은 그래프와 순환형 구조까지 모델링할 수 있는 방법론들이 제안되었다. 본 보고서에서는 Saelens et al.[104]의 위상 분류 방법을 참고하여 몇 가지 방법론들을 1)궤적 위상을 고정하여 모델링, 2)가지 친 트리 형태까지의 궤적을 자유롭게 모델링, 3)연결되지 않은 그래프 형태까지의 궤적을 자유롭게 모델링하는 방법론으로 나눈다.

​

 

​

가. 궤적 위상을 고정하여 모델링

​

Wanderlust와 Waterfall은 데이터의 선형 궤적을 가정하고 모델링한다. Wanderlust는 가장 먼저 고안된 궤적 추론 방법론 중 하나로서 데이터의 차원 축소나 군집화를 사용하지 않고, 모든 세포로 이루어진 k-nearest neighbor graph를 유전자 발현 변수 공간(feature space)에서 구축하여 궤적을 모델링한다. 모델링된 궤적을 바탕으로 사용자가 지정한 시작 세포(starting cell)로부터의 거리를 기준으로 세포의 위치를 설정하고, 이를 기반으로 의사 시간을 부여한다. 다음으로, waypoints라고 정의된 무작위로 선택된 세포들로부터의 거리 정보를 가중 합계하여 모든 세포들에 부여된 의사 시간을 수정(refine)한다. Waterfall은 PCA를 사용하여 데이터의 차원을 축소한 후, 저차원 변수 공간에서 k-means 클러스터링을 사용하여 세포를 군집화하고, 군집의 중심점들로 이루어진 MST (minimum-spanning tree, 최소신장트리)를 구축한다. 이어서 각 세포를 MST에 수직으로 투영한 뒤, 사용자가 지정한 시작 세포로부터 가장 가까운 세포를 찾아 그 거리로 해당 세포의 의사 시간을 계산하고, 그 세포로부터 다시 가까운 세포까지의 거리를 재귀적으로 더하는 방식으로 모든 세포에 의사 시간을 부여한다.

Wishbone[46]과 DPT[45]는 분기하는 형태의 궤적을 모델링한다. Wishbone은 Wanderlust에서 발전된 방법론으로, 먼저 diffusion maps[105] 방식으로 데이터의 차원을 축소한다. Diffusion maps는 Haghverdi et al.[105]이 기존의 diffusion maps[106] 방법을 단일세포 데이터에 맞게 개발한 차원 축소 방법으로, 고차원 상에서의 데이터의 분포 형태와 동적 변화 궤적을 반영하여 궤적 추론에 적합한 저차원 변수인 diffusion components를 계산한다. Wishbone은 diffusion maps 방법을 통해 short circuit problem(동적 과정상에서는 멀리 위치하나 데이터상에서 가깝게 존재하는 세포들이 그래프상에서 이어지는 오류)을 최소화한다. 다음으로 Wishbone은 Wanderlust와 유사한 방식으로 k-nearest neighbor graph와 무작위로 선택된 waypoints 세포들을 사용하여 궤적을 모델링하고 의사 시간을 부여한다. 추가로, Wishbone은 두 갈래로 갈라지는 형태(bifurcating)의 궤적을 가정하므로 분기하는 지점(branching point)에 대한 탐색이 추가로 필요하다. 무작위로 선택된 두 개의 waypoints 세포를 각각 w1, w2라 할 때, 시작 세포로부터 w2까지의 그래프상에서의 최소 거리 a와, 시작 세포로부터 w1까지의 거리와 w1으로부터 w2까지의 거리를 합한 거리 b와의 차이 (b-a)를 계산한다. 두 waypoints 세포들이 한 가지 또는 줄기 내에 존재한다면 이 차이는 0에 가까운 값을 갖지만, 두 세포들이 서로 다른 가지에 존재한다면 이 차이는 커지게 된다. Wishbone은 모든 waypoints 세포 간의 이러한 차이를 기록한 disagreement matrix Q를 계산하고, 이 행렬의 두 번째 고유벡터를 계산한다. 이 두 번째 고유벡터는 모든 waypoints들의 궤적 내 위치(궤적의 줄기, 분기된 두 개의 가지)에 대한 정보를 담게 된다. DPT는 diffusion maps 방법론과 유사한 방식으로 데이터를 가공하고, 가공된 데이터상에서 세포 간의 유클리드 거리로 diffusion pseudotime을 정의한다. DPT는 이 정의를 바탕으로 지정된 root cell(시작 세포)로부터의 diffusion pseudotime으로 각 세포의 의사 시간을 계산한다. 가지의 분기점은 시작 세포와 그로부터 가장 멀리 떨어진 세포로부터 한 세포씩 의사 시간을 정렬할 때, 의사 시간의 크기가 한쪽에서는 커지고 한쪽에서는 작아지지만 분기점에 도달하여 새로운 가지에 존재하는 세포에 도달하게 되면 둘 모두에서 의사 시간이 커지게 되고 이 지점을 분기점으로 결정한다.

​

나. 가지 친 트리 위상의 궤적을 모델링

​

Waterfall은 MST 알고리즘을 통해 가지 친 형태의 궤적을 모델링하지만, 선형의 줄기에 모든 세포를 투영하여 의사 시간을 부여한다. TSCAN과 Slingshot은 모두 MST 알고리즘을 통해 가지 친 트리 형태의 궤적을 데이터에 따라 자유롭게 모델링하고, 가지 친 경로들에 독립적으로 의사 시간을 부여한다. TSCAN은 유전자들을 hierarchical clustering으로 군집화한 다음, 군집 별로 유전자 발현을 평균하여 군집 단위의 발현 데이터로 변환시킨다. PCA를 사용하여 변환된 데이터의 차원을 축소한 후, 이 공간 내에서 Mclust 알고리즘으로 세포를 군집화한다. 다음 과정은 Waterfall과 유사한 방식으로 군집 중심점으로 MST를 구축하고, 의사 시간을 부여한다. Slingshot은 먼저 다양한 방법(PCA, ICA, diffusion maps 등) 중 하나를 선택하여 데이터의 차원을 축소할 것을 요구한다. 다음으로 Slingshot은 Waterfall과 유사하게 임의의 방법으로 세포를 군집화한 다음, 군집의 중심점들로 MST를 구축하여 전체적인 궤적을 모델링한다. 다음으로 Slingshot은 principal curve[107] 알고리즘을 변형한 simultaneous principal curve라는 고유한 방법론을 사용하여 MST로 모델링된 궤적을 부드러운 형태로 다듬어 궤적을 최적화하고, 세포 군집화 결과나 노이즈 데이터로부터 분석을 견고하게(robust) 한다. 모든 세포는 이렇게 만들어진 궤적 곡선 위에 수직으로 투영되고 의사 시간이 계산된다. TSCAN과 Slingshot에서 사용자는 의사 시간 계산을 위해 시작 또는 마지막 군집을 지정할 수 있고, Slingshot의 경우 기존의 데이터에 대한 지식에 따라 마지막 군집(terminal cluster)의 개수를 지정하여 분기점(branching point)의 탐색을 용이하게 할 수 있다.

​

다. 연결되지 않는 그래프 위상의 궤적을 모델링

​

단일세포 데이터는 때로 다수의 다른 세포와는 독립적인 동적 과정을 갖는 일부 세포 군집을 포함한다. 이러한 세포들의 궤적을 정확하고 독립적으로 모델링하기 위하여, PAGA와 Monocle 3는 순환형 고리와 여러 개의 연결되지 않은 그래프 형태의 궤적까지 모델링할 수 있는 방법으로 제안되었다. PAGA는 PCA 등의 방법으로 데이터의 차원을 한 번 축소하고, 선택 사항으로 diffusion maps를 이용하여 데이터의 차원을 한 번 더 축소한다. PAGA는 저차원 변수 공간에서 모든 세포로 이루어진 k-nearest neighbor graph를 구축하고, Louvain clustering 등의 알고리즘을 이용해 세포를 군집화한다. 다음으로 PAGA는 군집 사이의 유사성 척도를 도입하여 군집을 노드로 하는 그래프를 구축한다. 이러한 작업을 통해 PAGA는 신뢰도 낮은 세포 간의 edge를 폐기할 수 있고, 노이즈가 제거되고 연결되지 않은 부분까지 잘 나타낸 궤적의 형태를 제시할 수 있다. 구축된 그래프 위에서, PAGA는 DPT와 유사한 방식으로 diffusion maps로 계산된 변수들과 시작 세포를 사용하여 의사 시간을 부여한다. Monocle 3는 PCA로 데이터의 차원을 축소한 후, 선택 사항으로 t-SNE[108] 또는 UMAP[109] 방법으로 데이터의 차원을 한 번 더 축소한다. 이어서 Monocle 3는 PAGA에서와 똑같은 방법으로 세포들의 그래프를 구축하고, Louvain 알고리즘으로 세포를 군집화한다. 이어서, 트리 형태의 그래프를 구축하는 DDRTree[110], DDRTree와 유사하나 차원 축소를 동반하지 않는 SimplePPT[111], 순환 구조를 포함할 수 있는 L1Graph[112] 중 하나의 알고리즘을 선택하여 군집의 중심점들을 노드로 하는 그래프를 구축하고, 모든 세포를 그래프 위로 투영한 후 사용자가 지정한 시작 세포로부터의 그래프상의 거리로 의사 시간을 계산한다. 이렇게 추론된 궤적을 따라 발현량이 변하는 유전자군을 동정하고, 동일 세포에서 생산된 염색질 접근성 데이터를 통합하여 궤적 초기, 중기, 말기에 유전자 발현량을 높이는 전사인자를 동정한다. 특히 궤적 중기 특이적으로 발현량이 증가했다가 말기에 감소하는 유전자군의 발현을 조절하는 전사인자들은 해당 궤적이 대변하는 세포 분화의 스위치로 작용할 수 있으므로, 이들을 활성화하는 리간드-수용체를 동정하여 해당 세포군의 분화를 조절할 수 있는 치료타겟으로 제시하기도 한다.

​

5. 맺음말

​

15종의 암종에 대한 유전단백체(Proteogenomics) 분석들을 통해, bulk 다중오믹스 분석이 단일오믹스에 비해 질환 서브타이핑, 서브타입 특성 해석 등에서 가지는 장점들이 보고되어 왔다. 이러한 다중오믹스 분석을 통해 증진된 질환의 특성 분석을 통해, 진단 및 치료를 위한 다양한 분자 타겟들이 제시되어왔다. 단일세포 수준에서도 다중오믹스 분석을 통해, 보다 증진된 질환 관련 세포군 동정과 그들의 특성 해석이 가능하다. 본 보고서에서 정리한 단일세포 다중오믹스 실험 방법 및 데이터 분석 방법을 적용하여 단일세포 다중오믹스 데이터의 통합 분석을 한다면, 주요 질환 표현형(질환 발생‧이행, 약물 반응성 등) 관련 세포군을 동정하고, 이들 세포군에서 특이적인 발현, 변형을 보이는 분자 정보(mRNA‧단백질 발현량, 후성유전학적 변형‧염색질 접근성 등)를 동정하며, 이들 분자 정보를 기반으로 해당 세포군의 특성을 규명하고, 통합된 데이터를 기반으로 이들 특성을 조절하는 전사인자, 리간드-수용체를 동정할 수 있다. 더 나아가 이들을 네트워크 모델로 통합하여 질환 관련 표현형을 해석 등을 수행할 수 있다. 이러한 연구를 통해 구축된 네트워크 모델을 분석하면, 주요 질환 표현형을 진단할 수 있는 고신뢰성 진단마커, 또는 표현형을 조절할 수 있는 새롭고, 정확한 세포군과 치료 타겟들을 제시할 수 있을 것으로 기대한다.

 

 

출처 : 바이오인