Introduction
기존 약물 발굴은 주로 질병 관련 단백질/효소를 먼저 정의하고 해당 타겟을 억제하는 화합물을 설계하는 one drug - one gene - one disease 방식의 표적 기반 약물 발굴이었다. ( imatinib ▶️ BCR-ABL 억제, trastuzumab ▶️ HER2 억제 )
하지만 실제 FDA에 승인된 많은 first-in-class 약물이 특정 약물 표적 가설 없이 발견되었으며, 이는 타겟을 먼저 정하는 접근이 항상 효과적이지 않음을 시사한다. 이는 기존의 표적 기반 약물 발굴이 복잡한 질병 메커니즘을 충분히 설명하지 못함을 시사한다.
이에 대한 대안으로 표현형 기반 약물 발굴이 등장하였다. 사전에 정의된 표적 없이 수행되는 분석을 통해 측정된 질병 표현형을 되돌리는 화합물 또는 치료 표적 조합(perturbagen)을 식별하는 데 초점을 두었다.
하지만 최근 딥러닝 기반 표현형 접근은 다음과 같은 한계가 존재한다.
- 화학적 및 유전적 라이브러리에 의존하므로 사전에 정의된 라이브러리 내에서만 perturbagen을 선택할 수 있고, 새로운 약물 표적 조합으로서의 perturbagen을 생성 할 수 없다.
- 라이브러리 내 모든 perturbation에 대한 표현형 변화만 예측하고, 원하는 반응을 보이는 perturbagen을 탐색하는 방식으로 perturbagen을 식별한다. 하지만 우리가 원하는 것은 표현형을 바꾸기 위해서 어떤 perturbagen이 필요한지를 예측하는 것으로 inverse problem을 해결하지 못한다
Methods
문제 정의를 하면 다음과 같다 $x^d$ (diseased gene expression), $x^t$(treated gene expression)이 입력으로 들어오고$\mathcal{U}'$이 출력된다. 즉, diseased 상태를 treated 상태로 이동시키는 perturbagen 예측 문제이다. (inverse perturbation prediction)
기존 방식 : perturbagen ▶️ 표현형 변화 예측 (forward problem)
PDGrapher : diseased 상태 ▶️ treated 상태 : 어떤 perturbagen이 필요? (inverse problem)
두 모듈로 구성되면 두 모듈은 같은 학습 과정( 병렬로 학습)에서 최적화되지만, perturbation discovery module 학습 시 response predictor는 frozen 상태로 사용한다.
Preturbagen discovery module
$$f_p(x^d,x^t) \rightarrow \mathcal{U}'$$
- 각 노드 구성 : $\left[x^d_i, x^t_i\right]$ 각각 현재 상태(diseased)와 목표(treated) 상태를 의미한다
- 임베딩 구성 : diseased embedding + treated embedding + positional embedding = 3d - dimention
- diseased embedding, treated embedding은 B-bin embedding을 사용하는데, expression 범위를 b개의 bin으로 구간을 나눠서 bin id로 embedding을 한다. (0.0 ~ 0.2는 bin 1, 0.3 ~ 0.4는 bin 2...)
- 같은 gene이면 샘플이 달라도 embedding이 동일함, 서로 다른 gene이면 샘플이 같고 bin이 같아도 다른 임베딩
- GNN Message Passing : $h'_{i} = \phi(h_i \oplus_{j \in \mathcal{N}(i)} \psi(h_i, h_j))$
- K번 message passing 후 최종 node embedding $z_i$를 구함
- $z_i$를 Multilayer Feedforward로 점수 $s_i$를 생성한다.
- 이렇게 생성된 각 노드 별 perturbation score $s_i$를 기반으로 높은 상위 P개의 노들를 선택하여 perturbation 집합 $\mathcal{U}$를 구성한다.
Response predictor module
위 이미지는 Inference 단계 이미지로 $x^d$만 입력으로 들어가지만, 학습 과정에서 $f_r(x^h,\mathcal{U}) \rightarrow x^d$와 $f_r(x^d,\mathcal{U}') \rightarrow x^t$ 두 가지를 모두 사용한다. ➡️ disease를 만드는 perturbation도 배우고, disease를 고치는 perturbation도 배우면서 intervention response function이 더 일반화됨
여기서 $\mathcal{U}$는 disease-association gene set, $\mathcal{U}'$은 chemical/genetic intervention이 적용된 gene set을 의미한다
$$f_r(x^d, \mathcal{U}') = \hat{x}^t$$
이때 예측한 treated expression을 정답 expression과 가까워지도록 최적화 한다.
- 각 노드 구 : $\left[x^h_i, x'_{\mathcal{U}}\right]$ : 정상 상태 gene expression, perturb flag(i번째 gene perturbation 여부)
- 임베딩 구성 : gene expression embedding + perturb embedding + positional embedding ▶️ 3d-dimension
- 이때 gene expression embedding은 B-bin embedding을 사용하는데, expression 범위를 b개의 bin으로 구간을 나눠서 bin id로 embedding을 한다. (0.0 ~ 0.2는 bin 1, 0.3 ~ 0.4는 bin 2...)
- GNN Message Passing : $h'_{i} = \phi(h_i \oplus_{j \in \mathcal{N}(i)} \psi(h_i, h_j))$
- K번 message passing 후 최종 node embedding $z_i$를 구함
- 이 $z_i$를 MLP에 넣어 scalar 값을 예측 $\hat{x}^t_i = MLP_r(z_i)$
- 모든 gene node에 대해 반복하면 전체 예측 treated expression vector가 된다 $\hat{x}^t = \left[\hat{x}^t_1, \hat{x}^t_2,...,\hat{x}^t_n \right]$
Model Optimization
response prediction module $f_r$
$$\mathcal{L}_{f_r} = CE(x^d, f_r(x^h, \mathcal{U})) + CE(x^t, f_r(x^d, \mathcal{U}'))$$
intervention discovery module $f_p$
$$\mathcal{L}_{f_p} = CE(x^t, f_r(x^d,f_p(x^d,x^t))) + CE(\mathcal{U}',f_p(x^d,x^t)) (with\ f_r \ frozen)$$
왜 $f_r$를 freeze할까?
cycle loss 계산 과정에서는 $f_r$가 freeze된 상태로 사용되며, 이는 perturbation discovery module이 response predictor 자체를 수정하는 대신 올바른 perturbation 집합을 찾도록 유도하기 위함이다.
- $ CE(x^t, f_r(x^d,f_p(x^d,x^t))) $ : Cycle loss, $f_p$가 예측한 perturbation 집합을 $f_r$에 입력했을 때 목표 상태 $x^t$가 잘 재현되도록 한다.
- $ CE(\mathcal{U}',f_p(x^d,x^t)) $ : Supervisor signal, diseased 상태 $x^d$와 treated 상태 $x^t$가 주어졌을 때, 실제로 사용된 perturbation set $\mathcal{U}'$를 맞추도록 $f_p$를 지도학습한다.
cycle loss를 통해 예측된 개입 집합 $\mathcal{U}'$에 대한 반응이 목표로 하는 치료 상태 $x^t$와 가깝도록 만들며, $\mathcal{\hat{U}}'$을 예측하도록 지도 신호를 제공하기 위해 cross-entropy 손실을 추가한다.
만약 cycle loss만 있게되면 정답 perturbation 집합을 맞추도록 학습되지 않고, 효과만 맞는 아무 perturbation을 찾도록 학습된다.
Example
여러 perturbation set이 존재한다. $\mathcal{U}'_1, \mathcal{U}'_2, \mathcal{U}'_3,...$
모두 동일한 phenotype 변화 유도 가능 $f_r(x^d, \mathcal{U}'_1) \approx f_r(x^d, \mathcal{U}'_2) \approx x^t$
이 경우 cycle loss만 있으면 모델은 아무거나 선택해도 된다.
정답 : $\mathcal{U} = \{ gene A, gene B\}$
모델 예측 : $\mathcal{U}' = \{geneC, geneD\}$
인데도 cycle loss가 작다. 즉, correct target recovery는 보장되지 않는다.
그래서 supervision loss가 추가되어 정답 perturbation set과 일치하도록 강제한다.
Result
Dataset
PDGrapher는 다양한 intervention type, cancer type, causal graph 조건에서 일반화 성능을 검증하기 위해 총 38개의 dataset 조합에서 평가되었다.
- intervention type : chemical / genetic
- cancer type : 11개(폐암, 유방암, 전립선암, 대장암, 피부암, 자궁경부암, 두경부암, 췌장암, 위암, 뇌종양, 난소암)
- causal graph : PPI / GRN
- PPI network는 BIOGRID database에서 가져왔으며,GRN은 각 dataset에 대해 GENIE3 알고리즘을 사용하여 gene expression 데이터로부터 추정하였다.
- cell line 조합
PDGrapher 학습 데이터는 두 가지 구성 요소로 이루어짐
1. disease intervention dataset
- (healty, diseased, associated disease genes)
- Lung, Breast, Prostate cell line에서만 가능함
2. treatment intervention dataset
- (diseased, treated, perturbagen)
- intervention : CRISPR, compund
모델의 Perturbagen 예측 성능 평가
성능은 nDCG, recall@k, exact-hit ratio, network distance를 기반으로 평가하였다.
모든 cell line과 모든 metric에서 기존 방법보다 높은 성능을 보인다.
- nDCG : 예측된 gene target ranking이 실제 therapeutic target ranking과 얼마나 일치하는가
- Extract-hit ratio : 모델이 예측한 상위 후보 중 실제 정답이 얼마나 포함되어 있는가
- recall@k : $\frac{$top-k$안에 \ 포함된 \ 실제\ target\ 수}{전체\ 실제\ target\ 수}$
- network distance : 예측이 틀렸어도, biological network 안에서 얼마나 가까운 위치의 gene을 맞췄는가를 평가
- x축은 PDGrapher - reference로 음수일 수록 PDGrapher가 예측한게 더 가까운 위치의 gene임을 의미한다.
- y축은 해당 distance 차이를 가진 샘플 비율을 의미한다
- 이 결과를 통해 PDGrapher가 단순히 랜덤하게 gene을 고르는 것이 아니라 실제 therapeutic target이 위치한 network neighborhood를 찾아냄을 보여준다.
Unseen Cell line generalization 검증
한 cell line에서의 random split만 보게 되면 해당 cell line의 패턴을 외울 수도 있다. 그래서 실제로 일반적인 규칙을 학습했는지를 보기 위해서 leave-cell-out ( 한 개 cell line을 제외하고 나머지로 학습한 뒤, 제외한 cell line으로 테스트 ) 을 수행한다.
unseen cell line에서도 기존 방법보다 높은 성능을 유지하여 cell-line specific pattern을 학습하지 않고 generalize 규칙을 학습했음을 보여준다
모델이 예측한 perturbagen이 실제로 생물학적으로 의미 있는 target인지 검증
(1) FDA-approved drug target recovery
training set에 없는 drug target(target gene)을 복원 가능하였고, 이는 generalizable mechanism을 학습했음을 의미한다
이 실험은 PDGrapher가 질병 상태의 세포주를 건강한 상태로 전환하기 위해 표적이 될 유전자들을 예측하고 상위 K개 gene 안에 FDA 승인 항암제의 실제 target gene이 얼마나 들어있는지를 recall로 평가한다.
(2) Open Targets association score 비교
random보다 disease relevance가 높은 gene 선택되었고, 이는 PDGrapher가 뽑은 target gene은 단순히 모델 내부 score가 높은 gene이 아니라, 외부 질병-유전자 데이터베이스에서도 실제 폐암 관련성이 확인되는 gene임을 보여준다.
A549 lung cancer cell line에서 PDGrapher가 예측한 top-10 target gene이 Open Targets에서 폐암과 관련성이 높은지를 확인한다. 이때 비교군으로 random gene 10개를 사용한다.
평가 기준은 두 가지이다.
- Evidence strength : Open Targets가 계산한 질병-유전자 association score
- Number of evidence sources : 해당 질병-유전자 관계를 지지하는 독립적 evidence source 수
(3) 실제 암관련 데이터베이스에서 암 관련 유전자로 확인되는지 검증
PDGrapher가 예측한 유전자들이 여러 evidence source에서 등장하는 것을 보여주며, 이는 여러 독립적인 데이터에서 동시에 암 관련 유전자임을 보여준다
PDGrapher가 예측한 target들이 실제로 질병 관련 evidence(Open Targets)에 의해 지지 되는가를 보여준다.
- 세로축 : PDGrapher가 예측한 candidate target genes (6개 cell line에서 top-10 gene만 모은 것)
- 가로축 : evidence source 종류
- 각 점
- 색 : 연관성 강도
- 크기: evidence 개수
(4) ranking과 biological evidence의 연관성
rank가 높을수록 disease 연관성이 높았다.
다음 두 가지 질문을 통해 연관성을 확인하다
- 상위 rank gene이 하위 rank gene보다 더 실제 치료 표적인가?
- 어디까지를 candidate gene으로 믿어도 되는가? (top 10? top 50?)
실험 결과는 랜덤하게 gene ranking이 나오지 않고 랭크가 높을수록 진한 색을 가지고 rank가 낮을수록 연한 색을 보였으며, top-100 안 gene 대부분이 강한 biological evidence를 가지고 있음을 확인할 수 있다.
Ablation Study
(1) Supervision loss & Cycle loss 중 무엇이 중요한가
supervision loss : $CE(\mathcal{U}', f_p(x^d,x^t)$
cycle loss : $CE(x^t,f_r(x^d,f_p(x^d,x^t)))$
다음 세 가지 모델을 비교한다
| 모델 | 사용 loss |
| PDGrapher-Cycle | Cycle only |
| PDGrapher_Super | Supervision only |
| PDGrapher-SuperCycle | 둘 다 사용 |
Perturbagen prediction에서는 PDGrapher-Super가 최고 성능을 보였지만, PDGrapher-Cycle은 최악의 성능을 보였다.
state reconstruction에서는 PDGrapher-Cycle이 최고 성능을 보였다. 이 두 실험 결과를 바탕으로 두 목적을 동시에 만족하는 PDGrapher-SuperCycle을 선택한다.
(2) PPI causal graph 품질에 민감한가?
(2 - 1) STRING confidence threshold 변화 실험
STRING database는 gene-gene interaction edge마다 confidence score를 제공한다. threshold를 높여 low-confidence interaction을 제거해 network를 sparse하게 만들어도 PDGrapher 성능이 안정적으로 유지됨을 보여준다
(2-2) bridge edge 제거 실험
그래프에서 영향이 큰 영역을 연결하는 핵심 edge로 이 edge를 제거하면 구조가 크게 변하게 된다. 하지만 실제 실험 결과 성능 변화는 거의 없었으며(특히 bridge edge를 100% 제거하여도 점수가 0이 되지 않음 ➡️ graph backbone이 붕괴돼도 모델이 완전히 실패하지 않음), 그 이유는 bridge edge는 구조적으로 중요하지만 개수가 적어서 전체 prediction 영향이 제한적이기 때문이다.
(2-3) random edge 제거 실험
edge를 무작위로 제거하여 중요한 edge도 제거되고 덜 중요한 edge도 제거하엿다. 기존 실험과 달리 성능이 많이 감소하였고, 이는 pDGrpaher는 그래프 구조 자체는 활용하지만 특정 edge 하나에 의존하지 않고 전체 connectivity 패턴을 사용함을 보여준다.
(3) Latent confounder에 robust 한가?
Latent confounder : 관측되지 않았지만 여러 gene expression을 동시에 바꾸는 숨겨진 요인
PDGrapher는 다음과 같은 기본 가정을 한다
- PPI network $\approx$ causal graph
- no unobserved confounders
즉, 관측되지 않은 숨겨진 변수가 gene expression에 영향을 주지 않는다고 가정한다.
하지만 현실에서는 다양한 latent confounder가 항상 존재한다. 그래서 실제 latent confounder가 존재할 때 PDGrapher의 성능이 얼마나 유지될지 실험하였다.
결과를 보면 confounder가 증가할수록 성능이 감소했다. 하지만 causal structure는 유지하고 expression만 망가뜨렸는데도 성능이 완전히 붕괴되지 않았다(전부 confounder가 적용돼도 0이 아님).
Latent Confounder를 어떻게 시뮬레이션했는가?
다음과 같은 인위적인 방식으로 생성하였다
- random subset의 gene module 선택 (gene module은 서로 기능적으로 관련되어 있고 network에서 서로 촘촘히 연결된 gene 집합)
- 서로 다른 평균과 분산을 가진 가우시안 노이즈 추가
즉, $x_i \rightarrow x_i + \epsilon$ 형태로 gene expression을 교란하였다.
- x축 : confounder 영향을 받는 gene module 비율 (오른쪽으로 갈수록 많은 데이터가 confounder됨)
즉, PDGrapher가 gene expression noise에만 의존하지 않고 causal graph 구조 정보를 함께 활용한다는 것을 의미한다.
(4) Healthy control 데이터(disease intervention data)가 반드시 필요한가?
healthy control 없이도 성능이 크게 떨어지지 않는다. 즉, PDGrapher는 disease intervention data에 크게 의존하지 않는다.
위에서 설명했듯이 논문에서 사용하는 학습 데이터는 두 종류이다. (response prediction module)
- disease intervention data : $<x^h, \mathcal{U}, x^d>$
- treatment intervention data : $<x^d, \mathcal{U}',x^t>$
데이터셋에서 설명했듯이 현실에서 많은 cell line에 대해 healthy control sample이 없다. 그래서 대부분 diseased, treated만 존재한다. 이러한 상황에서 PDGrapher가 healthy control 없이도 잘 동작되는지를 실험하였다.
Objective
논문은 표현형 기반 신약 발굴 문제를 다음처럼 정의한다
diseased sample이 주어졌을 때 treated 상태로 이동시키는 gene target set을 찾는 문제
즉, $x^d \rightarrow \hat{x}^t \approx x^t$ 이 변화를 만드는 gene set을 예측하는 문제이며, 기존처럼 single target이 아닌 $\mathcal{U} = \{g_1,g_2,...,g_k\}$ 형태의 combinatorial prediciton problem이다.
Contribution
- combinatorial target prediction 가능 : single-target 접근을 넘어 multi-gene intervention 설계 가능
- personalized therapy 가능성 : 입력이 $x^d$이므로 개별 환자 expression profile 기반 target 예측이 가능하다
- gene dependency 반영 : 기존 one drug - one gene 처럼 보지 않고 유전자 네트워크 상호작용을 함께 고려하여 multiple genes - disease phenotype으로 복합적인 질환 메커니즘에 맞는 perturbation 조합 예측 가능
Limitation
- Unobserved confounder 가정
- PDGrapher는 unobserved confounder(유전자 발현 변화는 우리가 알고 있는 gene interaction만으로 설명된다)가 없다고 가정하며, 이 가정은 경험적으로 검증하기 어려움 (실제로는 많은 biological 변수들이 영향을 미침)
- PPI/GRN을 causal graph의 proxy로 사용하는 한계
- PPI와 GRN을 causal gene network의 근사로 사용하지만 이들 네트워크는 본질적으로 noisy하고 incomplete
- 즉, 모델이 사용하는 gene network가 실제 causal 구조와 완전히 일치하지 않음
- Transcriptomics 데이터에 대한 의존성
- 모델이 transcriptomic 데이터에 기반하므로 다른 phenotype modality를 반영하지 못함 ( 표현형 변화의 일부만 보고 예측 = 정보 부족 상태에서 예측)
- protein activity, metabolite 상태, cell morphology, chromatin 구조 같은 정보도 중요함. 즉 gene expression은 같아도 protein activation 상태가 다를 수 있음
- Healthy control cell line mismatch 문제
- A549와 NL20은 서로 다른 anatomical origin(baseline이 다름)을 가지므로 baseline gene expression 차이에 따른 bias가 발생할 수 있다
- 모델이 차이를 disease effect로 생각하지만 실제로 이 차이는 tissue origin difference일 수도 있다. 즉, disease signal이 왜곡될 수 있다.