forked from marcoralab/MRcovid
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathformatedMan.R
More file actions
46 lines (40 loc) · 1.09 KB
/
formatedMan.R
File metadata and controls
46 lines (40 loc) · 1.09 KB
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
#!/usr/bin/env Rscript
```{r}
library(data.table)
library(dplyr)
library(qqman)
library(data.table)
getwd()
# ---- Input / Output ----
infile <- "/wynton/group/andrews/users/rakshyasharma/PRS/MR/data/formated/BelloyAll_AD/BelloyAll_AD_formated.txt.gz"
outfile <- "/wynton/group/andrews/users/rakshyasharma/PRS/MR/results/BelloyAll_AD_Manhattan.png"
# ---- Read data ----
message("Reading: ", infile)
df <- fread(infile)
df
# ---- Standardize chromosome ----
df <- df %>%
mutate(CHR = gsub("^chr", "", CHROM)) %>% # remove 'chr' prefix if exists
mutate(CHR = ifelse(CHR %in% c("X","Y","MT"), NA, as.numeric(CHR))) %>%
filter(!is.na(CHR)) %>%
rename(BP = POS)
# ---- Build qqman input ----
man_df <- df %>% select(SNP, CHR, BP, P)
# ---- Plot ----
png(outfile, width = 1600, height = 600)
manhattan(
man_df,
chr = "CHR",
bp = "BP",
snp = "SNP",
p = "P",
genomewideline = -log10(5e-8),
suggestiveline = -log10(1e-5),
cex = 0.6,
cex.axis = 0.8,
col = c("#4E79A7", "#E15759"),
main = "Belloy ALL Alzheimer's GWAS (formatted)"
)
dev.off()
message("Saved: ", outfile)
```