# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .r
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: R
#     language: R
#     name: ir
# ---

#install.packages('bigsnpr')
library(bigsnpr)

# +
system("dx pwd",intern=TRUE)   # View current directory in the DNAnexus project

# Create a test file and read it, as it were from bash terminal
# system("dx download PRS_reporting:/Input/geno/chr1.*",intern=TRUE)
# snp_readBed("chr1.bed")
# snp_readBed("PRS_reporting:/Input/geno/chr1.bed")

system("dx ls",intern=TRUE)   # View current directory in the DNAnexus project
# -

# Attach the "bigSNP" object in R session
obj.bigSNP <- snp_attach("chr1.rds")
# See how the file looks like
str(obj.bigSNP, max.level = 2, strict.width = "cut")

# Get aliases for useful slots
G   <- obj.bigSNP$genotypes
CHR <- obj.bigSNP$map$chromosome
POS <- obj.bigSNP$map$physical.pos
y   <- obj.bigSNP$fam$affection - 1
NCORES <- nb_cores()
# Check some counts for the 10 first variants
G <- G$copy(code=c(0,1,2,0, rep(NA,252)))
big_counts(G, ind.col = 1:10)

# +
# system("dx download PRS_reporting:/Input/gwas_catalog/trytrying_again_chr1.tsv.gz",intern=TRUE)
# -

sumstats <- read.table(gzfile("trytrying_again_chr1.tsv.gz"), stringsAsFactors=FALSE, header=TRUE)
str(sumstats)

names(sumstats) <- c("rsid", "chr", "pos", "a1", "a0", "beta", "p")
sumstats$p <- as.numeric(sumstats$p)
sumstats$pos <- as.numeric(sumstats$pos)
map <- obj.bigSNP$map[,-(2:3)]
head(map)
head(sumstats)
names(map) <- c("chr", "pos", "a0", "a1")
info_snp <- snp_match(sumstats, map, match.min.prop = 0)

info_snp

beta <- rep(NA, ncol(G))
beta[info_snp$`_NUM_ID_`] <- info_snp$beta
lpval <- rep(NA, ncol(G))
lpval[info_snp$`_NUM_ID_`] <- -log10(info_snp$p)

# The clumping step might take some time to complete
all_keep <- snp_grid_clumping(G, CHR, POS, ind.row = 1:length(y),
                              lpS = lpval, exclude = which(is.na(lpval)),
                              ncores = NCORES)
attr(all_keep, "grid")

multi_PRS <- snp_grid_PRS(G, all_keep, beta, lpval, ind.row = 1:length(y),
                          backingfile = "tmp-data/tmp", 
                          n_thr_lpS = 50, ncores = NCORES)
dim(multi_PRS)  ## 1400 C+T scores for 100000 individuals

ID <- obj.bigSNP$fam$family.ID
PRS_results <- data.frame(FID=ID, PID=ID)

prs_ <- multi_PRS[]
PRS_scores <- as.data.frame(prs_)

dim(PRS_results)
dim(PRS_scores)
final_output <- cbind(PRS_results, PRS_scores)
dim(final_output)
write.csv(final_output, 'PRS.csv', row.names=FALSE)

StarcoderdataJupyterScriptsDedupFiltered