전사체 분석을 이용한 유방암 내 허브 유전자, 단일 염기 다형성 및 잠재적 약물 표적 식별

Summary

본 연구는 HER2 양성 유방암에서 치료 내성과 관련된 미토콘드리아 산화 스트레스 관련 유전자를 조사하였습니다. 전사체 데이터셋, 통합 생물정보학, 임상 데이터(n = 4,929)를 사용하여 저자들은 MTHFD2와 PRDX3를 잠재적 예후적 가치를 가진 주요 차별 발현 유전자로 확인했습니다. 실리코에서. 분석 결과, 유전적 변이의 기능적 효과가 시사됩니다.

Abstract

HER2 양성(HER2⁺) 유방암은 라파티닙과 같은 치료법에 내성을 보이는 경우가 많으며, 이는 미토콘드리아 대사 및 산화환원 재프로그래밍과 관련된 가능성이 있습니다. 이 연구는 미토콘드리아 산화 스트레스 관련 유전자(MOS-DEGs)를 내성의 바이오마커로 식별하고, 이들의 기능적 nsSNPs와 구조적 영향을 규명하는 것을 목표로 했습니다. RNA-seq 데이터셋(GSE231524, GSE231525)은 DESeq2를 사용하여 DEG를 분석했으며, DEGs는 미토콘드리아 산화 스트레스 유전자와 교차하여 MOS-DEG를 얻었습니다. 기능 풍부화, PPI 분석, ROC 분석, 발현 프로파일링, 생존 분석이 수행되었습니다. nsSNP는 여러 예측 도구를 사용해 평가되었고, 구조적 영향은 2차 및 3D 모델링을 통해 평가되었습니다. 통합 전사체 및 임상 분석을 통해 MTHFD2와 PRDX3가 미토콘드리아 산화환원 대사에서 상반되는 조절 역할을 하는 중심 MOS-DEG로 확인되었습니다. MTHFD2는 유의미하게 상향 조절되어 예후가 좋지 않은 것으로 나타났습니다(HR = 1.53, p = 1.1 × 10⁻¹⁶), PRDX3는 보호 발현 패턴(HR = 0.73, p = 7.7×10⁻¹⁰)을 보였습니다. ROC 분석은 이들이 치료 반응의 예측 인자로서 잠재력을 시사했습니다. nsSNP 분석 결과 MTHFD2에서 5개의 해로운 변이, PRDX3에서 4개의 해로운 변이가 발견되었으며, 그중 rs1471336772(MTHFD2)와 rs747786383(PRDX3.)이 가장 병원성이 높게 확인되었습니다. 이 변이들은 SIFT ≤ 0.05, PolyPhen-2 ≥ 0.85, 해로운 CADD 점수, 높은 REVEL 점수 등 여러 계산 점수 체계에 근거해 손상을 초래할 것으로 예측되었으며, 촉매 또는 산화환원 활성 도메인 내에 위치했습니다. 구조 모델링은 이러한 치환이 구조를 불안정하게 만들고, 금속 및 보조인자 결합 부위를 교란하며, NADPH 재생 또는 티오레독신 의존성 과산화효소 활성에 영향을 줄 수 있음을 시사했습니다. 분자 역학 예측은 구조적 안정성 상실과 유연성 변화 가능성을 시사하여 미토콘드리아 산화환원 조절의 기능 장애 가능성을 시사했습니다. 본 연구는 HER2⁺ 유방암에서 MTHFD2와 PRDX3가 미토콘드리아 산화 스트레스 조절 인자임을 확인한다. 이 유전자들 내 해로운 nsSNP는 산화환원 균형 변화에 기여하여 대사 적응(MTHFD2)과 항산화 방어(PRDX3)에 영향을 미칠 수 있습니다.

Introduction

유방암은 전 세계 여성에서 가장 흔하게 진단되는 암이며, 조기 발견과 표적 치료가 상당한 진전을 이루었음에도 불구하고 암 관련 사망의 주요 원인입니다. ^1,2,3. 유방 종양 내 분자 이질성은 다양한 임상 결과와 치료 반응의 근간이며, 정밀 종양학에 지속적인 도전 과제를 제기합니다 ^4,5. 알려진 분자 아형 중 인간 표피 성장인자 수용체 2(HER2) 양성 유방암은 전체 유방암의 약 15–20%를 차지하며, HER2 유전자 증폭과 HER2 수용체 티로신 키나아제의 과발현이 특징입니다 ^6,7. 트라스투주맙과 라파티닙과 같은 HER2 표적 치료법이 환자의 예후를 현저히 개선했음에도 불구하고, 이들에 대한 1차 및 후천성 내성이 자주 나타나 종양 재발과 질병 진행으로 이어집니다 ^8,9. 따라서 HER2 치료 내성을 유발하는 분자 기전을 이해하는 것은 임상 종양학에서 여전히 주요 미충족 과제로 남아 있습니다.

새로운 증거들은 미토콘드리아 기능 장애와 산화 스트레스가 유방암에서 약물 내성을 매개하는 데 중요한 역할을 한다는 것을 시사합니다. 미토콘드리아는 ATP 생성에서의 정규적 역할을 넘어, 세포자살자사, 산화환원 신호, 대사 가소성의 핵심 조절자로, 암세포가 치료 압력을 견디기 위해 이용하는 과정들입니다. 특히, 내성 HER2⁺ 암세포는 산화 인산화(OXPHOS), 향상된 활성 산소종(ROS) 완충, 미토콘드리아 생합성 변화로 대사적으로 전환하는 모습을 보입니다. 이러한 적응은 세포가 HER2 신호 전달 경로가 지속적으로 억제되어도 생존 신호를 유지하고 세포자멸사를 회피할 수 있게 합니다. 따라서 미토콘드리아 복합체 소단위체와 산화환원 조절 효소를 암호화하는 유전자는 저항성 극복을 위한 잠재적 바이오마커이자 치료 표적이 됩니다.

이와 병행하여, 미토콘드리아 또는 산화 스트레스 관련 유전자의 비동의어 단일 염기체 다형성(nsSNPs)은 단백질 구조와 기능을 변형시켜 효소 활성, 약물 반응, 질병 감수성에 영향을 미칠 수 있습니다¹⁰. 실 리코를 이용한 nsSNP의 해로운 영향에 대한 계산적 예측.SIFT, PolyPhen-2, CADD와 같은 도구는 암의 이질성과 치료 내성에 기여할 수 있는 기능적 변이를 신속히 식별할 수 있는 수단을 제공합니다. 따라서 전사체 및 돌연변이 데이터를 통합함으로써 전사 및 유전체 수준 모두에서 유전자 조절 이상과 구조 변이를 다차원적으로 관찰할 수 있습니다.

생물정보학과 시스템 생물학의 최근 발전은 대규모 발현 프로파일링과 네트워크 기반 분석을 통해 잠재적 드라이버 유전자를 고처리량으로 식별할 수 있게 했습니다. 유전자 발현 옴니버스(GEO)와 같은 공개 저장소는 실험 모델과 환자 코호트에 걸친 분자 변화를 포착하는 귀중한 RNA 시퀀싱 데이터셋을 제공합니다. DESeq2, clusterProfiler, STRING–Cytoscape 네트워크 분석과 같은 분석 도구와 결합하여, 이러한 자원들은 차별발현 유전자(DEG), 경로 풍부화, 허브 유전자 발견의 포괄적인 특성 분석이 가능합니다. 중요하게도, 미토콘드리아 산화 스트레스 유전자 서명을 DEGs와 통합하면 HER2 표적 치료 내성의 새로운 분자 기전을 밝히고 치료 개입을 위한 후보 바이오마커를 식별할 수 있습니다.

HER3 구동 저항축은 최근 HER2 억제 이후의 주요 보상 경로로 주목받고 있습니다. HER3 활성화는 하류 PI3K/AKT 신호를 복원하여 세포 생존과 약물 내성을 촉진합니다(그림 1). ERK 신호의 음성 조절자인 DUSP6(이중 특이성 인산효소 6)가 이 적응 반응을 조절하는 데 관여하는 것으로 알려져 있습니다. DUSP6의 억제는 ERK 신호를 재활성화하여 HER3 매개 내성을 상쇄할 수 있으나, 미토콘드리아 산화 스트레스 및 대사 적응과의 상호작용은 아직 충분히 밝혀지지 않았습니다. 따라서 DUSP6 억제와 만성 라파티닙 노출 하에서 HER2⁺ 세포주(BT474 및 MDA-MB-453)의 비교 전사체 분석은 약물 내성의 분자적 결정 인자를 해부하는 이상적인 모델을 제공합니다.

그림 1. 미토콘드리아 기능 장애, 산화 스트레스, 유방암 진행을 연결하는 분자 지형의 개략적 표현. 이러한 통합된 분자 메커니즘은 유방암에서 미토콘드리아 산화 스트레스 매개 암 발생의 기저입니다. 특히 DUFS3, UQCRC1, COX4I1, SDHA, ATP5PO 와 관련된 미토콘드리아 전자전달 사슬 복합체의 조절 이상은 과도한 활성 산소종(ROS)의 생성을 초래하여 산화성 DNA 손상, 지질 과산화, 미토콘드리아 막 전위 저하를 초래합니다. 이러한 산화적 교란은 PI3K/AKT, MAPK, NF-ΚB와 같은 종양성 신호 경로를 활성화하고, p53과 BAX를 포함한 세포자멸사 조절인자를 억제하여 종양 세포의 생존, 증식, 면역 회피를 촉진합니다. 도표는 또한 미토콘드리아 산화 스트레스가 암 대사와 종양 미세환경에 미치는 이중 영향을 강조하며, HER2+ 유방암 상태에서 치료 저항성과 면역 조절에 기여한다는 점을 강조합니다. 이 통합 모델은 이후 미토콘드리아 드라이버 유전자와 잠재적 치료 표적을 규명하기 위한 전사체 및 SNP 기반 분석의 기전적 기초를 형성합니다. 이 그림의 더 큰 버전을 보시려면 여기를 클릭해 주세요.

본 연구는 HER2⁺ 유방암에서 허브 유전자, 기능적 nsSNP, 그리고 미토콘드리아 산화 스트레스와 관련된 잠재적 약물 표적을 체계적으로 규명하도록 설계되었습니다. GSE231524 및 GSE231525의 전사체 데이터를 엄선된 미토콘드리아 및 산화 스트레스 관련 유전자 집합과 통합함으로써, 이번 연구는 치료 내성에 기여하는 미토콘드리아 산화 스트레스 관련 차별발현 유전자(MOS-DEGs)를 정의하는 것을 목표로 했습니다. 후속 분석에는 유전자 온톨로지(GO) 및 KEGG 경로 농축, 단백질-단백질 상호작용(PPI) 네트워크 구축, 허브 유전자 우선순위 지정, Kaplan–Meier 면역치료 데이터베이스를 이용한 생존 상관관계 분석이 포함되었습니다. 더 나아가, 돌연변이 유발 기능 결과를 구분하기 위해 nsSNP 변이와 잠재적 구조적 영향에 대한 주요 유전자가 추가로 조사되었습니다.

전사체 프로파일링, 네트워크 생물학, 그리고 실리코 돌연변이 분석을 결합함으로써, 본 연구는 HER2⁺ 유방암에서 잠재적 미토콘드리아 바이오마커와 치료 표적을 발견하기 위한 통합적 틀을 제공합니다. 산화성 인산화 연관 허브 유전자와 그 기능성 SNP의 확인은 HER2 표적 약물 내성을 극복하고 환자 치료 결과를 개선하기 위한 정밀 치료 전략 개발의 길을 열 수 있습니다.

Protocol

데이터 검색 및 전처리

이 연구는 공개된 전사체 및 유전 데이터를 사용하여 수행되었으며; 인간이나 동물 대상은 직접적으로 참여하지 않았습니다. HER2⁺ 유방암 치료 내성과 관련된 전사체 데이터셋은 NCBI 유전자 발현 옴니버스(GEO) 데이터베이스에서 검색되었습니다 (https://www.ncbi.nlm.nih.gov/geo/)¹¹. 두 개의 RNA-seq 데이터셋인 GSE231524와 GSE231525는 HER3 유도 내성과 HER2+ 유방암 세포주(BT474 및 MDA-MB-453)에서 DUSP6 억제에 구체적으로 초점을 맞추고 있기 때문에 선정되었습니다. 이 데이터셋에는 라파티닙(1 μM) 노출과 DUSP6 눌러다운에서 유래한 부모, 약물 내성, 약물 내성 표현형이 포함되었습니다. 원시 카운트 매트릭스와 해당 메타데이터 파일은 RStudio(v4.3.2)의 GEOquery(v2.70.0)와 Biobase(v2.62.0) 패키지를 통해 접근되었습니다^. 각 데이터셋에 대해 두 가지 주요 대조를 정의하기 위해 메타데이터를 선별했습니다: GSE231524 대조군(BT474 부모, 0일차)과 약물 내성 및 약물 저항성 샘플(9일째–9개월)을 비교했고, GSE231525 대조군(스크램블드 siRNA)과 DUSP6 노크다운(DUSP6-KD)을 비교했습니다. 품질 관리와 데이터 정규화는 DESeq2 (v1.42.0) 프레임워크를 사용하여 수행되었으며, 이 프레임워크는 분산 안정화 변환(VST)을 적용하여 이종편차성을 줄이고 샘플 간 비교 가능성을 보장했습니다. GGPLOT2(v3.5.0)와 pheatmap(v1.0.12)을 사용하여 데이터 균일성을 확인하고 차별 표현 분석 전에 잠재적 이상치를 식별하기 위해 데이터 분포 및 군집화 패턴을 시각적으로 평가하였다^13,14.

본 연구에서는 세포주 전사체 데이터셋에서 얻은 결과와 환자 유래 임상 데이터셋에서 얻은 결과가 명확히 구분되었습니다. 세포주 데이터는 주로 탐색적 분석에 사용되었으며, 여기에는 차별 발현 유전자의 식별과 통제된 실험 모델에서의 예비 기전적 통찰 도출이 포함되었습니다. 반면, 환자 유래 데이터셋은 유전자 발현 패턴의 외부 검증과 임상적 관련성 평가, 예후 평가를 위해 사용되었습니다. 따라서 세포주 모델과 임상 코호트의 결과는 과도한 일반화를 피하고 모든 결과에 적절한 번역적 맥락을 보장하기 위해 별도로 해석됩니다.

차등 유전자 발현 분석

대조군과 치료 조건 간에 유의하게 조절된 유전자를 식별하기 위해 차별 발현 분석이 수행되었습니다. 정규화된 카운트는 DESeq2에 통합된 조정 회귀 모델(ARM)을 사용하여 log₂ 폴드 변화와 통계적 유의성을 정확히 추정하였습니다. 설계 공식은 ~조건으로 정의되었으며, 대조군과 처리된 군을 나타냈다. 조정 p-값(FDR)이 0.05이고 절대 log₂ 폴드 변화가 1<≥ 있는 유전자는 유의하게 차이 있게 발현된 것으로 간주되었습니다. 효과 크기 추정의 견고성을 높이기 위해 apeglm 방법을 사용하여 Log₂ 폴드 변화 축소를 수행하였습니다. 분석 결과는 EnhancedVolcano (v1.22.0)¹⁵ 와 ggplot2¹⁶을 사용하여 시각화되었으며, 이들은 발현 규모와 통계적 신뢰도 간의 관계를 보여주는 화산 플롯과 MA 플롯을 생성했습니다. DESeq2 내에서 확산 추정치도 평가하여 생물학적 복제체 간 정확한 분산 모델링과 일관된 정규화를 보장하였다¹⁷.

미토콘드리아 산화 스트레스 관련 차별 발현 유전자(MOS-DEGs)의 검색 및 확인

에너지 대사, 산화 스트레스, 약물 내성 간의 연관성을 조사하기 위해, Human MitoCarta3.0¹⁸ (https://personal.broadinstitute.org/scalvo/MitoCarta3.0/human.mitocarta3.0.html), 유전자 온톨로지(GO:0006979, 산화 스트레스 반응)(http://geneontology.org/), 교토 유전자 및 유전체 백과사전(KEGG) 산화 인산화 경로 등 여러 데이터베이스에서 미토콘드리아 및 산화 스트레스 관련 유전자 목록을 종합적으로 수집했습니다 (https://www.genome.jp/kegg/), 인간 산화 스트레스 유전자 데이터베이스(HOSGDB)(http://hosgdb.com/). 모든 회출된 유전자는 org를 사용하여 HGNC 승인 유전자 심볼로 표준화되었습니다. Hs.eg.db(v3.18.0)와 주석Dbi(v1.64.0)는 중복된 항목, 의사유전자 및 비암호화 RNA는 주석 정확성을 보장하기 위해 제거되었습니다. 그 결과 선별된 미토콘드리아 산화 스트레스 유전자 패널(MOS 유전자)은 두 전사체 데이터셋에서 확인된 차등 발현 유전자와 통합하는 참고 자료로 사용되었습니다.

선별된 MOS 유전자 목록과 GSE231524 및 GSE231525에서 얻은 DEG의 교차는 r에서 dplyr(v1.1.3)¹⁹ 와 기본 R intersect() 함수를 사용해 수행되었습니다. 이러한 통합적 접근법은 미토콘드리아 대사, 산화환원 조절, 산화 스트레스 적응과 기능적으로 연관된 유전자를 대표하는 MOS-DEGs의 확인을 가능하게 했습니다. 데이터셋 간 겹침은 R의 VennDiagram(v1.7.3) 패키지를 사용하여^{실험 모델 전}반에 걸쳐 공유되고 고유한 유전자를 시각화했습니다. 정제된 MOS-DEG 목록은 후속 분석에 사용되어 HER2 표적 치료 내성의 전사 및 대사 재프로그래밍에 대한 메커니즘적 통찰을 제공했습니다.

MOS-DEG의 표현 프로파일링 및 시각화

확인된 MOS-DEGs의 발현 프로파일링은 RStudio의 ComplexHeatmap (v2.18.0)²¹ 및 pheatmap (v1.0.12) 패키지를 사용하여 부모, 약물 내성 및 내성 조건에 걸친 전반적 발현 패턴을 시각화하였습니다. 정규화된 개체 수는 z-점수 스케일링을 사용하여 샘플 간 유전자 발현 행렬을 표준화하였습니다. 클러스터링은 유클리드 거리 및 완전 연관 방법을 사용하여 공동발현 패턴을 검출하고 조건 특이적 전사 프로필을 구별하기 위해 수행되었습니다. ggplot2를 사용하여 히트맵과 군집 배치 플롯이 생성되어 조건 간 시각적 명확한 구분을 보장하였습니다. 이 시각화 접근법은 미토콘드리아 활동, 산화 스트레스 조절, 약물 저항성 상태에서의 대사 재프로그래밍과 관련된 유전자 그룹의 식별을 용이하게 했습니다.

기능 풍부화 및 경로 주석

확인된 MOS-DEGs의 생물학적 중요성과 조절 메커니즘을 탐구하기 위해, 유전자 온톨로지(GO)와 교토 유전자 및 유전체 백과사전(KEGG) 풍부화 분석이 R Studio(버전 4.3.1)를 사용하여 수행되었습니다. 분석은 여러 바이오컨덕터 패키지를 활용하여 재현 가능한 계산과 시각화를 위해 tidyverse 환경에서 수행되었습니다. 유전자 주석 및 식별자 매핑은 이 조직을 사용하여 수행되었습니다. Hs.eg.db 호모 사피엔스 참조 게놈(GRCh38)을 기반으로 한 데이터베이스(https://bioconductor.org/packages/org.Hs.eg.db/). GO 농축 분석은 clusterProfiler 패키지(버전 4.8.1; https://bioconductor.org/packages/clusterProfiler/)를 사용하여 수행되었으며, 이 패키지는 유전자를 생물학적 과정(BP), 세포 구성 요소(CC), 분자 기능(MF)의 세 가지 주요 온톨로지로 분류합니다. enrichGO²² 함수는 매개변수를 p-값을 0.05로 설정하고 p-값을 조정한 < 사용되었습니다. (FDR) < 0.05, 벤자미니-호흐버그 보정법을 적용했습니다. 막대 플롯, 점플롯, 코드 다이어그램 등 시각화는 enrichplot(https://bioconductor.org/packages/enrichplot/), ggplot2¹³ (https://cran.r-project.org/web/packages/ggplot2/), GOplot(https://cran.r-project.org/web/packages/GOplot/)을 사용해 생성되었습니다. 이 도구들은 풍부화된 GO 용어와 그 유전자 연관성에 대한 구조화된 관점을 제공했습니다.

KEGG 경로 풍부화는 clusterProfiler 패키지 내 enrichKEGG() 함수를 사용하여 수행되었으며, 이는 인간 KEGG 데이터베이스(https://www.genome.jp/kegg/)를 참조합니다. KEGGREST 패키지(https://bioconductor.org/packages/KEGGREST/)는 경로 데이터 검색 및 주석 작성에 사용되었습니다. 조정된 p-값을 가진 경로. (q-값) < 0.05는 유의한 것으로 간주되었습니다. 시각화와 경로 매핑은 pathview(https://bioconductor.org/packages/pathview/), ggplot2, enrichplot을 사용하여 수행되었으며, 네트워크 표현에는 igraph와 ggraph가 사용^{되었습니다.} 모든 농축 분석 및 시각화는 R Studio(v4.3.1)에서 재현 가능한 코드와 표준화된 Bioconductor 워크플로우를 사용하여 구현되어, MOS-DEGs와 관련된 농축 기능 범주와 생물학적 경로를 신뢰성 있게 식별할 수 있었습니다.

ROC 기반 유방암 예측 바이오마커 검증

MOS-DEG의 임상 예측력을 검증하기 위해 ROCplotter 온라인 도구(https://www.rocplot.org/)²⁴를 사용하여 수신자 작동 특성(ROC) 곡선 분석이 수행되었습니다. ROCplotter는 화학요법, 호르몬 요법, 항HER2 치료를 받은 3,104명의 유방암 환자의 유전자 발현 데이터와 임상 주석이 달린 치료 반응 데이터셋을 결합한 통합 웹 기반 플랫폼입니다.

분석은 결과 변수로 '병리학적 완전 반응'을, 치료 범주로 '모든 화학요법'을 사용하여 수행되었습니다. Affymetrix 마이크로어레이 데이터셋에서 파생된 유전자 발현 값은 플랫폼 내 임상 주석을 바탕으로 자동으로 반응자와 비응답자 그룹으로 층화되었습니다.

수신자 작동 특성(ROC) 곡선(AUC), Mann–Whitney U 검정, 접힘 변화, 카이제곱 검정을 적용하여 각 유전자가 반응자와 비반응자를 구별하는 능력을 평가하였습니다. 곡선 아래 면적(AUC)이 구별적 성과 평가의 주요 지표로 사용되었습니다. AUC 값이 0.55를 초과하고 ROC p값이 0.05< 유의미하게 간주되었으며, 이는 전사체 바이오마커에서 전형적인 중간 정도의 예측 성능을 나타내며, 분석적 엄밀성을 유지하기 위해 거짓 발견률(FDR) 보정이 적용되었습니다.

선택된 모든 MOS-DEG는 해당 Affymetrix 프로브 ID로 조회되었습니다. 각 유전자의 식별 잠재력은 임상 유방암 코호트에서 평가되었으며, 발현 데이터는 반응자와 비반응자군으로 계층화되었습니다. ROC 곡선, 박스플롯 및 관련 통계 출력은 ROCplotter 플랫폼에서 직접 생성되어 하위 시각화 및 비교를 위해 내보냈습니다. 분석은 미토콘드리아 산화 스트레스에 관여하는 산화환원 및 대사 조절 인자의 예측 가치를 정량화하였습니다. 임상 샘플 전반에 걸쳐 일관된 예측 유의성을 보인 유전자는 최종 예측 패널에 포함되도록 유지되었습니다.

종양, 정상 및 전이 조직에서의 차별 발현 분석 (TNMplot 분석)

정상, 종양, 전이 유방 조직에서 정상 유방 조직의 발현 패턴을 TNMplot 웹 도구 v2 (https://tnmplot.com/analysis/)²⁵를 사용해 분석하였다. RNA-Seq(TCGA + GTEx + MET500) 및 유전자 칩 데이터셋 모두를 검토하여 크로스 플랫폼 검증을 보장하였습니다. "다중 유전자 분석" 모듈은 유방 침습성 암과 함께 선정된 조직 유형²⁶으로 사용되었습니다. 발현 값은 log₂ 변환하여 종양 vs. 정상(TvsN), 전이성 대 종양(MvsT), 전이 대 정상(MvsN) 그룹 간 비교하였다. TNMplot은 Mann–Whitney U 검정을 사용하여 통계적 유의성을 평가하기 위해 폴드 변화(FC)와 p-값을 자동으로 계산했습니다. 발현 분포는 TNMplot 인터페이스에서 직접 생성된 박스플롯과 밀도 플롯으로 시각화되었으며, 녹색, 빨간색, 회색은 각각 정상, 종양, 전이 조직을 나타냅니다. 모든 수치는 결과 섹션에 통합하기 위해 고해상도로 내보내졌습니다. 이 이중 플랫폼 분석은 유방암 진행과 관련된 주요 미토콘드리아 산화환원-대사 조절인자를 견고하게 식별하고 검증할 수 있게 했습니다²⁷.

카플란-마이어 플로터를 이용한 생존 및 예후 분석

유방암에서 MOS-DEGs의 예후적 관련성을 평가하기 위해 Kaplan–Meier 온라인 플로터 도구(https://kmplot.com/analysis/)²⁸을 사용하여 생존 분석이 수행되었습니다. 이 데이터베이스는 여러 GEO, EGA, TCGA 데이터셋에서 파생된 4,900명 이상의 유방암 환자의 유전자 발현 및 생존 데이터를 통합합니다. 분석은 우선순위가 정해진 유전자에 대응하는 개별 Affymetrix 프로브 ID를 사용하여 재발 없는 생존율(RFS)을 위해 수행되었습니다: 225609_at(GSR), 201761_at(MTHFD2), 201619_at(PRDX3/AOP1), 201128_s_at(ACLY). 환자들은 중앙 발현 기준에 따라 고발현군과 저발현군으로 나뉘었고, 생존 확률은 Kaplan–Meier 방법으로 추정되었습니다. 로그 랭크 검정은 생존 곡선 간 통계적 유의성을 평가하는 데 사용되었으며, 95% 신뢰구간(CI)을 가진 위험비(HR)는 도구가 자동으로 계산했습니다. 모든 분석은 RFS 평가 지표를 사용하여 수행되었으며, 호르몬 수용체나 HER2 상태(ER, PR, HER2 = 모두)에 따른 제한은 없었다. 중복 샘플을 제거하고, 통계적 견고성을 보장하기 위해 비례 위험 가정을 검증하였습니다. 품질 관리 필터는 편향된 마이크로어레이를 제외했습니다. KM 플로터 기본 설정에 따라 다중 테스트에 대해 수동 프로브 선택이나 p-값 보정은 적용되지 않았습니다. 통계적 유의성은 p. < 0.05로 정의되었습니다. 생존 플롯을 시각화하여 고해상도로 다운로드하여 추가 해석을 위해 각 후보 MOS 유전자^29,30의 고발현자와 저발현자 간의 결과를 비교하였다.

본 분석은 HER2⁺ 유방암 샘플만을 사용하여 차별 발현 유전자(DEGs)와 허브 유전자를 식별하기 위해 시작되었습니다. 이후 HER2 상태(ER, PR, HER2 = 전부)에 제한 없이 생존 분석이 수행되어 확인된 유전자의 광범위한 예후 관련성과 일반화 가능성을 평가하였다. 이 접근법은 연구 초점을 재정의하기 위한 것이 아니라 2차 검증 단계로 사용되었습니다. 따라서 확인된 허브 유전자의 예후 함의는 신중하게 해석되며, 주요 결론은 HER2⁺ 유방암에 특이한 것으로 남아 있습니다.

MTHFD2-201(ENST00000394053.7)과 PRDX3-201(ENST00000298510.4)의 표준 전사체 서열은 Ensembl Genome Browser(https://www.ensembl.org)^31,32에서 추출되었습니다. 변이 주석 및 분류는 Ensembl 변이 효과 예측기(VEP)(https://www.ensembl.org/vep)를 사용하여 수행되었으며, 이 장치는 각 확인된 변이체에 대해 상세한 유전체 맥락, 코돈 변형, 아미노산 치환 정보를 제공했습니다. 오센스 변이(비동의어 SNP)만이 후속 분석을 위해 선택되었습니다.

병인성 예측 및 변이 우선순위 지정

각 nsSNP의 기능적 결과는 여러 계산 예측 도구를 조합하여 평가되었습니다. SIFT(https://sift.bii.a-star.edu.sg)는 아미노산 보존 평가에 적용되었으며, 점수 0.05인 변이≤ 해로운³³번으로 분류하였습니다. PolyPhen-2(http://genetics.bwh.harvard.edu/pph2)는 치환의 구조적·진화적 영향을 추정했으며, 점수 0.85≥ 추정 피해를 나타냈다³⁴. CADD(https://cadd.gs.washington.edu)는 여러 주석을 통합한 복합 해악성 점수를 제공했으며, ≥ 20은 높은 병원성 잠재력³⁵를 의미하였다. MetaLR³⁶, Mutation Assessor, REVEL의 보완 지표가 VEP 인터페이스에서 통합되어 예측 신뢰도를 향상시켰습니다³⁷. MetaLR ≥ 0.70, 돌연변이 평가자 ≥ 3.5, REVEL ≥ 0.75 임계값을 충족하는 변이가 병원성 가능성이 높은 것으로 우선순위를 보였다.

구조적 및 기전적 영향 예측

아미노산 치환이 구조적 완전성과 생화학적 기능에 미치는 영향을 평가하기 위해, 상위 순위의 nsSNP는 MutPred2 (http://mutpred.mutdb.org)³⁸ 과 DynaMut (http://biosig.unimelb.edu.au/dynamut)³⁹를 사용해 추가로 분석하였습니다. MutPred2는 촉매 활성 변화, 금속 결합 잔기의 증가 또는 손실, 용매 접근성 변화, 알로스테릭 조절 등의 기능 장애 확률을 추정했으며, ≥0.80 점수는 매우 병원성으로 분류되었습니다. DynaMut은 야생형 단백질과 돌연변이 단백질 간의 깁스 자유 에너지 변화(ΔΔG)를 계산하여 안정성 변화의 방향과 크기를 평가하고, 원자 변위와 수소 결합 재배열의 시각화를 생성했습니다.

2차 및 3D 구조 모델링과 용매 접근성 프로파일링

MTHFD2와 PRDX3의 실험적으로 분리된 결정 구조는 단백질 데이터 뱅크(PDB)에서 얻어져 PyMOL V:3.1 (https://pymol.org)⁴⁰ 을 사용해 유해 잔기의 공간적 분포를 시각화하였습니다. 돌연변이 모델은 해당 아미노산 치환을 도입한 후 구조적 정제와 에너지 최소화를 통해 만들어졌습니다. 3D 비교 검사에서는 2차 원소의 변화, 원자 간 접촉 변화, nsSNP가 촉매 및 보동인자 결합 도메인과 공간적으로 가까이 있다는 점이 강조되어 산화환원 및 대사 기능에 잠재적 교란이 있음을 드러냈습니다.

2차 구조 및 용매 노출 분석은 PSIPRED V: 3.2 (http://bioinf.cs.ucl.ac.uk/psipred)^41,42와 NetSurfP 3.0 (https://services.healthtech.dtu.dk/service.php?NetSurfP-2.0)⁴³을 사용하여 수행되었습니다. 이 도구들은 α나선, β가닥, 코일, 무질서한 영역과 상대용매 접근성(RSA) 점수를 예측했습니다. 중간에서 높은 RSA 값과 구조 질서를 보이는 잔류물을 지도화하여 용매 노출 및 기능적 임계 위치를 식별하였습니다. 영향을 받은 부위들은 2D 위상수도에 시각화되어 강체 촉매 코어에서 해로운 돌연변이가 발생했는지 또는 유연한 루프 영역에서 발생했는지 확인하여, 단백질 접힘 역학과 효소 효율성에 미칠 확연한 영향을 예측하였다.

데이터베이스 교차 검증, 기능 통합 및 안정성 검증

우선순위가 정한 각 nsSNP는 dbSNP, 1000 Genomes, ExAC, gnomAD 등 집단 수준 유전체 데이터베이스와 교차 참조하여 변이 빈도, 전반적 대립유전자 분포, 이전에 보고된 임상 연관성을 확인하였습니다. 진화 보존, 구조 모델링, 머신러닝 기반 기능 예측의 통합을 통해 MTHFD2와 PRDX3에서 높은 신뢰도의 해로운 변이를 식별할 수 있었습니다. 이러한 고영향 돌연변이들은 이후 기능 영역에 매핑되어 미토콘드리아 산화 스트레스 불균형, 대사 신호 전달 변화, 유방암 치료 내성에 대한 잠재적 역할을 밝혀냈습니다. 각 해로운 치환의 열역학적 결과를 추가로 확인하기 위해, iMutant 3.0 (https://folding.biofold.org/i-mutant/i-mutant3.0.html)⁴⁴ 를 사용하여 염기서열 및 구조 데이터를 사용하여 돌연변이가 단백질 안정성에 미치는 영향을 예측하였습니다. 분석 결과 야생형 단백질과 돌연변이 단백질 간 자유 에너지 변화를 나타내는 ΔΔG 값(kcal/mol)이 계산되었습니다. 음의 ΔΔG 값을 보이는 변이체는 불안정화 돌연변이로 분류되어 단백질 안정성 감소와 전개 확률 증가를 나타냈다. iMutant 예측과 DynaMut 및 MutPred2 결과를 통합하여 산화환원 기능, 촉매 무결성 및 전반적인 단백질 구조 안정성에 영향을 미칠 가능성이 있는 구조적으로 중요한 잔기의 상호 검증을 제공하였습니다.

통합 기능적 해석과 치료적 관련성

확인된 모든 해로운 nsSNP는 dbSNP, gnomAD, ExAC 집단 데이터베이스와의 교차 검증을 통해 소규모 대립유전자 빈도와 이전에 보고된 암 표현형과의 연관성을 검증하였습니다. 진화 보존, 구조 모델링, 안정성 데이터의 통합 해석에 따르면, 고영향 돌연변이 rs1471336772(MTHFD2)와 rs747786383(PRDX3)이 단백질 구조와 촉매 효율에 가장 강한 해로운 영향을 미친다는 것을 보여주었습니다. 계산 결과는 MTHFD2의 돌연변이가 NADPH 의존적 산화환원 대사를 불안정하게 하는 반면, PRDX3의 돌연변이는 과산화효소 매개 산화 스트레스 방어를 저해하여 미토콘드리아 기능 장애와 종양 공격성에 기여함을 시사합니다. 이 nsSNP 기반 구조 및 기능 분석은 향후 치료적 선별과 돌연변이 검증을 위한 계산적 기반을 제공하며, MTHFD2와 PRDX3가 산화환원 표적 유방암 치료의 정밀 바이오마커임을 강조합니다. 분석 전략의 명확성을 높이고 포괄적인 개요를 제공하기 위해, 연구의 주요 단계를 요약한 도식적 워크플로우가 그림 2에 제시되어 있습니다. 이 워크플로우는 차등 유전자 발현 분석, 미토콘드리아 유전자 필터링, 단백질-단백질 상호작용 네트워크 구축, ROC 분석을 이용한 임상 검증, 그리고 nsSNP 기반 구조 특성 분석을 통합합니다. 이 단계별 프레임워크는 전사체 데이터 처리에서 바이오마커 식별 및 기능 해석으로의 논리적 진행 과정을 강조합니다.

그림 2. HER2⁺ 유방암에서 미토콘드리아 산화 스트레스 관련 바이오마커의 식별 및 검증을 위한 통합 다단계 워크플로우. 이 도 식은 연구에서 사용된 분석 파이프라인을 요약한 것입니다. 첫째, RNA-seq 데이터셋(GSE231524 및 GSE231525)에서 차등 유전자 발현(DEG) 분석을 수행하여 유의하게 변한 유전자를 식별하였습니다. 이 DEGs는 선별된 미토콘드리아 산화 스트레스 관련 유전자와 교차하여 MOS-DEG를 얻었습니다. 다음으로, STRING과 Cytoscape를 사용하여 단백질-단백질 상호작용(PPI) 네트워크 분석을 수행하여 허브 유전자와 기능 모듈을 식별하였습니다. 이후 ROCplotter 플랫폼을 사용하여 임상 코호트에서 선택된 유전자의 예측 성능을 평가하기 위해 수신자 작동 특성(ROC) 곡선 분석을 적용하였습니다. 마지막으로, 비동의어 SNP(nsSNP) 분석 및 구조 모델링을 수행하여 우선순위 유전자(MTHFD2 및 PRDX3)에서 주요 변이의 잠재적 기능적·구조적 영향을 평가하였습니다. 이 통합 워크플로우는 전사체, 네트워크, 임상, 구조 분석을 연결하여 잠재적 바이오마커와 치료 표적을 식별합니다. 이 그림의 더 큰 버전을 보시려면 여기를 클릭해 주세요.

Results

차등 유전자 발현 분석

HER2⁺ 유방암 진행 및 치료 내성과 관련된 전사 변이를 조사하기 위해 차별 발현 프로파일링을 수행했으며, 두 개의 독립적인 RNA-seq 데이터셋인 GSE231524와 GSE231525. 첫 번째 데이터셋(GSE231524)에서는 총 19,727개의 유전자가 처음에 정량되었습니다. 정규화, 필터링, 조정 회귀 모델(ARM) 적용 후 차별 발현에 대해 6,604개의 유전자가 유의미하게 발현된 상태로 유지되었습니다(보충 표 S1). 유전자 기호 표준화, 중복 제거 및 주석 정제를 거쳐 6,300개의 고유 유전자가 후속 분석을 위해 최종 확정되었습니다. 화산 플롯, MA 플롯, 분산 추정치(그림 3A-C)는 유전적으로 상향 및 하향 조절된 유전자의 분포를 명확히 보여줍니다(조정 p < 0.05, |log₂FC| ≥ 1)....

Discussion

미토콘드리아 산화 스트레스는 세포 생존, 대사 가소성, 치료 저항성을 관장하는 종양 적응의 결정적 특징입니다. 본 연구는 전사체 데이터, 기능 풍부화, 임상 검증, 실리코 돌연변이 및 구조 분석을 통합하여 HER2⁺ 유방암에서 미토콘드리아 산화 스트레스 관련 차별 발현 유전자(MOS-DEGs)를 체계적으로 탐구하였다. 차별 발현과 돌연변이 프로파일링의 통합은 산화환원 대사 조절이 유전체 불안정성과 어떻게 교차하는지에 대한 다차원적 이해를 제공했습니다. 확인된 MOS-DEG 중 MTHFD2와 PRDX3는 산화 스트레스 적응, 대사 조절, 치료 반응에서 서로 다른 상호보완적 기능을 나타내는 중추적 미토콘드리아 결정인자로 부상했다⁴⁸. 이 유전자 내에서 여러 해로운 비동의 SNP(nsSNP)가 확인되어, 단일 아미노산 치환이 유방암 세포에서 효소 활성, 단백질 접힘, 미토콘드리아 회복력을 조절할 ...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
AnnotationDbi	Bioconductor	v1.64.0	R 주석 패키지로 HGNC 유전자 기호 표준화 및 주석에 사용됩니다.
apeglm	Bioconductor	DESeq2 방법	차등 발현 분석에서 log2 폴드 변화 축소에 사용되는 방법입니다.
Biobase	Bioconductor	v2.62.0	GEO에서 파생된 발현 및 메타데이터 객체에 액세스하고 관리하는 데 사용되는 R 패키지입니다.
CADD	University of Washington / Kircher Lab	웹 도구	nsSNP 병원성 예측에 사용되는 결합 주석 종속적 고갈 점수입니다.
clusterProfiler	Bioconductor	v4.8.1	GO 및 KEGG 풍부화 분석에 사용되는 R 패키지입니다.
ComplexHeatmap	Bioconductor	v2.18.0	히트맵 생성 및 발현 클러스터링에 사용되는 R 패키지입니다.
Cytoscape	Cytoscape Consortium	소프트웨어	STRING에서 파생된 PPI 네트워크 및 CytoHubba 우선 순위 지정과 함께 사용되는 네트워크 시각화 플랫폼입니다.
dbSNP	NCBI	데이터베이스	우선 순위가 지정된 nsSNP의 교차 검증에 사용되는 인구 변이 데이터베이스입니다.
DESeq2	Bioconductor	v1.42.0	정규화, 분산 모델링 및 차등 발현 분석에 사용되는 R 패키지입니다.
dplyr	CRAN / tidyverse	v1.1.3	데이터 조작 및 유전자 세트의 교차에 사용되는 R 패키지입니다.
DynaMut	University of Melbourne / BioSig	웹 서버	단백질의 돌연변이와 관련된 안정성 및 유연성 변화를 추정하는 데 사용되는 도구입니다.
EnhancedVolcano	Bioconductor	v1.22.0	차등 발현 결과를 화산 플롯으로 시각화하는 데 사용되는 R 패키지입니다.
Ensembl Genome Browser	EMBL-EBI / Ensembl	데이터베이스	MTHFD2 및 PRDX3에 대한 표준 전사체 서열의 출처입니다.
Ensembl Variant Effect Predictor (VEP)	EMBL-EBI / Ensembl	웹 도구	미스센스 변이를 주석화하고 예측 점수를 통합하는 데 사용되는 도구입니다.
ExAC	Broad Institute	데이터베이스	nsSNP 교차 검증에 사용되는 인구 수준의 엑솜 데이터베이스입니다.
Gene Expression Omnibus (GEO)	NCBI	데이터베이스	HER2 양성 유방암 분석을 위한 전사체학 데이터세트를 검색하는 데 사용되는 저장소입니다.
Gene Ontology	Gene Ontology Consortium	GO:0006979	산화 스트레스 관련 유전자를 검색하고 풍부화 분석에 사용되는 온톨로지 리소스입니다.
GEOquery	Bioconductor	v2.70.0	GEO 카운트 행렬 및 메타데이터에 액세스하는 데 사용되는 R 패키지입니다.
ggplot2	CRAN / tidyverse	v3.5.0	데이터 시각화, 클러스터링 플롯 및 그래픽 출력에 사용되는 R 패키지입니다.
ggraph	CRAN	R 패키지	풍부화 및 경로 표현 중 그래프 및 네트워크 시각화에 사용되는 R 패키지입니다.
gnomAD	Broad Institute	데이터베이스	우선 순위가 지정된 nsSNP의 빈도 및 분포를 확인하는 데 사용되는 인구 변이 데이터베이스입니다.
GOplot	CRAN	R 패키지	코드 및 요약 플롯을 포함하여 GO 풍부화 시각화에 사용되는 R 패키지입니다.
GSE231524	NCBI GEO	액세스	부모, 약물 내성 및 저항성 BT474 표현형 분석에 사용되는 RNA-seq 데이터세트입니다.
GSE231525	NCBI GEO	액세스	HER2 양성 유방암 세포에서 DUSP6 결손 분석에 사용되는 RNA-seq 데이터세트입니다.
Human MitoCarta3.0	Broad Institute	데이터베이스	미토콘드리아 유전자 집합을 정의하는 데 사용되는 정리된 미토콘드리아 유전자 리소스입니다.
Human Oxidative Stress Gene Database (HOSGDB)	HOSGDB	데이터베이스	산화 스트레스 관련 유전자를 검색하는 데 사용되는 데이터베이스입니다.
igraph	CRAN	R 패키지	네트워크 표현 및 경로 시각화에 사용되는 R 패키지입니다.
iMutant 3.0	University of Bologna / Biofold	웹 서버	단백질 안정성에 대한 돌연변이 효과를 예측하는 데 사용되는 도구입니다.
Kaplan–Meier Plotter	KMplot	웹 도