Project Name: Frozen-GWAS (named by Dr. Canzian): Multi-ancestry GWAS meta-analysis to identify effector genes for pancreatic cancer
Description: This repository contains the scripts and pipeline used in the multi-ancestry pancreatic cancer GWAS manuscript. Please follow the scripts listed below to generate the analysis results described in the manuscript.
Citation:
"Güler, M., et al. A meta-analysis of genome-wide association studies identify fourteen new pancreatic cancer risk loci. XX (2026). https://doi.org/XX"
| Software | Version | Link/Citation |
|---|---|---|
| regenie | 3.2.1 | regenie |
| plink | 1.90b6.21 | plink |
| plink2 | 2.00a3LM | plink2 |
| R | 4.3.0 or above | R |
| METAL | 2020-05-05 | METAL |
| MR-MEGA | v0.2 | MR-MEGA |
| GCTA | 1.94.1 | GCTA |
| LDSC | v1.0.1 | LDSC |
| GWASLab | v3.4.31 | GWASLab |
| FLAMES | v1.1.2 | FLAMES |
| Open4Gene | 0.1.1 | Open4Gene |
| R Package | Version | Link/Citation |
|---|---|---|
| tidyverse | 2.0.0 | tidyverse |
| ggplot2 | 3.1.3.1 | gplots |
| data.table | 1.14.10 | data.table |
| parallel | 4.3.0 | parallel |
| magrittr | 2.0.3 | magrittr |
| colorspace | 2.1.0 | colorspace |
| susieR | 0.12.35 | susieR |
Note
In this study, the scripts are created for IBM LSF job scheduler. These scripts can easily be converted to commonly used schedulers such as SLURM or PBS. Without having access to the multi-node cluster envrioment, it would not be possible to run some analysis. Thank you to the DKFZ Cluster, operated by the Omics IT and Data Management Core Facility (ODCF), for providing computational infrastructure.
To reproduce the results and figures presented in the manuscript, nearly all scripts (some runner script or simple data handling scripts may not included, please let us know if you need help) used are provided in this repository.
| Step | Short Summary |
|---|---|
| 01 | Imputation and GWAS (PanScan I/II, PanC4) |
| 02 | Format and harmonize GWAS summary statistics |
| 03 | Fixed-effect meta-analysis (EUR and multi-ancestry) |
| 04 | Trans-ancestry meta-analysis with MR-MEGA |
| 05 | Diagnostic GWAS plots and preparation of FUMA inputs |
| 06 | GCTA-cojo |
| 07 | Finemapping-susieR |
| 08 | Finemapping-MR-MGA |
| 09 | LDSC |
| 10 | FLAMES |
| 11 | Open4Gene |
| 12 | caABC |
| 13 | TWAS |
| 14 | PRS |
| 15 | Plots |
| 16 | Others |
In this step you can find scripts used for prepare PANSCAN I/II and PC4 for TOPMed r3 imputation, download results from imputation server and QC them, and finaly run GWAS with regenie.
Key operations include:
- Creating vcf files for TOPMed imputation
- Processing imputation results and QC
- GWAS with regenie
Note
PanScan I/II and PanC4 data imputation strategy similar to the PMID: 34508597. Briefly, we follow two-step imputation. In the first step, we imputed seperately the three cohort and merged them using < 1% missingness, infoscore > 99, hwe 1e-10. In the second step, we merged three data set and re-imputated merged the merged dataset.
List of script and their functions for the Step 1:
Example commands:
bsub < 01_Imputation-GWAS-PanScan-PC4/01_plinkmakevcf.bsub -R "rusage[mem=10GB]"
bsub < 01_Imputation-GWAS-PanScan-PC4/02_downloadtopmed.bsub -R "rusage[mem=10GB]"
bsub < 01_Imputation-GWAS-PanScan-PC4/03_plinkmakebed.bsub -R "rusage[mem=64GB]"
# to fix psam files
for chr in {1..22}; do
awk 'BEGIN {OFS="\t"}
NR==1 {print "#FID","IID","SEX"; next}
NR>1 {print 0, $1, $2}' \
mg_PS12PC4.chr${chr}.psam > tmp.psam && mv tmp.psam mg_PS12PC4.chr${chr}.psam
done
bsub < 01_Imputation-GWAS-PanScan-PC4/04_regenie_step1.bsub -R "rusage[mem=10GB]"
bsub < 01_Imputation-GWAS-PanScan-PC4/05_regenie_step2_PS1.bsub -R "rusage[mem=10GB]"
bsub < 01_Imputation-GWAS-PanScan-PC4/06_regenie_step2_PS2.bsub -R "rusage[mem=10GB]"
bsub < 01_Imputation-GWAS-PanScan-PC4/07_regenie_step2_PC4.bsub -R "rusage[mem=10GB]"
# to combine regenie results
for prefix in mg_PC4 mg_PS1 mg_PS2; do
# Extract header from the first chromosome
zcat ${prefix}_chr1_PHENO.regenie.gz | head -n 1 | sed 's/ /\t/g' > ${prefix}.tab.regenie
# Append data from all chromosomes
for chr in {1..22}; do
zcat ${prefix}_chr${chr}_PHENO.regenie.gz | sed 1d | sed 's/ /\t/g' >> ${prefix}.tab.regenie
done
done
This step standardizes and harmonizes cohort-level GWAS summary statistics to enable meta-analysis across studies and ancestries.
Key operations include:
- Converting cohort GWAS outputs to a common summary statistics schema (required fields, consistent naming, and data types)
- Harmonizing alleles and effect directions
- Liftover between genome builds (GRCh37/hg19 and GRCh38/hg38), as needed
- JaPAN: aligning effect allele orientation and adding East Asian allele frequencies (AF) from gnomAD v3.1.2
- CKB: formatting and moving final outputs into the GRCh38 clean summary statistics folder
List of scripts and their functions for Step 2:
Example commands:
Rscript --slave --no-restore --no-save 02_format_harmonize_sumstats/08_PS1PS2PC4_sumstat_format.R
cd ./trans_ancestry/hg19_clean_sumstats
wget https://hgdownload.cse.ucsc.edu/goldenpath/hg19/bigZips/hg19.fa.gz
gunzip hg19.fa.gz
bsub < 02_format_harmonize_sumstats/09_harmonize_JaPAN_CKB.bsub -R "rusage[mem=32GB]"
# Download gnomAD 3.1.2 vcf files
bsub < 02_format_harmonize_sumstats/10_gnomad_getVCF.bsub -R "rusage[mem=1GB]"
# extract east asian AF from VCFs:
bsub < 02_format_harmonize_sumstats/11_extract_eas_AF.bsub -R "rusage[mem=8GB]"
#After script 11, each chr EAS_AF need to be filtered (\$5) and merged to one-file.
# The fastest way to do that is using awk!
#!/bin/bash
# Set input/output directories
EAS_DIR="/your/path/gnomAD/eas_af"
OUT_FILE="eas_all_merged.tsv"
# Initialize output with header
echo -e "CHR\tPOS\tREF\tALT\tAF_EAS\tSNP" > "$OUT_FILE"
# Loop over each chromosome file
for FILE in "$EAS_DIR"/eas_af_chr*.tsv; do
echo "Processing $FILE..."
# Process each line, skipping bad AF_EAS values and adding SNP column
awk -F'\t' 'BEGIN {OFS="\t"}
$5 != "." && $5 != "" && $5 > 0 {
gsub(/^chr/, "", $1);
snp = $1 ":" $2 ":" $3 ":" $4;
print $1, $2, $3, $4, $5, snp
}' "$FILE" >> "$OUT_FILE"
done
echo "All files merged into $OUT_FILE"
# In this script also format CKB data to move final hg38 folder.
Rscript --slave --no-restore --no-save 12_CKB_format_combine_JaPAN_gnomAD.RNote
All raw cohort summary statistics are stored in trans_ancestry/raw_sumstats/.
In this step, GWAS summary statistics from individual cohorts were combined using fixed-effects meta-analysis implemented in METAL.
Meta-analyses were performed separately for:
- European ancestry cohorts
- Trans-ancestry cohorts
The METAL parameter files define study-specific inputs, allele alignment rules, and sample size weighting.
Filtering and post-processing applied to METAL outputs included:
- Minimum number of contributing studies per variant (HetDf ≥ 1; at least two cohorts)
- Minor allele frequency threshold (Freq1 ≥ 0.01 unless otherwise specified)
- Removal of variants with inconsistent allele orientation (fixing allele flips)
- Optional filtering based on heterogeneity P-value (HetP > 0.05) for sensitivity analyses
List of scripts and their functions for Step 3:
Example commands:
bsub < 03_metal/15_metal.bsub -R "rusage[mem=8GB]"
# Hetdf \>= 1 (at least two studies), Freq1 \>= 0.01, HetPVal >= 0.05 (also fixing allele flips REF=NEA, ALT=EA)
bsub < 03_metal/16_format_metal_results.bsub -R "rusage[mem=35GB]"Note
METAL has internal behavior that may flip alleles depending on strand orientation; therefore, additional checks were applied to ensure the reported effect allele corresponds to the ALT allele encoded in the SNP identifier. For METAL results, allele orientation was cross-checked using the SNP identifier (CHR:POS:REF:ALT) to ensure the reported effect allele corresponds to ALT. When needed, alleles and effect directions were corrected consistently across all downstream steps.
To account for ancestry-associated heterogeneity in effect sizes, we performed trans-ancestry meta-regression using MR-MEGA.
Key steps included:
- Formatting cohort-level GWAS summary statistics into MR-MEGA–compatible inputs
- Determining the optimal number of principal components (PCs) describing ancestry structure
- Running MR-MEGA using the selected PCs
- Recalculating extremely small p-values beyond the numerical limits of the MR-MEGA C++ implementation
The number of PCs was determined empirically. An initial diagnostic run using the maximum allowable number of PCs (number of cohorts − 3) was performed, and study-specific PCs were visualized. Based on separation of cohorts and populations in PC space, three PCs were selected for the final analysis.
Note
MR-MEGA reports association p-values with a lower bound of ~1×10⁻¹⁴; p-values smaller than this threshold were recalculated in R using chi-square statistics and the corresponding degrees of freedom.
List of scripts and their functions for Step 4:
Example commands:
# Download MR-MEGA
readonly MR_MEGA_URL="https://tools.gi.ut.ee/tools/MR-MEGA_v0.2.zip"
readonly MR_MEGA_ZIP="$(basename ${MR_MEGA_URL})"
readonly MR_MEGA_FILE="${MR_MEGA_ZIP%.zip}"
mkdir -p ~/bin/"${MR_MEGA_FILE}"
cd ~/bin/"${MR_MEGA_FILE}"
echo ${MR_MEGA_FILE}
curl -O "${MR_MEGA_URL}"
echo $(ls)
echo ${MR_MEGA_ZIP}
unzip "${MR_MEGA_ZIP}"
make
cd ..
ln -s "${MR_MEGA_FILE}"/MR-MEGA MR-MEGA
echo "MR-MEGA is now installed."
fi
# Format sumstats for MR-MEGA
bsub < 04_mr_mega/17_mr_mega_input.bsub -R "rusage[mem=32GB]"
# crete .in file
ls $PWD/*hg38_MRMEGA.txt | grep -v "MVP_EUR_hg38_MRMEGA.txt" | sort > MR-MEGA.in
# Determine number of PC
bsub < 04_mr_mega/18_mr_mega_maxPC.bsub -R "rusage[mem=64GB]"
# After running MR-MEGA with maxPC we can obtain PCs for each studies in the log file.
grep -n 'omics' \*.MAX_PCs.log
# To plot PCs:
Rscript --slave --no-restore --no-save 04_mr_mega/18_plot_mrmega_PCs.R
# numPCs parameter in the MR_MEGA selected as *3*
# run MR-MEGA with PC1-PC3
bsub < 04_mr_mega/19_mr_mega.bsub -R "rusage[mem=64GB]"
# Fix P values
wget https://tools.gi.ut.ee/tools/fixP.r
mv fixP.r 20_mr_mega_fixP.R
# The original script modified for speed(eg- data.table for r/w.)
# Also new script removes P-value_association = NA
Rscript --slave --no-restore --no-save 04_mrmega/20_mr_mega_fixP.R
# Formatting results
Rscript --slave --no-restore --no-save 04_mrmega/21_format_mrmega_results.RThis step generates a comprehensive set of diagnostic plots to compare results across different GWAS analyses and prepares summary statistics for downstream tools that require specific genome builds.
Using GWASLab, the following visualizations were created:
- Stratified Manhattan and QQ plots for:
- European fixed-effects meta-analysis
- Trans-ancestry fixed-effects meta-analysis
- MR-MEGA meta-analysis
- Stacked Manhattan plots combining all analyses
- Miami plots comparing:
- European vs trans-ancestry fixed-effects results
- European fixed-effects vs MR-MEGA results
- Trans-ancestry fixed-effects vs MR-MEGA results
- Scatter plots comparing effect sizes between analyses
In addition, summary statistics were lifted over from GRCh38 (hg38) to GRCh37 (hg19) to enable compatibility with downstream tools such as FUMA, which require hg19 input.
Note
Liftover was performed using chain files and validated by checking chromosome, position, and allele consistency before downstream analyses.
List of scripts and their functions for Step 5:
| Script | Function |
|---|---|
| 22 | Generate diagnostic plots and liftover for downstream tools |
Example command:
bsub < 05_gwas_plots/22_diagplots_fuma.bsub -R "rusage[mem=64GB]"To identify conditionally independent association signals within genome-wide significant loci, we performed conditional and joint analysis (COJO) using GCTA.
We used the PS1–PS2–PC4 reference genotype dataset as the LD reference panel. Although this reference panel consists of individuals of European ancestry, it was applied consistently to both:
- European fixed-effects meta-analysis
- Trans-ancestry fixed-effects meta-analysis
- Trans-ancestry MR-MEGA analysis
For trans-ancestry results, the same reference was used to ensure interpretability and to provide a consistent set of conditional signals for downstream fine-mapping and gene prioritization.
Note
COJO requires an LD reference panel. Using a single reference panel across analyses ensures consistent interpretation of independent signals, but may reduce sensitivity for ancestry-specific LD patterns in multi-ancestry analyses.
List of scripts and their functions for Step 7:
| Script | Function |
|---|---|
| 23 | Prepare COJO inputs (download/prepare GCTA + .ma files) |
| 23 | Run GCTA-COJO per chromosome |
Example commands:
# download gcta tool and generate cojo input files
bsub < 06_cojo/24_cojo_input.bsub -R "rusage[mem=16GB]"
# run cojo
bsub < 06_cojo/25_cojo.bsub -R "rusage[mem=64GB]"
# Pool cojo results
head -n 1 EUR_PC_hg38_meta_results_chr1.jma.cojo > EUR_METAL_combined_1e6.jma.cojo
awk 'FNR>1' EUR_PC_hg38_meta_results_chr*.jma.cojo >> EUR_METAL_combined_1e6.jma.cojo
head -n 1 trans_ancestry_PC_hg38_meta_results_chr1.jma.cojo > TRANS_METAL_combined_1e6.jma.cojo
awk 'FNR>1' trans_ancestry_PC_hg38_meta_results_chr*.jma.cojo >> TRANS_METAL_combined_1e6.jma.cojo
head -n 1 trans_ancestry_PC_hg38_mrmega.result_chr1.jma.cojo > TRANS_MRMEGA_combined_1e6.jma.cojo
awk 'FNR>1' trans_ancestry_PC_hg38_mrmega.result_chr*.jma.cojo >> TRANS_MRMEGA_combined_1e6.jma.cojoTo identify causal variants at associated loci, we performed Bayesian fine-mapping using SuSiE (Sum of Single Effects) with summary statistics and external LD information.
Fine-mapping was performed separately for:
- European-only fixed-effects meta-analysis (using UKB LD reference)
- Trans-ancestry fixed-effects meta-analysis (using UKB LD reference)
- Prepare LD files: Convert UKB LD reference panel from
.npz(NumPy) format to R-compatible.featherformat - Run SuSiE: Execute fine-mapping to identify credible sets of variants
- Merge results: Combine per-locus fine-mapping outputs into summary tables
The LD loading pipeline was adapted from PolyFun to enable efficient processing of large-scale LD matrices in R.
Note
SuSiE requires: (1) summary statistics from GWAS, (2) linkage disequilibrium (LD) correlations, (3) locus-specific boundaries. We used a gene-centric approach to define locus boundaries for variant fine-mapping.
Note
The npz to feather format conversion adapted from polyfun. Thank you to Omer Weissbrod!
List of scripts and their functions for Step 7:
Example commands:
# Step 1: Prepare LD files from UKB reference panel
bsub < 07_susieR/25_ld_file2feather.bsub -R "rusage[mem=16GB]"
# Step 2: Run SuSiE fine-mapping
# EUR fixed-effects
bsub < 07_susieR/27_FE_EUR_susieR.bsub -R "rusage[mem=16GB]"
# Trans-ancestry fixed-effects
bsub < 07_susieR/27_FE_TA_susieR.bsub -R "rusage[mem=16GB]"
# Step 3: Merge and summarize SuSiE results
Rscript --slave --no-restore --no-save 07_susieR/28_merge_susie_res.RTo leverage the trans-ancestry information for variant fine-mapping, we extracted Bayes Factor (BF) statistics from MR-MEGA results and computed credible sets using Bayesian model selection.
List of scripts and their functions for Step 8:
| Script | Function |
|---|---|
| 29 | Fine-mapping pipeline for MR-MEGA results |
Example command:
Rscript --slave --no-restore --no-save 08_finemap_mr_mega/29_mrmega_finemap.RTo quantify the contribution of common genetic variants to pancreatic cancer risk, we performed LD Score Regression (LDSC) to estimate SNP-based heritability (h²) on both the observed and liability scales.
Heritability estimation:
- Computed per-SNP effective sample size (Neff) using: Neff = 4 × Ncase × Ncontrol / (Ncase + Ncontrol)
- Ran LDSC with European reference LD scores, using sample prevalence (--samp-prev 0.5) and population prevalence (--pop-prev 0.02)
Note
COJO, heritability (h²), and genetic correlations (rg) were estimated only for the EUR-only meta-analysis. Subsequent causal inference analyses (e.g., fine-mapping genes, SMR) used the EUR fixed-effects GWAS results.
List of scripts and their functions for Step 9:
| Script | Function |
|---|---|
| 30 | Prepare GWAS data for LDSC analysis (compute Neff, format) |
| 31 | Estimate SNP-based heritability (h²) |
Example commands:
# Prepare GWAS data for LDSC
bsub < 09_ldsc/25_gwaslab_ldsc.bsub -R "rusage[mem=16GB]"
# Estimate heritability
bsub < 09_ldsc/26_ldsc_h2.bsub -R "rusage[mem=8GB]"
Key Findings:
- Observed-scale h²: 0.0827 (95% CI: 0.057 – 0.108)
- Liability-scale h²: 0.0512 (95% CI: 0.035 – 0.067)
To identify genes with significant associations to pancreatic cancer, we performed gene-based analysis using MAGMA through FUMA web interface.
Significance threshold: P = 0.05 / 19,260 (Bonferroni-corrected) = 2.596 × 10⁻⁶
Note
MAGMA gene-based association testing aggregates SNP-level association signals within gene boundaries, reducing the multiple testing burden compared to SNP-level analysis.
FLAMES (Fine-mapping Locus Analysis via Multi-criteria Evaluation Scores) integrates multiple fine-mapping and prioritization resources to identify genes most likely to be causal at each locus.
For each locus with fine-mapped variants, we created credible set files at different PIP (Posterior Inclusion Probability) thresholds:
- Primary threshold: 95% credible set (cumulative PIP ≥ 0.95)
- Alternative threshold: 80% credible set for loci with >100 variants in the 95% set (to reduce computational burden)
- Fallback: Lead SNP only, if no credible set could be defined
FLAMES analysis was performed using:
- MR-MEGA-based fine-mapped SNP sets
- SuSiE-based fine-mapped SNP sets (EUR and trans-ancestry separately)
List of scripts and their functions for Step 10:
Example commands:
Rscript --slave --no-restore --no-save 10_flames/32_mrmega_flames_locusfinemap_files.R
bsub < 10_flames/33_mrmega_flames.bsub -R "rusage[mem=32GB]"
Rscript --slave --no-restore --no-save 10_flames/34_fe_eur_flames_locusfinemap_file.R
bsub < 10_flames/35_fe_eur_flames.bsub -R "rusage[mem=32GB]"
Rscript --slave --no-restore --no-save 10_flames/36_fe_ta_flames_locusfinemap_file.R
bsub < 10_flames/37_fe_ta_flames.bsub -R "rusage[mem=32GB]"To link fine-mapped GWAS variants with their potential target genes in pancreatic cells, we performed Open4Gene analysis on single-cell multiome (scRNA-seq + scATAC-seq) data from pancreatic tumor and normal tissue (GSE241895).
Open4Gene integrates matched gene expression (RNA) and chromatin accessibility (ATAC) at single-cell resolution to identify statistically significant peak-to-gene regulatory links across cell types using a zero-inflated negative binomial hurdle model.
- Data preprocessing: Load and QC single-cell multiome data (tumor and normal separately)
- Peak annotation and cell-type identification: Identify cell types and annotate ATAC peaks to genes
- Peak-gene pair chunking: Divide peak-gene pairs into manageable chunks for parallelization
- Open4Gene analysis: Test for peak-gene associations in each cell type using hurdle models
- Combine results: Merge results from all chunks and tissues (tumor + normal)
- Link to GWAS: Overlap fine-mapped SNPs with significant peak-gene links to prioritize candidate effector genes
Note
Open4Gene uses a two-part hurdle model to account for zero-inflation in single-cell data. Statistical significance was determined using FDR < 0.05 in either the zero-inflation or count component.
List of scripts and their functions for Step 11:
Example commands:
# Install Open4Gene and dependencies (script 38)
Rscript --slave --no-restore --no-save 11_open4gene/38_install.R
# Preprocessing and QC (tumor and normal separately) (scripts 39-40)
# These scripts load raw data, perform QC, cell-type annotation, peak annotation
# and create peak-gene pairs, then save processed objects
Rscript --slave --no-restore --no-save 11_open4gene/39_open4gene_tumor_part1.R
Rscript --slave --no-restore --no-save 11_open4gene/40_open4gene_normal_part1.R
# Run Open4Gene analysis in parallel (scripts 42-43)
# Array jobs run Open4Gene on chunked peak-gene pairs (100 chunks per tissue)
bsub < 11_open4gene/42_tumor.bsub -R "rusage[mem=32GB]"
bsub < 11_open4gene/43_normal.bsub -R "rusage[mem=32GB]"
# Combine results from all chunks (script 44)
Rscript --slave --no-restore --no-save 11_open4gene/44_combine_results.R
# Link fine-mapped SNPs to significant peak-gene pairs (script 45)
Rscript --slave --no-restore --no-save 11_open4gene/45_summary.RTo complement fine-mapping and identify candidate effector genes at PDAC risk loci, we integrated chromatin accessibility and Hi-C-based enhancer-promoter contact information using the Activity-By-Contact (ABC) model.
The ABC model predicts which genes are regulated by enhancers based on:
- Activity: Chromatin accessibility (ATAC-seq or DNase-seq) and histone marks (H3K27ac) at candidate regulatory elements
- Contact: 3D chromatin interaction frequency between enhancers and promoters (from Hi-C)
We used pancreatic-specific ABC predictions from published data (doi: 10.1016/j.pan.2025.06.001) that includes:
- Enhancer regions linked to target genes
- Promoter regions for annotated genes
- Gene bodies for comprehensive coverage
Fine-mapped SNPs from credible sets were overlapped with these regulatory elements to prioritize genes based on:
- Number of SNPs overlapping ABC enhancers linked to each gene
- Number of SNPs overlapping ABC promoters
- Number of SNPs within gene bodies
- ABC score: Composite score reflecting evidence strength (0-3)
Note
The ABC model integrates epigenomic and 3D genomic data to predict enhancer-gene connections. SNPs overlapping ABC-predicted enhancers are more likely to affect gene expression through disrupting regulatory elements.
List of scripts and their functions for Step 12:
Example commands:
# Run ABC overlap analysis (script 46)
Rscript --slave --no-restore --no-save 12_caABC/46_getABC_snp_set.RABC Score Interpretation:
- Score = 3: SNP(s) overlap ABC enhancer(s) + promoter(s) + gene body (strongest evidence)
- Score = 2: SNP(s) overlap two of the three feature types
- Score = 1: SNP(s) overlap one feature type only
To identify genes whose genetically predicted expression levels are associated with pancreatic cancer risk, we performed Transcriptome-Wide Association Study (TWAS) analysis using FUSION pipeline.
Also, we performed SMR portal analysis after formatting GWAS sumstats via GWASLab in order to meet the input requirements.
List of scripts and their functions for Step 13:
| Script | Function |
|---|---|
| 47 | 1.Clean_sumstats |
| 47 | 2.Munge_sumstats.py |
| 47 | 3.TWAS_Expression_imputation |
| 47 | 4.TWAS_Gene_disease_association |
| 47 | 5.TWAS_Joint_conditional_test |
| 47 | 6.SMR_input |
Example commands:
# Run TWAS-FUSION (Script 47)
# QC GWAS sumstats by GWASLab
bsub < 13_TWAS_SMR/47_1.Clean_sumstats.bsub -R "rusage[mem=96GB]"
# Munge GWAS sumstas for TWAS input
bsub < 13_TWAS_SMR/47_2.Munge_sumstats.py.bsub -R "rusage[mem=48GB]"
# Run gene-disease association test
bsub < 13_TWAS_SMR/47_3.TWAS_Expression_imputation.bsub -R "rusage[mem=48GB]"
bsub < 13_TWAS_SMR/47_4.TWAS_Gene_disease_association.bsub -R "rusage[mem=48GB]"
# Run cojo/conditional test
bsub < 13_TWAS_SMR/47_5.TWAS_Joint_conditional_test.bsub -R "rusage[mem=48GB]"
# Format GWAS sumstats for SMR input
bsub < 13_TWAS_SMR/47_6.SMR_input.bsub -R "rusage[mem=30GB]"
To evaluate the predictive performance of GWAS-identified variants for pancreatic cancer risk, we developed and validated Polygenic Risk Scores (PRS) in two independent cohorts:
- All of Us Research Program (AoU) — using Jupyter notebooks with Hail/Python
- Estonian Biobank (EstBB) — using bash/PLINK2 and R scripts
PRS were calculated using different GWAS summary statistics:
- PRS-22: Based on 22 previously known PDAC risk loci
- PRS-34 (EUR): Based on 34 loci from European fixed-effects meta-analysis
- PRS-34 (MA): Based on 34 loci from multi-ancestry fixed-effects meta-analysis
- PRS-COMB (EUR): Combined score from European analyses
- PRS-COMB (TA): Combined score from trans-ancestry analyses
- Define cohort: Extract PDAC cases and cancer-free controls
- Prepare PRS weights: Format GWAS summary statistics as PRS weight tables (effect allele + beta)
- Calculate PRS: Compute weighted sum of risk alleles for each individual
- Fit logistic models: Evaluate PRS performance adjusting for age, sex, PCs, and ancestry
- Performance metrics: Calculate AUC, odds ratios by PRS quantiles, and risk stratification
Note
PRS analyses were performed in two separate cohorts with restricted data access. All of Us analyses used cloud-based Jupyter notebooks. Estonian Biobank analyses used PLINK2-based bash scripts and R for downstream evaluation.
List of scripts and their functions for Step 14:
| Script | Function | Cohort |
|---|---|---|
| 48 | Define case-control cohort and covariates (notebook) | All of Us |
| 49 | Prepare PRS weight tables and calculate PRS (notebook) | All of Us |
| 50 | Fit logistic regression models and evaluate (notebook) | All of Us |
| 51 | PRS calculation workflow (bash/PLINK2) | Estonian Biobank |
| 52 | PRS downstream analysis and visualization (R) | Estonian Biobank |
Example workflows:
# Step 1: Define cohort (in All of Us Researcher Workbench) (script 48)
# Run 48_covar_cohort.ipynb to extract:
# - Cases: pancreatic cancer (ICD-10-CM C25*)
# - Controls: cancer-free individuals
# - Covariates: age, sex, PCs, ancestry
# Step 2: Calculate PRS (in All of Us Researcher Workbench) (script 49)
# Run 49_prs_loop.ipynb to:
# - Prepare PRS weight tables for each score
# - Calculate PRS using Hail VDS and AoUPRS package
# Step 3: Model fitting and evaluation (script 50)
# Run 50_fit_model.ipynb to:
# - Fit logistic regression (PRS + covariates)
# - Calculate AUC for each PRS
# - Generate quantile-based ORs (quintiles)
# - Create ROC curves# Step 1: Calculate PRS using PLINK2
bash "51_bash estbb_prs.sh"
# This script runs PLINK2 --score for each PRS weight file
# Step 2: Downstream analysis and visualization
Rscript 52_estbb_prs_downstream.R <PRS.sscore> <phenotype.txt> <phenotype_col> <output_dir>
# Generates AUC, quantile ORs, and ROC curvesThis step generates publication-quality figures summarizing key results from multiple analyses, including PRS performance, gene prioritization, and tissue enrichment.
- PRS OR plots: Odds ratios by quantile across cohorts (EstBB, All of Us) and ancestry groups
- Gene prioritization tile plots: Heatmap summarizing evidence scores across multiple gene prioritization methods (VEP, GREAT, FLAMES, MAGMA, Open Target L2G, TWAS, SMR, Open4Gene, caABC)
- WebCSEA enrichment plots: Tissue and cell-type enrichment results from WebCSEA showing developmental stage (adult vs fetal) specificity
Note
These plotting scripts integrate results from multiple upstream analyses (Steps 1-14) to create comprehensive visualization of PDAC risk loci, candidate genes, and tissue-specific regulatory mechanisms.
List of scripts and their functions for Step 15:
| Script | Function |
|---|---|
| 53 | PRS odds ratio forest plots |
| 54 | Gene prioritization tile heatmap |
| 55 | WebCSEA tissue enrichment |
Example commands:
# PRS OR plot (EstBB + All of Us, by ancestry)
Rscript --slave --no-restore --no-save 15_plots/prs_plot.R
# Gene prioritization tile plot
Rscript --slave --no-restore --no-save 15_plots/tile_plot.R
# WebCSEA tissue enrichment plots
Rscript --slave --no-restore --no-save 15_plots/webcsea_plot.RPlot Features:
- Multi-panel comparison of PRS performance across cohorts and ancestries
- Integrated evidence scores for candidate effector genes at each locus
- Tissue and cell-type enrichment with developmental stage stratification
- Publication-ready PDF and PNG outputs
This section contains additional supporting analyses and data processing scripts used in the manuscript.
- SMR-FDR correction: Apply false discovery rate (FDR) correction to Summary-based Mendelian Randomization (SMR) results from multiple QTL types (eQTL, sQTL, mQTL, pQTL)
- Data grouping utilities: Format and clean input files for downstream tools (e.g., GREAT interval merging)
Note
SMR analyses test for pleiotropic associations between GWAS signals and molecular QTLs to identify genes whose expression/splicing/methylation/protein levels mediate disease risk. Results are filtered using FDR < 0.05 and HEIDI test p > 0.01 to exclude linkage.
List of scripts and their functions for Step 16:
| Script | Function |
|---|---|
| 51 | SMR FDR correction and filtering |
| 52 | Group and format input files |
Example commands:
# Apply FDR correction to SMR results (script 51)
Rscript --slave --no-restore --no-save 16_others/51_SMR_FDR.R
# Format grouped input files (script 52)
Rscript --slave --no-restore --no-save 16_others/52_group_by.RKey Outputs:
- FDR-corrected SMR results for EUR and multi-ancestry analyses
- Significant gene-QTL associations (FDR < 0.05, HEIDI p > 0.01)
- Summary tables of significant findings per QTL type
Note
This repository is created to enhance reproducibility. It will not be actively maintained.
Tip
If you need assistance at any step, please feel free to email me.
Important
It is not possible to share the full datasets used in this pipeline, such as UKBB data. Additionally, some steps require manual verification and the creation of data files.
Warning
This pipeline, or any part of it, should not be used for medical or diagnostic purposes.
Caution
If you use any part of this pipeline in your work, please make sure to cite the original work software/package and if possible our work.