MetaXcan 實作教學 | GUAN-WEI Aston LI

使用流程

Harmonization：將 PLINK 輸出的 GWAS 結果轉換為 MetaXcan 所需的格式。
Imputation：對 GWAS 數據進行插補（如果必要），填補缺失的基因型數據。
Merge：整合插補後的數據，生成一個完整的 GWAS 結果文件。
SPrediXcan：使用預訓練的基因表達模型，結合 GWAS 數據，計算每個基因的 eQTL 值，並識別出顯著相關的基因。
SMulTiXcan：對多個組織的 eQTL 結果進行整合分析，找出跨組織顯著相關的基因。

Harmonization

將plink輸出的GWAS結果整理成MetaXcan後面分析需要的格式。

python ./gwas_parsing.py \
 -gwas_file /plink_output.phenotype.glm.logistic.hybrid \ #plink 輸出
 -liftover $liftover_hg19_hg38 \ #liftover 的chain file (是需求而定)
 -snp_reference_metadata ./variant_metadata.txt.gz METADATA \ # refrence panel 可從MetaXcan github下載
 -output_column_map ID variant_id \ #定義 plink 輸出的columns 名字
 -output_column_map REF non_effect_allele \ #定義 plink 輸出的columns 名字
 -output_column_map ALT1 effect_allele \ #定義 plink 輸出的columns 名字
 -output_column_map OR effect_size \ #定義 plink 輸出的columns 名字
 -output_column_map P pvalue \ #定義 plink 輸出的columns 名字
 -output_column_map CHROM chromosome \ #定義 plink 輸出的columns 名字
 --chromosome_format \
 -output_column_map POS position \ #定義 plink 輸出的columns 名字
 -output_column_map A1_FREQ frequency \ #定義 plink 輸出的columns 名字
 --insert_value sample_size 15492 --insert_value n_cases 1137 \ #sample size
 -output_order variant_id panel_variant_id chromosome position effect_allele non_effect_allele frequency pvalue zscore effect_size standard_error sample_size n_cases \ # output 的順序
 -output ./output_gwas.txt.gz # output 檔案名稱

Imputation

如果做GWAS的時候沒有imputation則可以在此步驟imputation，但是這步驟有缺點 (作者使用MetaXcan做imputation之後的結果effect size都變成NA，原因不明知道的歡迎留言)。

for i in range(1,22):
    for j in range(10)
python ./gwas_summary_imputation.py \
 -by_region_file ./eur_ld.bed.gz \ # region file 可從MetaXcan github下載
 -gwas_file ./output_gwas.txt.gz \ #上一步驟輸出
 -parquet_genotype ./chr$i.variants.parquet \ # refrence panel 可從MetaXcan github下載
 -parquet_genotype_metadata ./variant_metadata.parquet \ # refrence panel 可從MetaXcan github下載
 -window 100000 \
 -parsimony 7 \
 -chromosome $i \
 -regularization 0.1 \
 -frequency_filter 0.01 \
 -sub_batches 10 \
 -sub_batch $j \
 -standardise_dosages \
 -output ./summary_imputation/output_gwas_chr{1}_sb{2}_reg0.1_ff0.01_by_region.txt.gz

Merge

將上步驟的imputation結果整合成為一個檔案。

python ./gwas_summary_imputation_postprocess.py \
 -gwas_file ./output_gwas.txt.gz \
 -folder ./summary_imputation \
 -pattern output_gwas.* \
 -parsimony 7 \
 -output ./processed_summary_imputation/imputed_output_gwas.txt.gz

SPrediXcan

做重要的步驟，利用Database計算eQTL。利用bonferroni correction找出顯著的基因，或是利用P 和 effect size 做出火山圖。

parallel --jobs 56 \
    python ./SPrediXcan.py \
    --gwas_file ./processed_summary_imputation/imputed_output_gwas.txt.gz \
    --snp_column panel_variant_id --effect_allele_column effect_allele --non_effect_allele_column non_effect_allele --zscore_column zscore \
    --model_db_path ./models/eqtl/mashr/mashr_{1}.db \
    --covariance ./models/eqtl/mashr/mashr_{1}.txt.gz \
    --keep_non_rsid --additional_output --model_db_snp_key varID \
    --throw \
    --output_file ./spredixcan/eqtl/imputed_output_gwas__PM__{1}.csv \
    ::: $(cat ${target})

SMulTiXcan

將上一步驟產生的多組織eQTL統整起來並利用bonferroni correction找出顯著的基因。

python ./SMulTiXcan.py \
--models_folder ./models/eqtl/mashr \
--models_name_pattern "mashr_(.*).db" \
--snp_covariance ./models/gtex_v8_expression_mashr_snp_smultixcan_covariance.txt.gz \
--metaxcan_folder ./spredixcan_heart/eqtl/ \
--metaxcan_filter "imputed_CAD_TWB_gwas__PM__(.*).csv" \
--metaxcan_file_name_parse_pattern "(.*)__PM__(.*).csv" \
--gwas_file ./processed_summary_imputation/imputed_output_gwas.txt.gz \
--snp_column panel_variant_id --effect_allele_column effect_allele --non_effect_allele_column non_effect_allele --zscore_column zscore --keep_non_rsid --model_db_snp_key varID \
--cutoff_condition_number 30 \
--verbosity 7 \
--throw \
--output ./smultixcan/eqtl/imputed_output_gwas_smultixcan.txt

Ref: https://github.com/hakyimlab/MetaXcan