From 00aaa4e8bed6ffed56ff4939ade4937cc2b2b8d3 Mon Sep 17 00:00:00 2001
From: MichelaTr <90008345+MichelaTr@users.noreply.github.com>
Date: Fri, 24 Nov 2023 11:35:35 -0800
Subject: [PATCH] Add files via upload

---
 GWAS_PRS/scripts/step03_plink_GWAS.sh        |  3 ++-
 GWAS_PRS/scripts/step04_GWAS_visualization.R | 15 +++++++++++----
 2 files changed, 13 insertions(+), 5 deletions(-)

diff --git a/GWAS_PRS/scripts/step03_plink_GWAS.sh b/GWAS_PRS/scripts/step03_plink_GWAS.sh
index 034d080..9f3c5f3 100644
--- a/GWAS_PRS/scripts/step03_plink_GWAS.sh
+++ b/GWAS_PRS/scripts/step03_plink_GWAS.sh
@@ -2,6 +2,7 @@
 #GWAS 
 ########################################################################
 plink_dir=~/Desktop/handson_plink
+cd ~/Desktop/handson_plink/GWAS
 mkdir plink_results/
 #linear model --glm or --linear
 #no covariates
@@ -31,7 +32,7 @@ $plink_dir/plink2 --pfile genotypes/unrelated_hg38_auto_maf_0.05_500k_snp_filter
 #logistic model --glm or --logistic
 $plink_dir/plink2 --1 --pfile genotypes/unrelated_hg38_auto_maf_0.05_500k_snp_filtered  \
 --covar    covariates/covariates_with_10PCs.csv    --covar-variance-standardize   \
---glm omit-ref cols=+beta,+orbeta,+ci hide-covar  \
+--glm omit-ref cols=+orbeta,+ci hide-covar  \
 --pheno  phenotypes/phenofile.csv  --pheno-name obesity --ci 0.95  \
 --out plink_results/sex.age.PCs.adj
 
diff --git a/GWAS_PRS/scripts/step04_GWAS_visualization.R b/GWAS_PRS/scripts/step04_GWAS_visualization.R
index a9d6e40..7f9c4aa 100644
--- a/GWAS_PRS/scripts/step04_GWAS_visualization.R
+++ b/GWAS_PRS/scripts/step04_GWAS_visualization.R
@@ -3,7 +3,7 @@
 #Visualization in R
 ##################
 #manhattan plot
-install.packages("qqman")
+#install.packages("qqman")
 library(qqman)
 
 setwd("~/Desktop/handson_plink/GWAS/")
@@ -44,15 +44,22 @@ qq(gwas_BMI_sex_age_adj$P, main="quantile-quantile plot - sex, age adjusted BMI"
 qq(gwas_BMI_sex_age_PCs_adj$P, main="quantile-quantile plot - sex, age, PCs adjusted BMI")
 
 dev.off()
-#include
-#pvals <- read.table("AMH_all_FVG_rs_pvalue", header=T)
 
-#observed <- sort(pvals$pgc)
+
+#alternative to calculate x and y axes for qqplot
+
+#observed <- sort(gwas_obesity$P)
 #lobs <- -(log10(observed))
 
 #expected <- c(1:length(observed)) 
 #lexp <- -(log10(expected / (length(expected)+1)))
 
+#pdf("R_plots/obesity_qqplot.pdf", width=480, height=480)
+#plot(c(0,8), c(0,8), col="red", lwd=3, type="l", xlab="Expected (-logP)", ylab="Observed (-logP)", main="quantile-quantile plot - sex, age, PCs adjusted obesity",  xlim=c(0,8), ylim=c(0,8), las=1, xaxs="i", yaxs="i", bty="l")
+#points(lexp, lobs, pch=23, cex=.4, bg="black") 
+#dev.off()
+
+
 
 ########################################################################
 #open and run step05_Plink_PRS.sh