...
BMI_1kgph3_chr16_snps_summarystat_2024.txt
3.1- Base Data
We will use as the base data part of GWAS Anthropometric 2015 BMI summary statistics ( https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4382211/), made available by the GIANT consortium and were extracted from their online portal
...
Some additional codes to try (need to download and install required items following instructions):
| Code Block | ||
|---|---|---|
| ||
###### PRSice2 module load lang/R/3.5.1-Python-3.8.5-Anaconda3-2020.11 R --vanilla library(ggplot2) library(optparse) library(method) library( #### ##Some quality control steps before using summary statistics file in PRSice2 ## *** make sure to review any documents that come with the summary statistics to understand how the summary statistics were derived and what the column names stand for *** # In R: ss=read.table("BMI_1kgph3_chr16_snps_summarystat_2024.txt",header=T) nrow(ss) #[1] 68385 head(ss) summary(ss$INFO) #imputation quality summary(ss$N) #N for each SNP in source GWAS summary(ss$Freq1.Hapmap) #minor allele frequency for each SNP nrow(ss[ss$INFO<0.8,]) #[1] 199 nrow(ss[ss$N<0.5*max(ss$N),]) #[1] 2629 nrow(ss[ss$Freq1.Hapmap<=0.01 | is.na(ss$Freq1.Hapmap),]) #[1] 1590 #Some quality control for imputation quality, N, and minor allele frequency: ss1=ss[ss$INFO>0.8 & ss$N>0.5*max(ss$N) & ss$Freq1.Hapmap>0.01 & !is.na(ss$Freq1.Hapmap),] nrow(ss1) #[1] 64933 head(ss1) write.table(ss1,"BMI_1kgph3_chr16_snps_summarystat_2024_short.txt", col.names=T, row.names=F, quote=F, sep='\t') ##Now we are ready to try PRSice2! ##Need to set up the following modules in R for PRSice2 to run properly: module load lang/R/3.5.1-Python-3.8.5-Anaconda3-2020.11 R --vanilla library(ggplot2) library(optparse) library(method) library(tools) library(data.table) library(grDevices) library(RColorBrewer) q() Rscript ##This is sample script (***Need to change the actual file names***): #Rscript workdir/PRSice.R --dir workdir --prsice ./workdir/PRSice_linux --base workdir/source_gwas_sumstat_info_9_prsicein.txt --target indir/target_gwas_cleaned_plink --bar-levels 5e-8,0.00001,0.00005,0.0001,0.0005,0.001,0.005,0.01,0.05,0.1,0.2,0.3,0.4,0.5,1 --seed 1234 --perm 10000 --fastscore --all-score --no-regress T --out outdir/target_gwas_cleaned_source_gwas_prs ## PRSice2 within R on Windows (from tutorial): system("PRSice_win64.exe PRSice.R --prsice ./PRSice.exe --base TOY_BASE_GWAS.assoc --target TOY_TARGET_DATA --thread 1 --stat OR --binary-target T") |
| Code Block | ||
|---|---|---|
| ||
## PRS-CS
#need effective sample size of source GWAS (neff =2*N_cases*N_controls/(N_cases+N_controls))
python workdir/PRScs/PRScs.py --ref_dir=workdir/ldblk_1kg_eur --bim_prefix=indir/target_gwas_cleaned_plink --sst_file=workdir/source_gwas_sumstat_prscsin.txt --n_gwas=[neff] --seed=1234 --out_dir=outdir/target_gwas_source_gwas_sumstat_prscsout
# Merge sumstat outputs from PRS-CS chr1-22 in R:
R
mci<-c()
mcall<-c()
for (i in 1:22) {
mci<-read.table(paste("outdir//target_gwas_source_gwas_sumstat_prscsout_pst_eff_a1_b0.5_phiauto_chr",i,".txt",sep=""),header=F)
mcall<-rbind(mcall,mci)
}
colnames(mcall)<-c("CHR","SNP","BP","A1","A2","B")
mcalla=subset(mcall,select=c("SNP","A1","B"))
colnames(mcalla)<-c("SNP","A1","Score")
write.table(mcalla,"outdir//target_gwas_source_gwas_sumstat_prscsout_pst_eff_a1_b0.5_phiauto_chr1-22.raw", row.names=F, col.names=T, quote=F, sep='\t')
# Run polygenic scoring in PLINK:
plink --bfile indir/target_gwas_cleaned_plink --score outdir/target_gwas_source_gwas_sumstat_prscsout_pst_eff_a1_b0.5_phiauto_chr1-22.raw --out outdir/target_gwas_source_gwas_sumstat_prscsout |
...