$$\rightleftharpoonup{xx}$$
$$\longleftharp{xx}$$,
$$\longrightharp{xx}$$,
이 연구는 유전자 발현 옴니버스(GEO) 데이터베이스에서 공개된 식별 해제 데이터셋만을 사용했습니다. 이 작업은 기존 공개 데이터의 2차 분석이었고 직접 참여자 접촉, 개입, 식별 가능한 개인 정보 접근이 포함되지 않았기 때문에 추가적인 윤리 위원회 승인과 사전 동의가 필요하지 않았습니다.
데이터 소스 및 전처리
모든 유전자 발현 및 단일 세포 데이터셋은 GEO 데이터베이스에서 얻었습니다24. 주요 우울장애의 경우, 128명의 환자와 64명의 건강한 대조군의 말초혈액 샘플을 포함하는 데이터셋 GSE98793가 사용되었습니다. 피부근염의 경우, 호모 사피엔스 발현 프로파일링, 명확히 식별 가능한 질병 및 대조군, 프로브-유전자 매핑을 위한 플랫폼 주석, 발견 또는 검증 분석 적합성 등 사전 정의된 기준에 따라 데이터셋이 선정되었습니다. GEO 시리즈에 여러 염증성 근병증 아형이 포함된 경우, 본 연구에서는 피부근염과 정상 대조군 샘플만 추출하였다. GSE1551, GSE46239, GSE128470가 발견/훈련 데이터셋으로 사용되었고, GSE5370, GSE39454, GSE11971가 독립적인 검증 데이터셋으로 사용되었습니다. 이 연구에서 분석된 피부근염 데이터셋은 주로 말초혈보다는 영향을 받은 근육이나 피부 조직에서 도출되었습니다. 피부근염의 단세포 데이터는 데이터셋 GSE190510에서 수집되었습니다.
GEO 데이터베이스에서 원시 표현 행렬과 해당 플랫폼 주석 파일이 다운로드되었습니다. 프로브 ID는 제조사에서 제공한 GPL 주석에 따라 공식 유전자 기호에 매핑되었습니다. 단일 공식 유전자 기호에 명확하게 매핑할 수 없는 프로브는 제거되었다. 여러 탐침이 동일한 유전자에 매핑될 때, limma 패키지의 'avereps' 함수가 구현한 평균 발현값을 사용하여 유전자 수준에서 압축하여 유전자별 샘플별 발현 행렬을 생성했습니다.
강도 의존 편향을 줄이고 분산을 안정화하기 위해, 표현 값 분포에 따라 적절한 경우 log2 변환이 적용되었습니다. 그 후 limma 패키지의 'normalizeBetweenArrays' 함수를 사용하여 배열 사이 정규화가 수행되었습니다. 누락된 값이 있을 경우 K-최근접 이웃 보충을 사용해 보충했습니다. 통합 피부근염 훈련 데이터셋의 경우, sva 패키지의 'ComBat' 기능을 사용하여 배치 교정이 수행되었으며, 데이터셋/플랫폼 출처를 배치 변수로 간주하고 샘플 그룹(피부근염 대 건강한 대조군)을 설계 행렬에 포함시켜 배치 조정 중 관심 있는 생물학적 변이를 보존했습니다.
모든 분석은 데스크톱 운영체제에서 R을 위한 통합 개발 환경을 사용하여 R으로 수행되었습니다. 림마 패키지는 프로브 요약 및 정규화에 사용되었습니다. SVA 패키지는 ComBat 배치 교정에 사용되었습니다. 누락된 값은 k = 10인 K-최근접 이웃 보정을 사용해 보정되었습니다.
가중 유전자 공동 발현 네트워크 분석
가중 유전자 공동 발현 네트워크 분석(WGCNA)은 WGCNA R 패키지25,26을 사용하여 주요 우울장애 및 피부근염 데이터셋에 대해 별도로 수행되었습니다. 샘플은 flashClust를 사용해 계층적으로 군집화하여 이상치를 식별했습니다; 수목도 높이가 100을 초과하고 분산 하위 25%에 속하는 유전자는 제외하였습니다. 각 네트워크에 대해 pickSoftThreshold를 사용하여 소프트 임계값 파워(β)를 선택하여 근사적인 스케일 프리 토폴로지(R2 > 0.8)를 달성했습니다. 인접 행렬은 위상 중첩 행렬(TOM)으로 변환되었고, 모듈은 최소 모듈 크기가 60이고 병합 절단 높이가 0.2527인 동적 트리 커팅을 통해 식별되었습니다. WGCNA R 패키지는 계층적 클러스터링을 위해 flashClust와 함께 사용되었습니다. 무작위 시드는 재현성을 위해 12345로 설정되었습니다. 모듈 고유유전자는 Pearson 상관관계를 사용하여 질병 상태와 상관관계 관계가 부여되었으며, P-값을 벤자미니-호흐버그 방법으로 조정했습니다. 각 질병에 대해 질병 상태와 가장 강하고 중요한 연관성을 보이는 모듈이 핵심 질병 관련 모듈로 유지되었습니다. 주요 우울장애 데이터셋과 피부근염 데이터셋의 핵심 모듈 유전자 간 중복이 하위 분석을 위한 후보 공유 유전자 집합으로 정의되었습니다. 통합 피부근염 코호트의 차별 발현 분석은 별도로 수행되어 피부근염 관련 전사 변화를 특성화하였습니다.
기능 풍부 분석
유전자 온톨로지(GO) 풍부화 분석은 R을 사용하여 수행되었으며, 유전자 심볼은 org를 사용하여 Entrez ID로 변환되었습니다. Hs.eg.db 그리고 clusterProfiler에서 enrichGO를 사용해 유의미하게 풍부해진 GO 항(p < 0.05)이 확인되었습니다. 결과를 다차원적으로 시각화하기 위해 enrichplot 패키지를 사용해 막대 플롯과 버블 플롯을 생성했고, Circlize 패키지를 사용해 GO 범주, 유전자 개수, 풍부도 인자를 표시하는 원형 플롯을 구성했습니다. ComplexHeatmap 패키지와 함께 전설이 추가되었습니다. 교토 유전자 및 유전체 백과사전(KEGG)의 차별 발현 유전자에 대한 경로 풍부화 분석도 R에서 수행되었으며, 유전자 심볼은 조직을 기반으로 Entrez ID로 변환되었습니다. 데이터베이스Hs.eg.db 그리고 유의미하게 풍부해진 경로(FDR < 0.05)는 clusterProfiler 패키지 28,29,30,31의 enrichKEGG 함수를 사용해 확인되었습니다. 농축 결과는 막대 및 버블 플롯을 사용하여 시각화하였습니다.
GeneMANIA 기반 기능 연합 네트워크 분석
이전에 확인된 공유 유전자를 바탕으로, GeneMANIA 기반 기능적 연관 네트워크가 구축되어 이 유전자들과 관련 파트너 간의 상호작용 맥락을 탐구하였습니다. 유전자 목록은 호모 사피엔스를 기준종으로 GeneMANIA에 제출했습니다. GeneMANIA는 공발현, 물리적 상호작용, 경로, 공동위치, 유전적 상호작용, 공유 단백질 도메인 등 여러 근거 유형을 통합합니다. 생성된 네트워크는 시각화 및 분석을 위해 네트워크 시각화 플랫폼으로 내보내고 가져왔습니다. 네트워크 시각화 플랫폼에서 위상적 분석을 수행하여 고도로 연결된 후보 노드32, 33, 34를 식별하기 위한 시각화 및 분석을 수행했습니다.
머신러닝 기반 진단 모델 구축
진단 분류에는 랜덤 포레스트(RF), 서포트 벡터 머신(SVM), 선형 판별 분석(LDA), 나이브 베이즈, 그라디언트 부스트 머신(GBM), XGBoost, glmBoost, 탄성 넷(Enet), 능선, 최소 절대 축소 및 선택 연산자(LASSO), 단계별 일반화 선형 모델(Stepglm), 편소 최소제곱 회귀 일반화 선형 모델(plsRglm) 등 여러 머신러닝 알고리즘이 사용되었습니다.. 2단계 모델링 프레임워크를 적용하여 113개의 후보 모델 조합을 생성하였습니다. 1단계에서는 초기 알고리즘이 훈련 코호트에서 가변 선별에 사용되었고; 두 번째 단계에서는 보존된 변수들을 진단 분류 모델에 적합시키는 데 사용되었습니다. ≤ 5개의 선택된 변수를 가진 모델은 추가 비교에서 제외되었습니다. 결합된 피부근염 데이터셋이 훈련 코호트로 사용되었으며, 라벨은 피부근염과 건강한 대조군으로 구분되었고, 독립적인 검증 코호트는 외부 성과 평가에 사용되었습니다. 내부 재샘플링과 튜닝은 알고리즘별로 이루어졌으며, glmnet 기반 모델(LASSO, Ridge, Elastic Net)은 lambda.min을 선택하기 위해 10배 교차 검증을 사용했습니다; GBM은 최적의 트리 수를 결정하기 위해 10배 내부 교차 검증을 사용했습니다; XGBoost는 최소 테스트 로그 손실에 따라 최종 부스터팅 라운드를 5배 재샘플링으로 선택했으며; glmBoost는 cvrisk 기반 내부 교차 검증을 사용하여 종료 반복을 결정했습니다; LDA는 캐릿 교차 검증 프레임워크 하에 적용되었습니다. 현재 구현에 명시적인 튜닝 단계가 없는 알고리즘의 경우, 고정 또는 패키지 기본 설정이 사용되었습니다. 정보 유출을 줄이기 위해 특징 선택, 모델 적합, 내부 조정은 훈련 코호트만을 사용하여 수행되었으며, 검증 코호트는 독립적인 예측 및 AUC 기반 성과 평가에만 사용되었습니다. caret 패키지는 glmnet, randomForest, e1071, gbm, xgboost, mboost, plsRglm , MASS 등 개별 알고리즘을 위한 머신러닝 워크플로우 관리에 사용되었습니다. SHAP 분석은 shapviz 패키지를 사용하여 수행되었습니다. 무작위 시드는 각 모델 피팅 전에 12345로 설정되었습니다. 선택된 특징이 5개 미만인 모델은 제외되었습니다. SHapley Additive exPlanations(SHAP)를 사용하여 모델 해석 가능성과 유전자 수준 기여도를 추가로 평가했으며, 가장 유익한 유전자가 하위 생물학적 해석을 위한 모델 선택 특징 후보로 우선순위를 정했습니다.
진단 성과 평가
후보 바이오마커의 진단 성능을 평가하기 위해 "pROC" R 패키지를 사용하여 수신자 작동 특성(ROC) 곡선을 생성하였습니다. 후보 마커의 발현 수준과 예측 정확도는 독립적인 데이터셋(GSE5370, GSE11971, GSE39454)에서 검증되었습니다. 모델 성능은 혼동 행렬을 사용하여 추가로 평가하였습니다. 주요 모듈 유전자의 차등 발현은 화산 및 박스 플롯을 사용해 시각화되었고, 개별 유전자의 진단 가치를 평가하기 위해 ROC 곡선을 구축하였습니다.
유전자 집합 풍부화 분석
후보 공유 전사체 신호와 관련된 조정된 기능 변화를 탐구하기 위해 clusterProfiler36,37을 사용하여 유전자 집합 풍부 분석(GSEA)을 수행하였습니다. 피부근염과 대조군 샘플의 유전자 발현 데이터는 차등 발현에 따라 순위가 매겨졌습니다. KEGG 경로(c2.cp.kegg.Hs.symbols.gmt)에 해당하는 미리 정의된 유전자 집합을 사용하여 각 경로 내 유전자가 상향 또는 하향 조절의 조정된 경향을 보이는지 평가하였습니다. 통계적 유의성은 P < 0.05로 정의되었습니다.
면역 세포 침투 분석
정규화되고 log2 변환되었으며 배치 교정된 피부근염 매트릭스가 면역 디컨볼루션에 사용되었습니다. CIBERSORT 알고리즘을 적용하여 LM22 참조 행렬38을 사용하여 면역 세포 아형의 상대적 풍부도를 추정하였습니다. 디컨볼루션 P 가 0.05< 샘플은 하위 분석을 위해 보관되었습니다. 그룹 간 추정된 면역세포 비율 차이는 박스플롯을 사용하여 시각화되었고, 면역세포 하위 집합과 후보 공유 유전자 간의 연관성을 평가하기 위해 스피어먼 상관분석을 수행하였다.
세포 맥락화를 위한 단일 세포 RNA 시퀀싱 분석
단세포 RNA-seq 분석은 R에서 Seurat를 사용해 수행되었습니다. 배치 보정에는 Harmony, 더블릿 검출에는 DoubletFinder, 주변 RNA 추정에는 celda/decontX, 의사시간 궤적 분석에는 Monocle, 세포 간 통신 분석에는 CellChat, 유전자 집합 활성 점수는 AUCell, ssGSEA 점수는 GSVA가 사용되었습니다. 원시 카운트 행렬은 min.cells = 5, min.features = 300의 매개변수를 가진 Seurat 객체에 가져왔습니다. 각 세포별로 미토콘드리아, 리보솜, 헤모글로빈 유전자 비율을 포함한 품질 관리 지표가 계산되었습니다. 세포는 다음 기준을 모두 충족할 때만 보유했습니다: nFeature_RNA > 500, nCount_RNA < 5,000, percent_mito < 25, percent_ribo > 3, percent_hb < 1. 3개 미만의 세포에서 검출된 유전자는 제외하였다. 또한, 후속 분석 전에 MALAT1과 미토콘드리아 유전자를 제거하였습니다. 초기 필터링 후, DoubletFinder를 사용해 각 샘플에서 이중항을 식별했으며, PCS = 1:30, pN = 0.25; 예상 이중복 비율은 샘플별 세포 수에 따라 설정되었습니다(<4,000 세포: 2.5%; 4,000–8,000 세포: 5%; >8,000 세포: 6.5%). 싱글렛만 남겨졌다. 주변 RNA 오염은 decontX를 사용해 추가로 추정했으며, 오염 점수 0.2< 세포는 유지되었습니다.
필터링된 데이터는 LogNormalize 방법으로 10,000의 척도 인수를 기준으로 정규화되었고, 이어서 가변 유전자 확인, 데이터 척도화, 주성분 분석이 이루어졌습니다. 샘플 간 배치 효과는 orig.ident를 배치 변수로 사용하여 Harmony를 보정했습니다. 처음 15개의 하모니 차원은 UMAP 시각화와 이웃 그래프 구성에 사용되었습니다. 클러스터링은 FindNeighbors와 FindClusters를 사용하여 수행되었으며, 최종 클러스터링 결과는 해상도 0.05로 정의되었습니다. 세포 유형은 정통 마커 유전자와 FindAllMarkers 결과39에 따라 수동으로 주석을 달았습니다.
하위 기능적 맥락화를 위해 후보 유전자 활성은 단일 세포 수준에서 평가되었고, 관련 면역 세포 하위 집합은 궤적 및 세포 간 소통 분석을 받았다. 의사시간 분석은 단안경을 사용해 DDRTree 기반 차원 감소를 사용한 후 세포 순서 분석을 수행하였습니다. 세포 간 통신 분석은 CellChat을 사용해 인간 리간드-수용체 데이터베이스를 사용했으며, 분비 신호 범주로 제한되었으며, 10개 미만의 세포가 포함된 통신은 걸러냈습니다.
각 세포에 대해 후보 유전자 활성은 AUCell, ssGSEA, AddModuleScore 세 가지 상호 보완적 접근법을 사용하여 정량화되었습니다. AUCell 점수는 유전자 순위 행렬을 기반으로 계산되었고, ssGSEA 점수는 GSVA 프레임워크를 사용해 생성되었습니다. AddModuleScore는 Seurat 내장 함수를 사용해 계산되었습니다. 결과적으로 생성된 AUCell, ssGSEA, AddModuleScore 값은 단일 점수 행렬로 결합되었습니다. 각 점수 유형은 먼저 Z-점수 변환으로 표준화되었고, 이후 최소-최대 정규화를 통해 0–1 범위로 재척도되었습니다. 각 셀의 최종 종합 점수("점수")는 세 개의 정규화된 점수의 합으로 정의되었습니다:
점수 = 정규화된 AUCell + 정규화된 ssGSEA + 정규화된 AddModuleScore.
하위 그룹 분석을 위해 CD8⁺ T세포 하위 집합을 추출하고, 이 하위 집합 내 중앙값 점수 값에 따라 세포를 이분화하였습니다. 점수 점수가 중앙값보다 높은 셀은 High_Hub_genes 그룹에, 나머지 셀은 Low_Hub_genes 그룹에 할당되었습니다.
