PrediXcan , SPrediXcan，MetaXcan是近些年基於GWAS后續分析開發出來的工具。
主要功能是在組織和表達的層面預測影響表型的基因，彌補了GWAS只能在基因組層面解釋表型的不足。
下面是這幾個工具的工作流程：
今天給大家介紹一下如何使用SPrediXcan和MetaXcan工具進行全轉錄組分析（Transcriptome-Wide Analysis）
該工具最大的優點是不需要個體水平的基因型數據和表型數據，只需要提供GWAS的summary統計文件就可以預測影響表型的靶基因。
缺點是對電腦的配置要求比較高。
1、配置MetaXcan環境
下載、解壓MetaXcan：
wget https://github.com/hakyimlab/MetaXcan/archive/master.zip
unzip master.zip
cd MetaXcan-master/software/
使用conda配置MetaXcan的環境（沒安裝conda的請自行安裝后再做這一步）：
conda env create -n MetaXcan -f conda_env.yaml
激活MetaXcan
conda activate MetaXcan
2、配置gwasimputation環境
下載、解壓gwasimputation：
wget https://github.com/hakyimlab/summary-gwas-imputation/archive/master.zip
unzip master.zip
配置環境：
conda env create -n gwasimputation -f conda_env.yaml
conda activate gwasimputation
conda install -c conda-forge pyliftover=0.4
conda install -c conda-forge pandas=0.25.3
conda install -c conda-forge scipy=1.4.1
conda install -c conda-forge numpy=1.18.1
conda install -c conda-forge bgen_reader=3.0.2
conda install -c conda-forge cyvcf2=0.20.0
conda install -c conda-forge pyarrow=0.9.0
各位仔細觀察的話，會發現我指定了安裝包的版本。
建議各位嚴格安裝這個教程的版本安裝，不然后面會出現各種因版本過高或過低的報錯。
3、下載、解壓測試數據
nohup wget https://zenodo.org/record/3657902/files/sample_data.tar?download=1 &
tar -xvpf sample_data.tar\?download\=1
解壓以后，會出現如下的數據：
gwas文件夾包含GWAS結果文件的測試數據cad.add.160614.website.txt.gz；
完成以上工作以后，我們就可以開始分析數據啦！
以下是全轉錄組分析的流程
4、協調GWAS文件的坐標（hg19坐標統一為hg38）
原始GWAS結果文件cad.add.160614.website.txt.gz的格式如下所示：
我們可以看到，這個文件的坐標是hg19版本的，但MetaXcan的最新版是hg38，因此我們需要先將GWAS結果文件cad.add.160614.website.txt.gz的坐標轉為hg38。
這里面，我們需要先定義幾個文件夾：
GWAS_TOOLS=/chenwenyan/software/gwas-imputation-master/src #步驟二的gwas-imputation-master的路徑
DATA=/chenwenyan/anaconda2/envs/gwasimputation/data/data #步驟三下載的測試數據的路徑
OUTPUT=/chenwenyan/test #這個自己定義：輸出文件的路徑
注意，以上路徑用的是我自己的路徑（比如
/chenwenyan/software/gwas-imputation-master/src），不要照搬。請修改成你自己的路徑。
定義完路徑以后，在輸出文件夾$OUTPUT（我的示例是/chenwenyan/test）下新建一個文件夾創建文件夾harmonized_gwas，輸入命令：mkdir harmonized_gwas
隨后，創建一個test.sh的文件，輸入命令：
vi test.sh
在test.sh文件輸入以下內容：
#!/bin/bash
GWAS_TOOLS=/chenwenyan/software/gwas-imputation-master/src #gwas-imputation-master的路徑
DATA=/chenwenyan/anaconda2/envs/gwasimputation/data/data #下載的測試數據的路徑
OUTPUT=/chenwenyan/test #輸出文件的路徑
python $GWAS_TOOLS/gwas_parsing.py \
-gwas_file $DATA/gwas/cad.add.160614.website.txt.gz \
-liftover $DATA/liftover/hg19ToHg38.over.chain.gz \
-snp_reference_metadata $DATA/reference_panel_1000G/variant_metadata.txt.gz METADATA \
-output_column_map markername variant_id \
-output_column_map noneffect_allele non_effect_allele \
-output_column_map effect_allele effect_allele \
-output_column_map beta effect_size \
-output_column_map p_dgc pvalue \
-output_column_map chr chromosome \
--chromosome_format \
-output_column_map bp_hg19 position \
-output_column_map effect_allele_freq frequency \
--insert_value sample_size 184305 --insert_value n_cases 60801 \
-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/harmonized_gwas/CARDIoGRAM_C4D_CAD_ADDITIVE.txt.gz
GWAS結果文件的染色體用數值表示，不要用
chr+數值的方式GWAS文件格式必須是gz壓縮格式，分隔符為tab
運行：sh test.sh
順利的話，會輸入以下信息：
生成的新文件如下所示：
如截圖的紅框所示，此時坐標已變為hg38了。
5、對gwas結果文件進行imputation
先在輸出文件夾$OUTPUT（我的示例是/chenwenyan/test）下新建一個文件夾summary_imputation,輸入命令：mkdir summary_imputation
隨后輸入如下命令進行imputation：
GWAS_TOOLS=/chenwenyan/software/gwas-imputation-master/src
DATA=/chenwenyan/anaconda2/envs/gwasimputation/data/data
OUTPUT=/chenwenyan/test
python $GWAS_TOOLS/gwas_summary_imputation.py \
-by_region_file $DATA/eur_ld.bed.gz \
-gwas_file $OUTPUT/harmonized_gwas/CARDIoGRAM_C4D_CAD_ADDITIVE.txt.gz \
-parquet_genotype $DATA/reference_panel_1000G/chr1.variants.parquet \
-parquet_genotype_metadata $DATA/reference_panel_1000G/variant_metadata.parquet \
-window 100000 \
-parsimony 7 \
-chromosome 1 \
-regularization 0.1 \
-frequency_filter 0.01 \
-sub_batches 10 \
-sub_batch 0 \
--standardise_dosages \
-output $OUTPUT/summary_imputation/CARDIoGRAM_C4D_CAD_ADDITIVE_chr1_sb0_reg0.1_ff0.01_by_region.txt.gz
出現如下頁面說明順利運行了：
上面示例只是對一號染色體上的第一個批次進行imputation，我們實際上要進行220次（22條染色體*10個批次）循環imputation。
因此跑完上面流程以后，需要繼續跑剩下的119次imputation。
需要修改的地方有
$DATA/reference_panel_1000G/chr1.variants.parquet、-chromosome 1、-sub_batch 0、$OUTPUT/summary_imputation/CARDIoGRAM_C4D_CAD_ADDITIVE_chr1_sb0_reg0.1_ff0.01_by_region.txt.gz
$DATA/reference_panel_1000G/chr1.variants.parquet修改的數值是chr1，chr2，chr3，……，chr22；
-chromosome 1修改的數值是：1,2,3，……，22；
-sub_batch 0修改的數值是：0,1,2,……，9；
$OUTPUT/summary_imputation/CARDIoGRAM_C4D_CAD_ADDITIVE_chr1_sb0_reg0.1_ff0.01_by_region.txt.gz是輸出文件，建議修改chr1_sb0為對應的染色體和批次。
舉個例子，我現在想跑3號染色體的第十個批次，那么腳本應該如下所示：
GWAS_TOOLS=/chenwenyan/software/gwas-imputation-master/src
DATA=/chenwenyan/anaconda2/envs/gwasimputation/data/data
OUTPUT=/chenwenyan/test
python $GWAS_TOOLS/gwas_summary_imputation.py \
-by_region_file $DATA/eur_ld.bed.gz \
-gwas_file $OUTPUT/harmonized_gwas/CARDIoGRAM_C4D_CAD_ADDITIVE.txt.gz \
-parquet_genotype $DATA/reference_panel_1000G/chr3.variants.parquet \
-parquet_genotype_metadata $DATA/reference_panel_1000G/variant_metadata.parquet \
-window 100000 \
-parsimony 7 \
-chromosome 3 \
-regularization 0.1 \
-frequency_filter 0.01 \
-sub_batches 10 \
-sub_batch 9 \
--standardise_dosages \
-output $OUTPUT/summary_imputation/CARDIoGRAM_C4D_CAD_ADDITIVE_chr3_sb9_reg0.1_ff0.01_by_region.txt.gz
6、合並220個批次的數據
輸入如下命令：
GWAS_TOOLS=/chenwenyan/software/gwas-imputation-master/src
DATA=/chenwenyan/anaconda2/envs/gwasimputation/data/data
OUTPUT=/chenwenyan/test
python $GWAS_TOOLS/gwas_summary_imputation_postprocess.py \
-gwas_file $OUTPUT/harmonized_gwas/CARDIoGRAM_C4D_CAD_ADDITIVE.txt.gz \
-folder $OUTPUT/summary_imputation \
-pattern CARDIoGRAM_C4D_CAD_ADDITIVE.* \
-parsimony 7 \
-output $OUTPUT/processed_summary_imputation/imputed_CARDIoGRAM_C4D_CAD_ADDITIVE.txt.gz
和上面一樣，路徑不要照搬我的，請修改成你自己的路徑。
正常運行的話，會出現如下界面：
7、使用單個組織預測基因表達
在輸出文件夾$OUTPUT（我的示例是/chenwenyan/test）下新建一個文件夾spredixcan,輸入命令：mkdir spredixcan
在spredixcan文件夾下新建文件夾eqtl，輸入命令：mkdir eqtl
建完后，輸入如下命令：
DATA=/chenwenyan/anaconda2/envs/gwasimputation/data/data
METAXCAN=/chenwenyan/software/MetaXcan-master/software
python $METAXCAN/SPrediXcan.py \
--gwas_file $OUTPUT/processed_summary_imputation/imputed_CARDIoGRAM_C4D_CAD_ADDITIVE.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 $DATA/models/eqtl/mashr/mashr_Whole_Blood.db \
--covariance $DATA/models/eqtl/mashr/mashr_Whole_Blood.txt.gz \
--keep_non_rsid --additional_output --model_db_snp_key varID \
--throw \
--output_file $OUTPUT/spredixcan/eqtl/CARDIoGRAM_C4D_CAD_ADDITIVE__PM__Whole_Blood.csv
正常運行的話，會出現如下界面：
8、使用多個組織預測基因表達
在輸出文件夾$OUTPUT（我的示例是/chenwenyan/test）下新建一個文件夾smultixcan，輸入命令：mkdir smultixcan
在smultixcan文件夾下新建文件夾eqtl，輸入命令：mkdir eqtl
python $METAXCAN/SMulTiXcan.py \
--models_folder $DATA/models/eqtl/mashr \
--models_name_pattern "mashr_(.*).db" \
--snp_covariance $DATA/models/gtex_v8_expression_mashr_snp_covariance.txt.gz \
--metaxcan_folder $OUTPUT/spredixcan/eqtl/ \
--metaxcan_filter "CARDIoGRAM_C4D_CAD_ADDITIVE__PM__(.*).csv" \
--metaxcan_file_name_parse_pattern "(.*)__PM__(.*).csv" \
--gwas_file $OUTPUT/processed_summary_imputation/imputed_CARDIoGRAM_C4D_CAD_ADDITIVE.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 $OUTPUT/smultixcan/eqtl/CARDIoGRAM_C4D_CAD_ADDITIVE_smultixcan.txt
9、結果演示
生成的結果如下所示：
結果文件標題所代表的意思：
gene: Ensembl id;
gene_name: HUGO gene name;
Z-score: S-PrediXcan's association result for the gene;
pvalue: P-value of the aforementioned statistic;
var_g: variance of the gene expression, calculated as W' * G * W (where W is the vector of SNP weights in a gene's model, W' is its transpose, and G is the covariance matrix);
pred_perf_r2: R2 of tissue model's correlation to gene's measured transcriptome (prediction performance); pred_perf_pval: pval of tissue model's correlation to gene's measured transcriptome (prediction performance);
n_snps_used: number of snps from GWAS that got used in S-PrediXcan analysis;
n_snps_in_cov: number of snps in the covariance matrix;
n_snps_in_model: number of snps in the model.
結果文件的pvalue值超過Bonferroni correction的閾值就說明該結果有意義。
Bonferroni correction定義：0.05/（N個基因*M個組織）
