**게놈 전체 연관 연구(GWAS)**는 대규모 코호트에서 유전자형과 표현형 간의 관계를 연구하는 중요한 도구입니다. GWAS의 주요 한계점은 병인 변이와 관련된 연결 불균영(LD) 영역에서 생물학적 메커니즘을 식별할 수 있지만, 유전자 수준에서의 표현형 연관을 직접 얻을 수 없다는 점입니다. 이 문제를 해결하기 위해, 전사체 전체 연관 연구(TWAS) 기반 대안 방법이 유전자와 특성 간의 연결을 밝히기 위해 유전자 발현에 초점을 맞춥니다. 유전자형, 전사체학 및 표현형 데이터를 보유한 대규모 코호트에 대해 2단계 최소제곱 회귀를 적용할 수 있지만, 이 방법에도 여러 문제점이 존재합니다:
- i) 이러한 코호트를 구축하는 데 비용이 많이 듭니다;
- ii) 회귀 분석에서 유전자와 특성 간의 인과 방향을 직접 추출하기 어렵습니다.
이러한 문제점을 극복하기 위해 Porcu 등은 2단계 메델란드 무작위화 방법 TWMR를 제안했습니다. 이 방법은 GWAS 요약 통계 데이터와 공개된 eQTL 데이터를 연결하여 유전자형과 유전자 발현을 관련시키며, 메델란드 무작위화 방법 계열에 속합니다. 간단히 말해, TWMR는 유전적 변이를 도구 변수로, 유전자 발현을 노출 변수로, 관심 특성을 결과 변수로 사용합니다. 이를 바탕으로 Porcu 등은 TWMR 방법의 역방향 버전인 revTWMR도 제안했습니다. eQTL을 다른 유형의 xQTL 데이터로 대체하면 TWMR 계열 방법을 사용하여 단백체학 및 메틸화체학과 같은 다른 오믹스의 인과 효과를 탐색할 수도 있습니다.
▲ 위는 Porcu 등이 개발한 revTWMR 방법으로, 21년에 NC(Nature Communications)에 발표되었습니다. 가장 초기의 TWMR 방법은 2019년에 동일한 잡지에 발표되었습니다.
TWMR와 revTWMR 사용을 가속화하기 위해, 스위스 로잔 대학의 연구자들은 2024년 8월 10일에 **Bioinformatics [IF=4.4]**에 발표한 최적화되고 더 빠른 Python 라이브러리 **pyTWMR**를 개발했습니다. pyTWMR는 GPU 연산을 지원하면서 TWMR와 revTWMR 프로세스를 통합하며, 고도로 상관된 유전적 변이와 유전자 발현을 처리할 때 더욱 강력합니다.
▲ DOI: 10.1093/bioinformatics/btae505
본 글에서는 **pyTWMR**의 사용 예시를 소개하며, 코드 및 관련 파일은 후문에 pyTWMR라고 답변하면 얻을 수 있습니다. 관심 있는 독자들은 아래 주소에서 더 많은 정보를 확인할 수 있습니다:
1 pyTWMR 설치
.ipynb 파일이나 명령줄에서 다음 코드를 실행하여 설치합니다:
%%bash
pip3 install git+https://github.com/soreshkov/pyTWMR.git
네트워크가 불안정한 경우 로컬 설치를 권장합니다(pyTWMR-master 폴더 경로로 이동 필요):
%%bash
# Python 버전이 3.10 이상이면 성공적으로 실행 가능
pip install ./pyTWMR-master/
- 편집자는 Linux 환경에서 설정했으며, Windows 환경에서 설정 시 위의 bash를 cmd로 수정하십시오.
2 필요한 라이브러리 가져오기
다음 코드를 실행하여 필요한 라이브러리를 가져옵니다(사전 설치 필요):
from TWMR_engine import TWMR_analysis
from revTWMR_engine import revTWMR_analysis
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')
3 TWAS 분석
① TWAS 분석의 입력은 GWAS와 eQTL 데이터이며, 다음 코드를 사용하여 예시 데이터를 읽어옵니다:
genetic_effects = pd.read_csv('datasets/TWMR/ENSG00000000419.matrix', sep='\t', index_col=0)
genetic_effects.head()
▲ 코드에서 genetic_effects 표 설명: 각 행은 SNP를 나타내고, 각 열은 유전자를 나타내며, 채워진 값은 QTL 결과에서 SNP가 유전자 발현에 대한 효과이며, 마지막 열은 GWAS 결과에서 각 SNP가 연구 중인 특성에 대한 효과입니다.
또한 genetic_effects 표의 모든 SNP 간의 LD 상관 행렬이 필요하며, 위에는 총 3개의 SNP가 있으므로 행렬 형태는 3*3입니다:
linkage = pd.read_csv('datasets/TWMR/ENSG00000000419.ld', sep='\t', header=None).to_numpy().astype(np.float32)
linkage
#array([[1. , 0.0487694, 0.151668 ],
# [0.0487694, 1. , 0.118453 ],
# [0.151668 , 0.118453 , 1. ]], dtype=float32)
② TWMR 분석을 실행하고 결과를 하나의 표로 통합합니다:
gwas_sample_size = 239087 # GWAS 연구의 표본 수
qtl_sample_size = 32000 # QTL 연구의 표본 수
analysis_result = TWMR_analysis(
beta_values=genetic_effects.drop('GWAS', axis=1).to_numpy(),
gamma_values=genetic_effects['GWAS'].to_numpy(),
eqtl_participants=qtl_sample_size,
gwas_participants=gwas_sample_size,
ld_matrix=linkage,
use_pseudoInverse=False,
compute_device='cpu')
result_table = pd.DataFrame()
result_table['Genetic_Element'] = genetic_effects.columns.drop('GWAS')
result_table['Effect_Size'] = analysis_result.Alpha
result_table['Standard_Error'] = analysis_result.Se
result_table['P_Value'] = analysis_result.Pval
result_table['Heterogeneity_PValue'] = analysis_result.HetP
display(result_table)
# Genetic_Element Effect_Size Standard_Error P_Value Heterogeneity_PValue
#0 ENSG00000000419 -0.005793 0.006057 0.338818 1.0
#1 ENSG00000101126 -0.013730 0.005323 0.009903 1.0
결과 표의 헤더 설명은 공식 설명서를 참조하십시오:
Effect_Size: TWMR로 추정된 인과 효과.Standard_Error: Effect_Size의 표준 오차.P_Value: Effect_Size와 Standard_Error로 계산된 p값.D: QTL에 대한 Cohrian Q 통계량.Heterogeneity_PValue: 이질성 검정의 p값.
③ TWMR 결과를 숲 그림(forest plot)으로 시각화합니다:
# 신뢰 구간 계산
result_table['CI_lower'] = result_table['Effect_Size'] - 1.96 * result_table['Standard_Error']
result_table['CI_upper'] = result_table['Effect_Size'] + 1.96 * result_table['Standard_Error']
# 숲 그림 시각화
fig, ax = plt.subplots(figsize=(5, 2))
y_positions = [yi+1 for yi in range(len(result_table))]
plt.errorbar(result_table['Effect_Size'], y_positions,
xerr=[result_table['Effect_Size'] - result_table['CI_lower'], result_table['CI_upper'] - result_table['Effect_Size']],
fmt='o', color='blue', ecolor='black', capsize=5)
plt.yticks(y_positions, result_table['Genetic_Element'])
plt.ylim(0, max(range(len(result_table['Genetic_Element'])))+2)
plt.axvline(x=0, linestyle='--', color='blue')
plt.xlabel('인과 효과 (Effect_Size)')
plt.title('숲 그림 (Forest Plot)')
plt.show()
4 RevTWMR 분석
① RevTWMR 분석의 입력은 차이가 있지만, 여전히 GWAS와 eQTL 데이터에서 옵니다:
genetic_data = pd.read_csv("datasets/revTWMR/effect.matrix.tsv", sep='\t', index_col=0)
genetic_data.iloc[0:5, -5:]
# 표본 크기
gwas_effects = genetic_data.BETA_GWAS.values
gwas_sizes = genetic_data.N.values
sample_data = pd.read_csv("datasets/revTWMR/genes.N.tsv", sep='\t', index_col=0).N.to_dict()
sample_data
#{'ENSG00000000003': 9047.50387596899,
# 'ENSG00000000419': 30201.3203125,
#...
# 'ENSG00000070759': 30987.2519685039,
# 'ENSG00000070761': 30780.3671875,
# ...}
▲ 코드에서 genetic_data 표 설명: 각 행은 SNP를 나타내고, 각 열은 유전자를 나타내며, 채워진 값은 QTL 결과에서 SNP가 유전자 발현에 대한 효과입니다. TWMR와 다른 점은 TWMR 입력 표의 마지막 열 GWAS를 BETA_GWAS로 변경하고, 표준 오차 SE와 연구 표본 수 N을 추가해야 합니다.
또한 각 유전자에 대한 sample_data 파일을 제공해야 하며, 각 유전자에 해당하는 QTL 노출 표본 크기를 설명합니다:
- 각 유전자에 대한 표준화된 연구 크기
② RevTWMR 분석을 실행하고 결과를 하나의 표로 통합합니다:
analysis_results = {}
# 주의: 아래 예시 코드는 첫 5개 유전자만 분석합니다
for genetic_element in genetic_data.drop(['BETA_GWAS', 'SE', 'N'], axis=1).columns[:5]:
effect_table = genetic_data.drop(['BETA_GWAS', 'SE', 'N'], axis=1)[genetic_element].values
qtl_size = sample_data[genetic_element]
analysis_result = revTWMR_analysis(
effect_table,
gwas_effects,
snp_labels=genetic_data.index.values,
gwas_sample_counts=gwas_sizes,
qtl_exposure_size=qtl_size,
p_value_threshold=0.05,
use_pseudoInverse=False, compute_device='cpu')
analysis_results[genetic_element] = {
'Original_Effect_Size': analysis_result.Alpha,
'Original_Standard_Error': analysis_result.Se,
'Original_P_Value': analysis_result.Pval,
'Original_Sample_Size': analysis_result.N,
'Original_Heterogeneity_P': analysis_result.HetP,
'Iterative_Effect_Size': analysis_result.AlphaIterative,
'Iterative_Standard_Error': analysis_result.SeIterative,
'Iterative_P_Value': analysis_result.PvalIterative,
'Iterative_Heterogeneity_P': analysis_result.HetPIterative,
'Final_Sample_Size': len(analysis_result.rsname)
}
final_result = pd.DataFrame.from_dict(analysis_results, orient='index')
final_result['Genetic_Element'] = final_result.index.tolist()
final_result
# Original_Effect_Size Original_Standard_Error Original_P_Value Original_Sample_Size Original_Heterogeneity_P Iterative_Effect_Size Iterative_Standard_Error Iterative_P_Value Iterative_Heterogeneity_P Final_Sample_Size Genetic_Element
#ENSG00000000003 -0.024712 0.068517 0.718346 101.0 1.000000 -0.024712 0.068517 0.718346 1.000000 101 ENSG00000000003
#ENSG00000000419 -0.042308 0.035751 0.236646 123.0 0.957932 -0.042308 0.035751 0.236646 0.957932 123 ENSG00000000419
#ENSG00000000457 0.007607 0.035183 0.828817 123.0 0.942955 0.007607 0.035183 0.828817 0.942955 123 ENSG00000000457
결과 표의 헤더 설명은 TWMR와 거의 동일하며, 이름도 기본적으로 일치합니다.
③ RevTWMR 결과를 숲 그림(forest plot)으로 시각화합니다:
# 신뢰 구간 계산
final_result['CI_lower'] = final_result['Iterative_Effect_Size'] - 1.96 * final_result['Iterative_Standard_Error']
final_result['CI_upper'] = final_result['Iterative_Effect_Size'] + 1.96 * final_result['Iterative_Standard_Error']
# 숲 그림 시각화
fig, ax = plt.subplots(figsize=(5, 2))
y_positions = [yi+1 for yi in range(len(final_result))]
plt.errorbar(final_result['Iterative_Effect_Size'], y_positions,
xerr=[final_result['Iterative_Effect_Size'] - final_result['CI_lower'], final_result['CI_upper'] - final_result['Iterative_Effect_Size']],
fmt='o', color='blue', ecolor='black', capsize=5)
plt.yticks(y_positions, final_result['Genetic_Element'])
plt.ylim(0, max(range(len(final_result['Genetic_Element'])))+2)
plt.axvline(x=0, linestyle='--', color='blue')
plt.xlabel('인과 효과 (Effect_Size)')
plt.title('숲 그림 (Forest Plot)')
plt.show()
데이터 준비는 비교적 명확합니다
사용법도 비교적 간단합니다
오늘은 여기까지 공유하겠습니다