면역관문억제제(immune checkpoint inhibitors, ICIs)는 진행성 비소세포폐암(non-small cell lung cancer, NSCLC)에 대한 표준 치료로 채택되고 있으나, 반응 이질성에 대한 유전적 결정인자는 여전히 명확하지 않다[1]. 대부분의 조혈계 계열(hematopoietic lineages)은 종양 병인과 면역치료 과정에서 역동적으로 변화하므로[2], 생식세포(germline) 변이가 그들의 전사체(transcriptome)를 어떻게 조절하는지 규명하는 것은 중요하다. 특히 단일세포 RNA 염기서열분석(single-cell RNA sequencing, scRNA-seq)과 통합된 eQTL(expression quantitative trait loci) 분석은 단일세포 수준에서 유전자 조절 지도를 작성할 수 있게 해준다[3, 4]. 구체적인 방법론은 부록(Supplementary Materials)에 기술되어 있다. ICI 치료 중 생식세포 변이가 면역 유전자 조절에 어떤 영향을 주는지 조사하기 위해, 우리는 단일세포-eQTL(single-cell-eQTL, sc-eQTL) 분석과 전사체 네트워크 프로파일링을 수행하였다. 말초혈액단핵세포(peripheral blood mononuclear cells, PBMCs)는 항-프로그램드 세포사멸 단백질 1(anti-programmed cell death protein-1, PD-1) 또는 프로그램드 데스 리간드 1(programmed death-ligand 1, PD-L1) 치료를 받은 NSCLC 환자 73명에서 치료 전(baseline)과 치료 후 1~5주 시점에 수집하였다(그림 1A, 부록 표 S1). scRNA-seq와 SNP array 데이터를 통합하여, 세포 유형별(cell-type-resolved) sc-eQTL과 유전자 네트워크를 분석하였다(그림 1A-B). sc-eQTL-조절 전사체 네트워크는 NSCLC 환자의 ICI 치료 결과를 설명한다. (A) 실험 및 분석 개요. 치료 전과 치료 중 73명의 ICI 치료 NSCLC 환자의 PBMC에서 sc-eQTL 매핑과 유전자 네트워크 분석을 수행하여, 폐암 및 치료 특이적 eQTL, 전사체 네트워크, GNN(refined) 유전자 세트에 기반한 임상적 관련성, 그리고 유전자 네트워크 내 eQTL의 위치를 평가하였다. (B) 이 코호트의 단일세포 면역 지형과 치료 관련 구성. 왼쪽: 559,714개의 면역 세포에 대한 UMAP 시각화는 여덟 개의 면역 클러스터를 드러냈다. 오른쪽: 세포 유형 비율은 치료(BL: baseline vs. TRT: treatment) 및 임상 반응(DCB vs. NCB)에 따라 달랐다. (C) 세포 유형 및 조건에 따른 sc-eQTL의 분포. 업셋 플롯은 공유 패턴에 따른 eQTL 개수를 나타낸다. 점 플롯에서 “Shared”는 두 개 이상 클러스터에서 공유되는 eQTL을 의미하고, “Cell type-specific”는 해당 클러스터에 제한된 eQTL을 의미하며 검은 점은 각 특정 클러스터를 나타낸다. 막대 플롯은 치료 특이적 분포를 표시한다. 주목할 점은 CD8+ T 세포와 단핵구가 상당한 양의 고유한 eGenes에 기여했다. (D) 대표적인 치료 특이적 eQTL. 치료 특이적 eQTL은 baseline과 치료 조건에서 유전자형(genotype) 의존적 유전자 발현을 비교하여 확인하였다(추가 세부사항은 부록 그림 S2C 참조). 박스플롯은 치료 조건 전반에 걸친 TNF(단핵구), TNFRSF1A, PRF1(CD8+ T 세포)의 유전자형 의존적 발현을 보여준다. 각 연관은 lfsr < 0.05(mashr) 기준을 충족했으며, posterior β 값이 표시된다. (E) 폐암 특이적 eGenes의 확인. 베ン 다이어그램은 GTEx v8 및 OneK1K 혈액 eQTL 데이터셋과의 eGenes 중복을 보여준다. 총 702개의 eGenes(4.6%)가 본 코호트에만 고유하여 질병 특이적 전사체 조절을 시사한다. 본 연구에서만 확인된 eQTL은 빨간색으로 강조하였다. (F) ICI 반응을 설명하는 공발현(co-expression) 모듈의 확인. 모듈 활성 차이에 대한 화산(Volcano) 플롯(기준: DCB의 baseline)을 WGCNA 모듈 전반에서 제시하였다. 모듈은 WGCNA를 사용하여 개별 면역 클러스터로부터 구성하였다. CD8-brown 모듈은 반응 그룹 간 baseline에서 유의한 차등 활성( |log2FC| > 1, FDR < 0.05)을 보였다. (G) CD8-brown 모듈의 조절적 관련성. CD8+ T 세포에서 CD8-brown 모듈과 EOMES 및 TBX21 조절자(regulons)의 상관 분석을 수행했다. 색 척도는 세포 수의 밀도를 나타낸다. (H) 종양 미세환경(tumor microenvironment, TME) CD8+ T 세포에서의 CD8-brown 모듈 활성 풍부화. 18명의 NSCLC 환자로부터 얻은 종양 및 전이성 림프절(scRNA-seq) 데이터셋에서 NK/T 세포를 확인하고, CD8+ T 특이 WGCNA 네트워크를 이 데이터셋에 투영하였다. 바이올린 플롯은 NCB 종양의 effector CD8+ T 아형(CD8+ TEM, CD8+ TEMRA)에서 CD8-brown 모듈 활성이 상승했음을 보여준다. ****P < 0.0001, Wilcoxon rank-sum test. (I) CD8-brown+ CD8+ T 세포 분율의 임상적 관련성. OS 및 PFS는 CD8+ T 세포 중 CD8-brown 활성이 높은 세포의 비율(>75%; CD8-brown 활동은 CD8-brown core gene score가 0.36 이상으로 정의됨)에 따라 분층한 Kaplan-Meier 곡선을 제시하였다. 아래 표의 ‘Ratio-high’에 해당하는 CD8-brown+ CD8+ T 세포 비율이 더 높은 환자들은 유의하게 더 나쁜 예후를 보였다(OS: P = 0.035, PFS: P = 0.018). 아래 표는 각 군의 환자 수를 나타낸다. (J) eGenes는 CD8+ T 세포에서 더 높은 중심성을 보인다. 바이올린 플롯은 baseline 및 ICI 치료 후 CD8+ T 세포에서 eGenes와 비-eGenes 간의 rank-normalized WGCNA 모듈 연결성(kME)을 비교한다. eGenes는 더 중심적인 위치에 더 많이 분포했다(P < 0.001). (K) 모듈 소속에 따른 eQTL 효과 크기. WGCNA 모듈에 할당된 유전자에 대한 eQTL의 효과 크기(|β|)는 더 작았으며, 이는 고도로 연결된 유전자에 기능적 제약이 있음을 시사한다. P < 0.001. (L) 모듈 중심성과 eQTL 효과 크기 사이의 역상관. kME rank 분위수(quintiles) 전반의 박스플롯은 baseline과 치료 모두에서 중심성이 증가할수록 eQTL 효과 크기가 감소함을 보여준다. 유의한 경향(P < 0.001, 선형 회귀)은 핵심 네트워크 유전자의 제한된 유전적 조절을 나타낸다. eQTL 유의성은 mashr의 local false sign rate(lfsr < 0.05)를 이용해 정의하였다. 군 간 비교는 Wilcoxon rank-sum 또는 선형 회귀로 평가했으며, 다변량 생존 모델은 임상 공변량으로 보정하였다. 약어: BL, baseline; CD8+ TCM, CD8+ central memory T cell; CD8+ TEM, effector memory CD8+ T cells; CD8+ TEMRA, CD8+ terminally differentiated effector memory T cells with expressing CD45RA; DCB, durable clinical benefit; DC, dendritic cell; eGenes, eQTL-regulated genes; EOMES, eomesodermin; eQTL, expression quantitative trait loci; FC, fold change; FDR, false discovery rate; GNN: graph neural network; GRN, gene regulatory network; GTEx, Genotype-Tissue Expression project; ICI, immune checkpoint inhibitor; kME, module eigengene; lfsr, local false sign rate; LN, lymph node; NCB, non-durable clinical benefit; NK, natural killer cell; NSCLC, non-small cell lung cancer; OneK1K, OneK1K sc-eQTL study; OS, overall survival; PBMC, peripheral blood mononuclear cells; PFS: progression-free survival; PRF1, perforin 1; sc-eQTL, single-cell eQTL; TBX21, t-box transcription factor 21; TNF, tumor necrosis factor; TNFRSF1A, TNF receptor 1A; TRT, treatment; UMAP, uniform manifold approximation and projection; WGCNA, weighted gene co-expression network analysis. 품질 관리와 pseudobulk 집계를 수행한 후, 9,147개의 eQTL 쌍(발현 조절 SNP: expression-regulating SNPs [eSNPs])을 확인했으며, 이는 8개 면역 세포 클러스터와 치료 조건 전반에 걸쳐 3,616개의 혈액 발현-조절 유전자(eGenes)와 연관되어 있다(그림 1B-C, 부록 그림 S1A, 부록 표 S2). 선행 연구[3, 4]와 일관되게, eGene 개수는 세포 풍부도(cell abundance)와 상관했고, eSNP는 조절 요소에서 풍부화되어 있었다(부록 그림 S1B-D). Multiadaptive shrinkage[5]는 치료- 및 세포 유형 의존적인 조절이 뚜렷함을 보여주었으며, 여기에는 245개의 치료 특이적 eQTL이 포함되었다(부록 그림 S2A-D, 부록 표 S3-S4). 예를 들어, 종양 괴사 인자(tumor necrosis factor, TNF)는 치료 후 단핵구에서 조절되었고(posterior β = 1.17), TNF 수용체 1A(TNF receptor 1A, TNFRSF1A)는 CD8+ T 세포에서 baseline 조절(β = 1.10)을 보였는데, 이는 ICI 치료 중 면역 유전자 발현이 생식세포 변이에 의해 좌우될 수 있음을 시사한다(그림 1D). 추가 예로는 baseline의 CD8+ T 세포에서 세포독성 매개체인 perforin 1(PRF1)과 granzyme B(GZMB)가 포함된다(그림 1D, 부록 그림 S2E). 본 결과를 검증하기 위해 두 가지 상호보완적 분석을 수행하였다. 첫째, 자가면역 및 혈액 특성에 대한 전장유전체 연관 연구(genome-wide association study, GWAS) 좌위와의 공국소화(colocalization) 분석에서 중복이 관찰되었다(PP.H4 > 0.6). 이는 가능한 공유 조절 기전을 시사한다(부록 그림 S3, 부록 표 S5-S6). 둘째, 외부 eQTL 연구와의 비교에서 본 연구 특이적 eQTL(이하 폐암 특이적 eQTL로 지칭됨)은 암 및 면역 반응 관련 경로에서 풍부화되어 있었다(그림 1E, 부록 그림 S4A). 이는 만성 면역 적응을 반영한다. 그럼에도 불구하고, 1M-scBloodNL 연구[3]에서 건강 공여자에게서 관찰된 대부분의 조절 양상은 보존되어 있었다(부록 그림 S4B, 부록 표 S7-S8). eQTL 매핑을 넘어, 우리는 weighted gene co-expression network analysis(WGCNA)를 사용하여 각 주요 면역 세포 유형에 대한 공발현 네트워크를 구축했으며[6], 이를 통해 세포 유형 특이적 조절과 전사체를 비교할 수 있었다(부록 표 S9-S10). 확인된 모듈 중 CD8-brown 모듈은 비지속 임상적 유익(non-durable clinical benefits, NCB) 군에서 baseline 활성이 증가되어 있었다(그림 1F, 부록 그림 S5A). 이 모듈은 세포독성 유전자 풍부(PRF1, apolipoprotein B mRNA-editing enzyme, catalytic subunit 3G [APOBEC3G], 및 GZMB)였고 NCB 군에서 특히 분화된 CD8⁺ T 아형들에서 고발현되었다(부록 그림 S5B-D). 모듈 활성은 건강한 개인, 종양 환자, ICI 치료 종양 데이터셋의 외부 혈액 데이터로부터 지지되었다(부록 그림 S5E-F). 잠재적 조절자를 탐색하기 위해, 우리는 단일세포 조절 네트워크 추론 및 클러스터링(SCENIC)[7]을 적용하여 CD8+ T 세포 분화 및 소진(exhaustion)과 연결된 가설적 조절자로 eomesodermin(EOMES)과 t-box transcription factor 21(TBX21)을 확인했다[8]. 이들의 regulon 활성은 CD8-brown 모듈(그림 1G, 부록 그림 S6A-B)과 effector memory CD8+ T(CD8+ TEM) 세포의 풍부도(부록 그림 S6C-D) 모두와 일치했다. 핵심 유전자(EOMES, TBX21, 그리고 interleukin-2 receptor, beta subunit [IL2RB])는 DCB(durable clinical benefit) 군에 비해 NCB 군에서 상향 조절되었다(부록 그림 S6E). 특히 TBX21과 EOMES는 본 분석에서 eGenes가 아니었다. 이러한 결과를 바탕으로, 우리는 단백질-단백질 상호작용(PPI) 네트워크, 유전자 온톨로지(GO) 주석, 공발현 데이터를 통합하여 CD8-brown 모듈을 graph neural network(GNN) 기반 프레임워크로 정교화하였다(부록 그림 S7A-B). 정교화된 서브네트워크는 더 일관된 CD8⁺ T 세포 분화 프로그램을 포착했으며, STRING(Search Tool for the Retrieval of Interacting Genes/Proteins)에서 원래 WGCNA 유래 유전자 세트보다 더 강한 위상적 연결성을 보였다(노드 밀도: 0.333 vs 0.183; 부록 그림 S7C-D, 부록 표 S11)[9]. CD8+ T 세포의 canonical 경로에 대한 경로 농축 결과는 이러한 정교화를 지지했다(부록 그림 S7E-F). 이 모듈은 종양 미세환경(TME) 및 ICI 치료 NSCLC 환자의 외부 scRNA-seq 데이터셋(부록 그림 S8)[10]과 암 환자의 혈액 데이터셋(부록 그림 S9) 등 다양한 외부 scRNA-seq 데이터셋에서 잘 보존되었다. 이 모듈의 활성은 비반응자(non-responders)에서 유의하게 강화되었으며, 특히 CD8+ TEM과 같은 분화된 CD8+ T 세포 아형에서 두드러졌다(그림 1H, 부록 그림 S8-S9, 부록 표 S12). 이는 이러한 전신적 세포독성 시그니처가 TME에서 재현(recapitulate)되며 ICI 반응과 연관된 면역 상태를 반영할 수 있음을 시사한다. 이 모듈이 ICI 반응과 강하게 연관되므로, CD8-brown 모듈로부터 정교화된 핵심 유전자가 생존 결과를 설명할 수 있는지 평가하였다. 우리는 네트워크 중심성이 높은 15개의 GNN-우선화 CD8-brown 유전자로부터 핵심 유전자 점수(core gene score)를 도출했으며(부록 방법, 부록 그림 S7, 부록 표 S13), 이 점수는 모듈 활성과 강하게 상관했다(부록 그림 S10A). 그 후 그리드 서치와 부트스트래핑(부록 그림 S10B-C)을 사용하여 환자 분류 임계값을 설정했다. 즉 CD8+ T 세포 중 CD8-brown 핵심 유전자 점수가 높은 세포가 >75%인 환자(>0.36)는 유의하게 더 짧은 전체 생존(OS, P = 0.035)과 무진행 생존(PFS, P = 0.018)을 보였다(그림 1I, 부록 그림 S10B-C). 이 생존 연관성은 임상 공변량을 보정한 다변량 Cox 모델에서도 유의하게 유지되었다(부록 그림 S10D). 전반적으로, 본 분석은 TBX21–EOMES 조절축이 CD8-brown 모듈을 구동할 수 있으며, 이는 PRF1 및 GZMB와 같은 세포독성 eGenes를 포함한 CD8+ T 분화를 나타내고 결과적으로 ICI에 대한 불량 반응자를 계층화한다는 점을 시사한다. 마지막으로, 생식세포 변이가 ICI 반응과 연관된 면역 네트워크에 어떻게 영향을 주는지 평가하기 위해, 공발현 및 조절 네트워크 내에서 eQTL의 분포를 살펴보았다(그림 1J-L, 부록 그림 S11A-C). eGenes는 중심 네트워크 위치에서 풍부화되어 있었다(부록 그림 S11B). 그러나 PRF1, GZMB 또는 TNFRSF1A를 포함한 단일 eQTL이 CD8-brown 모듈에 지배적인 조절 효과를 나타내지는 않았다. 대신 중심성이 더 높은 유전자 또는 WGCNA 할당이 높은 유전자는 더 작은 eQTL 효과 크기와 연관되었다(그림 1K-L). 이는 핵심 면역 분화 경로에 기능적 제약이 있음을 시사한다. 그럼에도 불구하고 몇 가지 제한점이 존재한다. 첫째, 반응을 DCB와 NCB로 이분화하면 이질적인 결과가 가려질 수 있으며, 치료로 유도된 전사체 변화는 샘플 간 이질성보다 작은 것으로 보였다(부록 그림 S12). 둘째, 단핵구 특이적 eQTL은 검출되었고 단핵구-ICI 상호작용은 잘 알려져 있으나, ICI 관련 공발현 모듈에서의 기여가 미미하여 후속 분석에서 제외했다(부록 그림 S13). 셋째, 세포 유형 및 조건 특이적 eQTL과 공발현 네트워크에 초점을 맞춘 것은 면역 기능에 영향을 주는 세포 간 상호작용과 trans 효과를 간과할 수 있다. 마지막으로, CD8-brown 핵심 유전자 점수는 본 코호트에서 생존을 계층화했지만(c-index = 0.68-0.71 for PFS/OS, 부록 그림 S10D), 그 예측 가치는 독립 데이터셋에서 검증이 필요하다. 종합하면, 본 연구는 sc-eQTL 매핑을 통해 혈액-특이 3,616개 및 폐암-특이 702개의 eGenes를 확인했으며, ICI 반응의 비반응자와 연관된 CD8+ T 세포 분화와 관련된 전사체 모듈(CD8-brown)을 규명하였다. 공발현 네트워크에 대한 eQTL의 제한된 영향은 면역 전사체에 기능적 제약이 있음을 시사한다. 유전적 변이를 세포독성 네트워크 활성 및 임상적 결과와 연결함으로써, 본 분석은 전이성 NSCLC에서 ICI 치료 하의 전신 면역을 이해하기 위한 프레임워크를 제공한다. 향후 더 큰 코호트와 핵심 조절 변이—예를 들어 PRF1 및 GZMB에 영향을 주는 변이—에 대한 실험적 검증이 포함된 연구는 인과적 기전을 명확히 하고 개인맞춤형 ICI 전략을 정교화하는 데 기여할 수 있을 것이다. Hyungtai Sim: data curation, formal analysis, validation, investigation, writing-original draft. Geun-Ho Park: data curation, formal analysis, investigation. Woong-Yang Park: investigation. Se-Hoon Lee: conceptualization, investigation, writing-original draft. Murim Choi: conceptualization, investigation, writing-original draft. 본 연구에 참여해 준 환자들께 감사드린다. 저자들은 경쟁 이해관계가 없다고 선언한다. 본 연구는 부분적으로 한국연구재단(Korean Research Foundation)의 보조금(NRF-2021R1A2C3014067, NRF-RS-2023-00207857 to Murim Choi; 및 NRF-2020R1A2C3006535, NRF-RS-2025-00519956 to Se-Hoon Lee)과, 한국보건산업진흥원(Korea Health Industry Development Institute)을 통한 한국 보건기술 R&D 프로젝트(HR20C0025 to Se-Hoon Lee), 그리고 삼성서울병원(Samsung Medical Center)의 Future Medicine 20*30 프로젝트(SMX1230041 to Se-Hoon Lee)의 지원을 받아 수행되었다. 본 연구는 사람 유래 시료를 대상으로 하였으며, 삼성서울병원 임상시험심사위원회(Institutional Review Board)에서 허가를 받았다(No. 2018-04-048, 2022-01-094). 모든 참가자는 연구에 등록되기 전 informed consent 과정을 거쳤다. 행렬 형태의 발현 데이터와 scRNA-seq에 대한 메타데이터는 K-BDS에 accession number KAP241028로 예치하였다. 본 원고에서 분석한 추가 데이터셋은 다음과 같이 공개되어 있다. Gene Expression Omnibus(GEO; www.ncbi.nlm.nih.gov/geo/)에서 accession number GSE200996(HNSCC 데이터셋), GSE123814(TNBC 데이터셋), GSE205335(NSCLC 데이터셋) 및 CellxGene의 AIDA Phase 1 데이터셋(https://cellxgene.cziscience.com/collections/ced320a1-29f3-47c1-a735-513c7084d508)이다. 이 연구의 결과를 뒷받침하는 기타 모든 데이터(파생 데이터셋 포함)는 원고 및 그 부속 파일에 제공되거나 GitHub 저장소(https://github.com/snu-mchoi-lab/LungCancer_ICI_sceQTL_2024)에 이용 가능하다. 출판사는 저자들이 제공한 부가 정보의 내용 또는 기능에 대해 책임지지 않는다. 누락된 콘텐츠 이외의 문의사항은 해당 논문의 교신저자에게 전달되어야 한다.

변이 마우스는 면역억제 효과를 유지하면서, 글루코코르티코이드가 유발하는 근력 약화 및 섬유화(섬유형성) 효과를 상쇄함으로써 근육 표현형이 개선되었다. 우리의 분석은 DMD에 대해 현재 사용 가능한 치료 옵션으로도 용이하게 표적화할 수 있는 병인 기전을 밝혀낸다.