Introduction
최근 single-cell RNA-seq은 약물 처리가 세포 하나하나의 유전자 발현을 어떻게 바꾸는지 볼 수 있게 만들었다. 이 덕분에 기존 bulk 실험에서 놓치기 쉬웠던 세포 간 이질성과 미세한 반응 차이를 포착할 수 있었다. 하지만 비용과 규모의 문제가 존재한다. single-cell HTS는 정밀하지만 수천 개 약물을 실험하기에는 너무 비싸고, 실제 데이터셋도 대부분 수백 개 이하의 약물만 포함한다.
그래서 모든 약물을 single cell로 직접 실험할 수 없다면, 모델이 관측하지 않은 perturbation을 예측해야 한다고 주장한다. 특히 중요한 목표는 학습 중 보지 못한 새로운 약물의 반응을 예측하는 것으로 이는 drug repurposing이나 후보 약물 발굴에 직접적으로 연결되지만, 현재 single-cell 데이터만으로는 다양한 화학 구조와 세포 맥락을 충분히 학습하기 어렵다.
chemCPA는 이러한 한계를 해결하기 위해 제안되었으며, 기존 CPA처럼 세포의 basal state, drug effect, covariate effect를 latent space에서 조합하되, 약물을 단순 ID embedding으로 보지 않고 molecular structure에서 drug embedding을 만든다. 그래서 학습 데이터에 없던 약물이라도 분자 구조 정보가 있으면 perturbation effect를 예측할 수 있다.
[논문 리뷰] Compositional perturbation autoencoder for single-cell response modeling
perturbation이 세포의 유전자 발현에 어떤 영향을 미치는지를 정밀하게 분석하는 것은 실제 실험 환경에서 조합이 매우 다양하므로 불가능하다. 따라서 관측되지 않은 새로운 perturbation조건에서의
velog.io
또한 bulk RNA HTS 데이터를 활용한다. 비록 bulk 데이터는 single cell 해상도는 없지만 훨씬 많은 약물을 포함하므로, 약물 구조와 유전자 발현 변화의 관계를 먼저 배우기에 적합하다. 즉, 제한된 single-cell 데이터를 대규모 bulk perturbation 데이터로 보완해, unseen drug response를 single-cell 수준에서 예측한다.
Chemical Compositional Perturbation Autoencoder
논문은 데이터셋을 다음과 같이 정의한다.
$$\mathcal{D} = \{(x_i,y_i)\}^N_{i=1} = \{(x_i,(d_i,s_i,c_i))\}^N_{i=1}$$
- $x_i$ : $n$차원의 gene expression
- $y_i$ : 속성 집합
- $d_i$ : 약물에 대한 정보
- $s_i$ : 약물의 용량에 대한 정보
- $c_i$ : cell line에 대한 정보
chemCPA가 풀고 싶은 문제는 단순한 reconstruction이 아니라 관측하지 않은 조건의 유전자 발현을 예측하는 것이다. 더 나아가 chemCPA는 학습 중 아예 등장하지 않은 새로운 약물에 대해서도, 그 약물의 분자 구조를 이용해 세포 반응을 예측하려 한다. 이런 예측을 논문에서는 counterfactual prediction이라고 부른다.
이러한 counterfactual prediction을 위한 한 가지 가능한 접근법은, 세포의 gene expression $x_i$을 속성$y_i$과 무관한 latent vector $z_i$(basal state)로 인코딩한다. 즉, basal state 안에는 drug 정보나 cell line 정보가 최대한 없어야 한다. 이러한 basal state가 주어지면 $z_i$는 속성 표현 $z_d$ 및 $z_c$와 결합되어 임의의 속성 조합 $y'_i$를 인코딩할 수 있으며, 이후 새로운 속성 집합이 반영된 유전자 발현 상태 $\hat{x}_i$로 다시 디코딩될 수 있다.
이를 위해서 chemCPA를 세 부분으로 세분화한다.
1. Encoder-Decoder
Encoder($E_{\theta}$)는 파라미터 $\theta$를 가지는 MLP로 $n$차원 gene expression을 $l$차원 latent vector $z_i$로 바꾼다. 이 과정에서 adversairal classifiers를 통해 $z_i$에 속성 $y_i$에 대한 정보가 포함되지 않도록 학습된다. 이렇게 생성된 $z_i$에 우리가 원하는 속성 임베딩을 더해 $z'_i = E_{\theta}(x_i) + z_{attribute}$를 얻는다
Decoder($D_{\psi}$) 또한 MLP로 $z'_i$를 입력으로 받아 gene expression 기반 분포 $\mathbb{P}$에 대한 성분별 파라미터를 계산한다. 왜냐하면 gene expression 데이터는 noise가 많아, 같은 조건의 세포라도 expression 값이 완전히 같지 않으므로 gene expression 값을 정확하게 예측하는 것보다, 평균과 분산으로 각 gene의 expression이 따를 확률분포를 예측하는 것이 자연스럽기 때문이다.
평균과 분산으로 파라미터화한다고 가정하면, 두 경우 모두 디코딩된 gene expression 상태를 설명하기 위해
$$\mu = D^{\mu}_{psi}(z') \ 및 \simga^2 = D^{sigma^2}_{\psi}(z')$$
을 얻는다.
chemCPA는 Gaussian likelihood에서 더 나은 수렴을 관찰하였고 이 경우 reconstruction loss는 다음과 같다
$\mathcal{L}_rec(\theta,\psi) = N(x_i \mid \mu_i,\sigma_i) = \frac{1}{2} \left[ln(D^{\sigma^2}_{\psi}(z'_i) + \frac{(D^{\mu}_{\psi}(z'_i) - x_i)^2}{D^{\sigma^2}_{\psi}(z'_i)}\right]$
- $ ln(D^{\sigma^2}_{\psi}(z'_i)$ : 분산을 너무 크게 잡지 못한게 하는 penalty
- $ \frac{(D^{\mu}_{\psi}(z'_i) - x_i)^2}{D^{\sigma^2}_{\psi}(z'_i)} $ : 예측 평균이 실제값과 가까워지도록 만드는 penalty
2. Attribute embedding and additive latent space
latent space에서 perturbation 반응이 additive 구조를 가진다고 가정한다
$$z'_i = z_i + z_{attribute} = z_i + z_{c_i} + \hat{s}_i z_{d_i}$$
이러한 additive한 구조는 사용자가 모델을 해석할 수 있게 하며, 순열 불변성을 가지며, 미세조정 과정에서 새로운 covariate를 추가할 수 있다는 장점을 가진다.
drug 속성과 cell line 속성은 성격이 서로 다르므로, 이들을 latent space에서 별도로 인코딩한다. cell line은 CPA와 동일한 접근법을 사용하며, drug의 경우, 새로운 임베딩 네트워크 $P_{\mathcal{\phi}$를 제안한다.
Perturbation Network
$P_ φ $는 분자의 그래프나 SMILES 같은 분자 표현, 사용된 용량을 latent perturbation state로 매핑한다. 이 네트워크는 분자 인코더 $G$, perturbation 인코더 $M$, 용량 스케일러 $S$로 구성된다
$G$는 약물의 분자 구조를 고정된 크기의 임베딩 $h_d$으로 인코딩하고, $M$은 분자 임베딩 $h_d$를 입력으로 받아 drug perturbation vector $z_d$를 생성한다. $S$는 $h_d$와 용량 $s_i$를 받아 스케일링된 용량 값 $\hat{s}_i$를 계산한다. 여기서 $z_d$는 약물의 일반적인 효과 방향, $\hat{s}$는 용량에 따른 효과의 세기를 나타낸다.
$$\hat{s}_i \times z_{d_i} = P_ φ (g_i,s_i) = S(h_{d_i},s_i) \times M(h_{d_i})\ with \ h_{d_i} = G(g_i)$$
$G$는 SMILES, molecular graph, RDKit feature 등 어떤 분자 표현을 사용해도 되지만, scRNA-seq HTS에서 이용 가능한 약물 수가 제한적이므로, 사전학습된 인코딩 모델에 의존하고 학습 중에는 $G$를 고정한다. 특히 RDKit을 분자 인코더 $G$로 사용한 경우가 좋은 성능을 보여주었다.
3. Adversarial classifiers for invariant basal states
basal state를 생성하고 $z_i,z_{d_i},z_{c_i}$의 분리된 표현을 만들기 위해서 adversarial classifiers $A^{drug}_{\phi}, A^{cov}_{\phi}$를 사용한다. 두 adversary network는 $z_i$를 입력으로 받아, 예시 $i$에 적용된 dru$d_i$g와 cell line $c_i$를 예측하는 것을 목표로 한다.
이 classifiers는 분류 성능을 향상시키도록 학습되는 반면, $E_{\theta}$의 학습 목적 함수에는 classification loss를 부호를 반대로하여 추가한다. 따라서 인코더는 속성에 대한 정보를 포함하지 않는 latent representation $z_i$를 생성하려 한다.
두 분류기 모두에 대해 Cross-entropy를 사용한다
$$\mathcal{L}^{drugs}_{class} = CE(A^{drug}_{\phi}(z_i),d_i) \ and \ \mathcal{L}^{cov}_{class} = CE(A^{cov}_{\phi}(z_i),c_i)$$
CPA 구현을 따라, adversarial classifiers의 손실 함수에 gradient penalty를 추가하여 다음 항을 최소화 한다
$$\mathcal{L}^j_{pen} = \frac{1}{k}\sum_k \left\| \partial_{z_i} A_{\psi}^{j}(z_i)_k \right\|_2^2$$
이 gradient penalty를 discriminator를 nosie에 강건하게 만들고 localconvergence를 가능하게 하는 것으로 나타났다.
학습 중에는 다음 경쟁적 목적함수들 사이에서 업데이트 단계를 번갈아 수행한다
$$\mathcal{L}_{AE}(\theta,\psi, φ \mid \phi) = \mathcal{L}_{rec}(\theta, \psi, φ) - \lambda_{dis} \sum_j \mathcal{L}^j_{class}(\theta \mid \phi) \ and \ \mathcal{L}_{Adv}(\phi \mid \theta) = \sum_j \mathcal{L}^j_{class}(\phi \mid \theta) + \lambda_{pen} \mathcal{L}^j_{pen}(\phi)$$
Datasets and Transfer learning
- L1000
- data type : bulk RNA
- 관측치 수 : 약 130만 개
- gene 수 : 978개
- drug 수 : 약 2만개
- 특징 : FDA 승인 약물도 있고, 합성 화합물도 있음
- drug 수가 매우 많으므로 다양한 약물 구조와 gene expression 변화의 관계를 미리 배울 수 있음
- sci-Plex3
- data type : single-cell RNA-seq
- cell 수 : 649,340
- gene 수 : 7561
- cell lien : A549, MCF7, K562
- drug 수 : 188 개
- dosage : 10nM, 100nM, $1_{\mu}M, 10_{\mu}M$
- preturbation : single-compund perturbation
- 세포 하나하나의 반응을 볼 수 있음
- Transfer learning
- Gaussian likelihood loss로 gene expression count 값을 모델이 다루기 좋은 연속적인 값으로 바꿈
- 두 데이터셋의 gene set을 같게 맞춤, L1000, sci-Plex3에서 공통으로 대응되는 gene만 골라서 둘다 977개 gene으로 맞춤
- pretraining / finetuning의 input / output 차원이 같으므로 transfer learning이 쉬움
- L100 977개 gene 만으로는 single-cell 데이터의 세포별 다양성 설명이 어려움, sci-Plex3에 HVGs 1023개 추가
- input이 977, output이 2000으로 바뀜
- 두 개의 layer를 추가하여 문제를 해결함
- $h_{enc} : R^n \rightarrow R^n$ : 인코더 앞에 두어 2000 차원을 977차원으로 축소
- $h_{dec} : R^{2n} \rightarrow R^{2n}$ : 디코더 뒤에 두어 977 차원을 2000차원으로 변환
- 이때 $2n$인 이유는 디코더는 Gaussian likelihood를 위해 gene 마다 $\mu, \sigma^2$을 출력함
Experiments
Comparing chemCPA against existing methods on unseen drug-covariate combinations
학습 때 drug 와 cell line은 본적이 있지만 특정 drug-cell line 조합은 보지 못한 상황에서 예측할 수 있는가에 대한 실험을 진행하였다.
기존 scGen과 CPA는 학습때 보지 못한 drug를 처리하기 어려우므로 unseen drug-covariate combination을 사용하였다. 또한 scGen은 dosage를 명시적으로 구분해서 처리하지 못하므로 공정한 비교를 위해 dose 별로( 1 µM와 10 µM ) 실험을 따로 나눠서 실험하였다.
세 개의 cell line을 하나씩 돌아가며 test로 사용하였으며, 총 9가지 화합물을 테스트 대상으로 선택하였다. 이 약물들 대부분 후성유전학적 조절, tyrosine kinase ignaling, 세포주기 조절에 속한다.
성능은 실제 gene expression과 모델의 counterfactual prediction 사이의 r²로 평가했다. 전체 gene 기준 점수와, 약물 효과가 강하게 나타나는 DEGs 기준 점수를 모두 보고했다. baseline은 perturbation 정보를 사용하지 않는 모델이므로, baseline보다 성능이 높다는 것은 drug encoding이 실제로 예측에 기여했다는 의미다.
결과적으로 chemCPA는 scGen과 CPA보다 좋은 성능을 보였고, 특히 L1000으로 pretraining한 chemCPA가 가장 우수했다. 흥미롭게도 pretraining하지 않은 chemCPA도 CPA보다 좋은 성능을 보였는데, 이는 drug를 단순 ID embedding으로 외우는 대신 molecular representation을 통해 perturbation을 생성하는 구조가 regularization 효과를 주기 때문으로 해석된다.
Using chemCPA to predict single-cell responses for unseen drugs
chemCPA가 학습 중 보지 못한 약물의 single-cell 반응을 예측할 수 있는지 평가한다. 특정 약물을 학습에서 제외한 뒤 그 약물의 반응을 예측한다.
앞에서와 동일하게 9개의 약물을 unseen drug로 사용하고, 두 가지 gene set에서 실험한다
Shared gene sets
L1000과 sci-Plex3에 공통으로 존재하는 977개 shared genes
shared gene set 실험에서는 pretrained chemCPA가 baseline과 non-pretrained chemCPA를 일관되게 앞선다. 낮은 dose에서는 약물 효과가 거의 없어 baseline도 높은 점수를 보이지만, 높은 dose에서는 약물 효과가 뚜렷해지면서 chemCPA의 장점이 나타난다. 특히 DEGs 기준 성능은 전체 gene보다 낮지만, pretrained chemCPA는 약물로 인해 실제로 변한 gene expression도 어느 정도 설명할 수 있음을 보여준다.
Extended gene sets
sci-Plex3의 1023개 HVGs를 추가한 2000개 extended genes
extended gene set에서도 같은 경향이 유지된다. 2000개 gene을 사용하는 더 어려운 설정에서도 pretrained chemCPA는 가장 좋은 성능을 보였고, non-pretrained chemCPA는 baseline보다 약간 나은 수준에 머물렀다. 이는 L1000 같은 대규모 bulk RNA perturbation 데이터에서 배운 정보가 single-cell 데이터로 전이될 수 있으며, 심지어 두 데이터셋의 gene set이 다를 때도 도움이 된다는 점을 보여준다.
Measure uncertainty on the drug embedding
unseen drug 예측을 얼마나 믿을 수 있는지 평가하기 위한 uncertainty score를 제안하는 부분
chemCPA는 보지 못한 약물에 대한 반응을 예측하도록 설계되었지만, 그 일반화 능력은 결국 학습 데이터의 범위에 의해 제한된다. sci-Plex3 데이터에서는 전체 drug-dose 조합 중 20% 미만만이 control phenotype과 r2r^2 기준으로 35% 이상 차이를 보였고, 명확한 효과를 보이는 약물도 tyrosine kinase signaling, DNA damage and repair, cell cycle regulation, epigenetic regulation 등 일부 pathway에 집중되어 있었다. 저자들은 이러한 기술적 노이즈와 제한적인 약물 효과가 사전학습되지 않은 chemCPA가 확장 유전자 집합에서 baseline을 크게 넘어서지 못한 이유라고 해석한다. 반면 L1000 bulk RNA perturbation data로 사전학습한 모델은 더 다양한 약물 반응 정보를 활용할 수 있어 더 강건한 성능을 보였다.
이 한계는 perturbation latent space에서도 확인된다. 이상적으로는 비슷한 MoA를 가진 약물들이 잠재 공간에서 함께 클러스터링되어야 하지만, 실제로는 일부 약물에서 클러스터링이 불완전하게 나타난다. 특히 baseline 점수가 높은 경우, 즉 약물 효과가 control과 크게 다르지 않은 경우에는 chemCPA가 뚜렷한 교란 효과를 식별하기 어렵다.
이를 보완하기 위해 저자들은 약물 임베딩 공간에서의 불확실성 척도를 제안한다. 핵심은 KNN graph를 이용해 특정 약물 주변의 이웃 약물들이 어떤 MoA를 가지는지 확인하는 것이다. 이웃들이 동일한 pathway에 속하면 모델의 확신이 높고, 여러 pathway가 섞여 있으면 불확실성이 높다고 볼 수 있다. 여기에 이웃 약물들과의 거리 정보도 함께 반영하여, 약물이 주변 약물들과 얼마나 구별되는 교란 효과를 갖는지 평가한다.
$$u_i = \sum_{j \in N_i} \frac{1}{log(d(i,j))} \times H(X)$$
- $d(i,j)$ : 약물 간 유클리드 거리
- $H(X)$ : 이웃 약물들의 pathway 분포에 대한 Shannon entropy
결과적으로, unseen drug에 대한 불확실성 점수는 실제 예측 성능과 잘 상관되었다. 이는 chemCPA의 잠재 공간이 단순히 예측을 만드는 데 그치지 않고, 모델이 어떤 약물에 대해 신뢰할 만한 일반화를 수행할 수 있는지 해석하는 데에도 활용될 수 있음을 보여준다
Conclusion
이 논문은 보지 못한 약물(unseen drug) 에 대한 단일세포 유전자 발현 반응을 예측하기 위해 chemCPA를 제안한다. chemCPA는 약물의 분자 구조 정보를 활용해, 학습 데이터에 없는 약물의 교란 효과까지 예측할 수 있도록 설계되었다.
주요 기여는 세 가지다. 첫째, 분자 표현으로부터 단일세포 수준의 약물 반응을 예측하는 chemCPA를 제안했다. 둘째, 대규모 bulk RNA-seq HTS 데이터를 사전학습에 활용하는 전이학습 전략을 도입했다. 셋째, chemCPA가 CPA와 scGen보다 우수한 성능을 보이며, 기존 방법으로는 어려웠던 unseen drug 일반화 과제까지 다룰 수 있음을 보였다.
또한 저자들은 unseen drug 예측의 신뢰도를 평가하기 위한 불확실성 척도를 제안했다. 전체적으로 chemCPA는 단일세포 기반 약물 스크리닝과 신약 발견에 활용될 수 있는 가능성을 보여준다.