PYRAMA/analysis_code.Rmd at main · pbagos/PYRAMA · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152

---
title: "PYRAMA - Analysis Code"
author: Georgios A. Manios
output:
  html_document:
    toc: true
    toc_depth: 3
    toc_float: false
    theme: cosmo           # Bootstrap theme for better visuals

---

# 1. Use case scenario (Evangelou et al. 2007)
## Data availability
To replicate the results, all data and results used are provided in this [link](https://drive.google.com/file/d/1GLbLWvpG1LVYXQoV5qaQsbmUpMerxZk7/view?usp=sharing).

## 1.1 Table 3 meta-analysis

### Imports

```{r load-libraries, message=FALSE}
library(readr)
library(dplyr)
```
### Command for PYRAMA meta-analysis to replicate Evangelou et al. Table 3 result
```
python3 pyrama.py --input_file use_case_scenario/Table3_data_results/mayo_tier1_common_in_all.txt use_case_scenario/Table3_data_results/mayo_tier2_common_in_all.txt use_case_scenario/Table3_data_results/ninds_common_in_all_harmonized.txt  --o replicate_results_table3/table3_meta_analysis_results.txt

```
### Load the meta-analysis results of PYRAMA for Table 3

```{r load-meta-analysis-data}
# Load the meta-analysis results for Table 3
table3_meta_analysis_results <- read_delim(
  "replicate_results_table3/table3_meta_analysis_results.txt",
  delim = "	",
  escape_double = FALSE,
  trim_ws = TRUE
)

```

### Load the reported significant variants of Table 3 from the orignal publication and compare them to the results


```{r filter-significant-snps}
SNPs_PD_table3 <- read_delim("use_case_scenario/Table3_data_results/SNPs_PD_table3.txt",
     delim = "\t", escape_double = FALSE,
    trim_ws = TRUE)
print(SNPs_PD_table3, n =50)


# Filter significant SNPs that are also in the PD-related SNP list
significant_snps_in_PD <- table3_meta_analysis_results[
  table3_meta_analysis_results$`P(FE)` <= 0.05 &
  table3_meta_analysis_results$SNP %in% SNPs_PD_table3$SNP, ]

# Print number of intersecting significant SNPs
cat("Number of significant SNPs in PD list:", nrow(significant_snps_in_PD), "\n")

# Display the intersected SNPs
print( significant_snps_in_PD, n =50)

```

## 1.2 Meta-analysis in all studies

### Quality control in PYRAMA  in all studies of Evangelou et al.
```
python3 quality_control.py --input_file use_case_scenario/imputation_meta_analysis/mayo_tier1.txt use_case_scenario/imputation_meta_analysis/mayo_tier2.txt  use_case_scenario/imputation_meta_analysis/ninds.txt --o qc_report.txt
```

### Meta-analysis in all studies after quality control (Harmonization etc.)
```
python3 pyrama.py --input_file  mayo_tier2_harmonized.txt  mayo_tier1_harmonized.txt  ninds_harmonized.txt --o PD_meta_analysis.txt
```

#### Read meta-analysis results and see how many SNPs where significant and were in more than one study (N>1)
```{r}
PD_meta_analysis <- read_delim("PD_meta_analysis.txt",
     delim = "\t", escape_double = FALSE,trim_ws = TRUE)

sig <- PD_meta_analysis[
  PD_meta_analysis$`P(FE)` >= min(PD_meta_analysis$`P(FE)` ) & (PD_meta_analysis$`P(FE)` <= 0.05) &
  PD_meta_analysis$N > 1,
]

print(sig)
```


## 1.3 Imputation and meta-analysis with all studies


### PRED-LD commands to perform summary statistics imputation
```
python3 pred_ld.py  --file-path mayo_tier1.txt  --r2threshold 0.5 --pop EUR --maf 0  --ref all_panels  -- imp_list mayo_tier1.txt_missing_snps_from_all_others.txt


python3 pred_ld.py  --file-path mayo_tier2.txt  --r2threshold 0.5 --pop EUR --maf 0  --ref all_panels  --imp_list mayo_tier2.txt_missing_snps_from_all_others.txt


python3 pred_ld.py  --file-path ninds.txt  --r2threshold 0.5 --pop EUR --maf 0  --ref all_panels  --imp_list ninds.txt_missing_snps_from_all_others.txt

```
### Quality control with the imputed studies
```
 python3 quality_control.py --input_file mayo_tier1.txt_imputation_results_chr_all.txt   mayo_tier2.txt_imputation_results_chr_all.txt  ninds.txt_imputation_results_chr_all.txt --o qc_report.txt
```
### Post Imputation meta-analysis
```
python3 pyrama.py --input_file  mayo_tier2_harmonized_imp.txt  mayo_tier1_harmonized_imp.txt  ninds_harmonized_imp.txt --o PD_post_imp_meta_analysis.txt
```
### Post Imputation meta-analysis results
```{r}
PD_post_imp_meta_analysis <- read_delim("PD_post_imp_meta_analysis.txt",
     delim = "\t", escape_double = FALSE,trim_ws = TRUE)

# Random Effects P-value < 1e-8  (all variants were in >=2 studies)
print (PD_post_imp_meta_analysis[PD_post_imp_meta_analysis$P<1e-8,])
```

### Enrichment analysis (g:Profiler) with the significant variants of the post-imputation meta-analysis
```{r}
gProfiler_hsapiens_22_05_2025_17_35_10_intersections <- read_csv("use_case_scenario/enrichment_analysis_results/gProfiler_hsapiens_22-05-2025_17-35-10__intersections.csv")

print(gProfiler_hsapiens_22_05_2025_17_35_10_intersections,n =300)
```

# 2. Execution time comparison with PLINK, GWAMA and METAL

This section compares execution times of PLINK, GWAMA and METAL. For this part, we repeated the simulated studies of (100k, 500k and 1 million rows) 10, 15, 20 and 25 times as input files.

```
# PYRAMA
time python3 pyrama.py --input_file  1000000_rows.txt mayo_tier1_harmonized.txt [repeated_times n times] ninds_harmonized.txt --o test.txt

# PLINK
time ./plink --meta-analysis 1000000_rows.txt [repeated_times n times]   + logscale report-all


# METAL (reads metal.txt file to load studies)
./metal  metal.txt

# GWAMA (reads GWAMA.in file to load studies)
time   ./GWAMA -qt --random -o test.txt
```

---