lmodel =  lm(allResults[,i] ~ .^2, data = allPredictors)

sink(paste0(outputLoc,"_int_",method,"txt"))

print(summary(lmodel))

sink()  # returns output to the console

}

# go through each method but only look at the 1000 SNP runs

only1000SNPs = which(allPredictors$X.causals =="1000")

for (i in 1:length(modelNames)) {

method = modelNames[i]

# perform regression analyses for each method

lmodel =  lm(allResults[only1000SNPs,i] ~ p_current + shared_corr, data = allPredictors[only1000SNPs,])

sink(paste0(outputLoc,"_1000_",method,"txt"))

print(summary(lmodel))

sink()  # returns output to the console

}

for (i in 1:length(modelNames)) {

method = modelNames[i]

# perform regression analyses for each method

lmodel =  lm(allResults[,i] ~ p_current + shared_corr, data = allPredictors[,])

sink(paste0(outputLoc,"_all_",method,"txt"))

print(summary(lmodel))

sink()  # returns output to the console

}

sampleSizeAsN = allPredictors$sample_size

sampleSizeAsN[which(sampleSizeAsN == "_half")] = 8000

sampleSizeAsN[which(sampleSizeAsN == "_full")] = 16000

sampleSizeAsN[which(sampleSizeAsN == "_double")] = 32000

sampleSizeAsN = as.numeric(sampleSizeAsN)

# overall results

median(allResults$meta) # 0.05962532

median(allResults$subpheno) # 0.07367951

median(allResults$shaPRS) # 0.1139283

median(allResults$SMTPred) # 0.07585284

mean(allResults$meta) #  0.07426966

mean(allResults$subpheno) #  0.09124099

mean(allResults$shaPRS) #   0.1189231

mean(allResults$SMTPred) #   0.09280919

# Combined seems a bit low, check when it is better than subpheno

combined_better_than_subpheno_indices = which(allResults$meta > allResults$subpheno)

combined_better_than_subpheno = allResults[combined_better_than_subpheno_indices,]

SMTPred_better_than_shapRS_indices = which(allResults$SMTPred > allResults$shaPRS)

SMTPred_better_than_shapRS = allResults[SMTPred_better_than_shapRS_indices,]

# create heatmaps

i=1

# create rownames

coln = c()

rGs = c()

ps= c()

cors= c()

splits= c()

sample_sizes= c()

for (i in 1:length(allPredictors[,1])) {

print(i)

sample_size = gsub("_", "", allPredictors[i,7]) #remove underscores

rG= padTo_dec( allPredictors[i,2], 4)

p= padTo_dec( round(allPredictors[i,3], 2) ,4)

corre = padTo_dec( round(allPredictors[i,4], 2) ,4)

annotations = paste0("rG:",rG, " p:",p, " cor:",corre, " ",allPredictors[i,5], " split:",allPredictors[i,6], " size:",sample_size)

coln = c(coln,annotations)

splits = c(splits,allPredictors[i,6] )

ps = c(ps,p )

cors = c(cors,corre )

sample_sizes = c(sample_sizes,sample_size )

rGs = c(rGs,paste0(rG," ") ) # add an extra space for padding, otherwise legend will be cut off as pheatmap is shit

}

rownames(allResults) = coln

# want to sort the data rows by rG

orderByRG = order( as.numeric(rGs))

allResults = allResults[orderByRG,]

splits = splits[orderByRG]

rGs = rGs[orderByRG]

ps = ps[orderByRG]

cors = cors[orderByRG]

sample_sizes = sample_sizes[orderByRG]

rGs_DF <- data.frame( rGs,as.numeric(ps),as.numeric(cors),sample_sizes, splits,row.names=rownames(allResults)) # need to match the rownames to the data for the annotation to work

colnames(rGs_DF) <- c("rG", "p", "cor","N", "split") # this is the header for the annotation

# function to separate regular/extra results:

filterOutByTerm_all = function(allResults,splits,ps,cors,filterTerm = "rG:0.50") {

subsetResults = allResults

ps_subset = ps

cors_subset = cors

splits_subset = splits

indices_kept = c()

i=1

for (i in 1:nrow(subsetResults) ) {

rowname = rownames(subsetResults[i,])

# check if rowname includes extra/regular,

if ( grepl( filterTerm, rowname, fixed = TRUE) ) {

# if yes, we replace it with nothing, and keep it

rowname_new = gsub(filterTerm, "",rowname) #remove underscores

rownames(subsetResults)[rownames(subsetResults) == rowname] <- rowname_new

indices_kept = c(indices_kept, i)

} # discard it otherwise

}

subsetResults = subsetResults[indices_kept,]

# function to separate regular/extra results:

filterOutByTerm_all = function(allResults,splits,ps,cors,filterTerm = "rG:0.50") {

subsetResults = allResults

ps_subset = ps

cors_subset = cors

splits_subset = splits

indices_kept = c()

i=1

for (i in 1:nrow(subsetResults) ) {

rowname = rownames(subsetResults[i,])

# check if rowname includes extra/regular,

if ( grepl( filterTerm, rowname, fixed = TRUE) ) {

# if yes, we replace it with nothing, and keep it

rowname_new = gsub(filterTerm, "",rowname) #remove underscores

rownames(subsetResults)[rownames(subsetResults) == rowname] <- rowname_new

indices_kept = c(indices_kept, i)

} # discard it otherwise

}

subsetResults = subsetResults[indices_kept,]

ps_subset = ps_subset[indices_kept]

cors_subset = cors_subset[indices_kept]

splits_subset = splits_subset[indices_kept]

results = NULL

results$subsetResults = subsetResults

results$ps_subset = ps_subset

results$cors_subset = cors_subset

results$splits_subset = splits_subset

return(results)

}

)

# function to separate regular/extra results:

filterOutByTerm_all = function(allResults,splits,ps,cors,filterTerm = "rG:0.50") {

subsetResults = allResults

ps_subset = ps

cors_subset = cors

splits_subset = splits

indices_kept = c()

i=1

for (i in 1:nrow(subsetResults) ) {

rowname = rownames(subsetResults[i,])

# check if rowname includes extra/regular,

if ( grepl( filterTerm, rowname, fixed = TRUE) ) {

# if yes, we replace it with nothing, and keep it

rowname_new = gsub(filterTerm, "",rowname) #remove underscores

rownames(subsetResults)[rownames(subsetResults) == rowname] <- rowname_new

indices_kept = c(indices_kept, i)

} # discard it otherwise

}

subsetResults = subsetResults[indices_kept,]

ps_subset = ps_subset[indices_kept]

cors_subset = cors_subset[indices_kept]

splits_subset = splits_subset[indices_kept]

results = NULL

results$subsetResults = subsetResults

results$ps_subset = ps_subset

results$cors_subset = cors_subset

results$splits_subset = splits_subset

return(results)

}

# Filter to keep the main interesting scenarios, rG 0.5, regular, full

results_RG05 = filterOutByTerm_all(allResults,splits,ps,cors,filterTerm = "rG:0.50")

results_regular = filterOutByTerm_all(results_RG05$subsetResults,results_RG05$splits_subset,results_RG05$ps_subset,results_RG05$cors_subset,filterTerm = "regular")

results_full = filterOutByTerm_all(results_regular$subsetResults,results_regular$splits_subset,results_regular$ps_subset,results_regular$cors_subset,filterTerm = "size:full")

subsetResults = results_full$subsetResults

subset_DF <- data.frame( results_full$ps_subset, results_full$cors_subset,results_full$splits_subset,row.names=rownames(subsetResults)) # ,row.names=rownames(subsetResults) # need to match the rownames to the data for the annotation to work

colnames(subset_DF) <- c("p","cor","split") # this is the header for the annotation

#plotName = "shaPRS - rG:0.5, n:full, no extra" # no plotname for final publication

plotName =""

pheatmap(subsetResults, main = plotName , filename=paste(outputLoc,"_subset.png", sep="" ),annotation_row=subset_DF, show_rownames = F, height=5, width=5 , cex=1 ,cluster_rows=F, cluster_cols=F)

# Filter to keep the main interesting scenarios, rG 0.5, regular, full

results_RG05 = filterOutByTerm_all(allResults,splits,ps,cors,filterTerm = "rG:0.50")

results_regular = filterOutByTerm_all(results_RG05$subsetResults,results_RG05$splits_subset,results_RG05$ps_subset,results_RG05$cors_subset,filterTerm = "regular")

results_full = filterOutByTerm_all(results_regular$subsetResults,results_regular$splits_subset,results_regular$ps_subset,results_regular$cors_subset,filterTerm = "size:full")

subsetResults = results_full$subsetResults

subset_DF <- data.frame( results_full$ps_subset, results_full$cors_subset,results_full$splits_subset,row.names=rownames(subsetResults)) # ,row.names=rownames(subsetResults) # need to match the rownames to the data for the annotation to work

colnames(subset_DF) <- c("p","cor","split") # this is the header for the annotation

plotName =""

###

# Filter to keep the main interesting scenarios, rG 0.5, extra, full

results_RG05 = filterOutByTerm_all(allResults,splits,ps,cors,filterTerm = "rG:0.50")

results_regular = filterOutByTerm_all(results_RG05$subsetResults,results_RG05$splits_subset,results_RG05$ps_subset,results_RG05$cors_subset,filterTerm = "extra")

results_full = filterOutByTerm_all(results_regular$subsetResults,results_regular$splits_subset,results_regular$ps_subset,results_regular$cors_subset,filterTerm = "size:full")

subsetResults_extra = results_full$subsetResults

subset_DF_extra <- data.frame( results_full$ps_subset, results_full$cors_subset,results_full$splits_subset,row.names=rownames(subsetResults_extra)) # ,row.names=rownames(subsetResults_extra) # need to match the rownames to the data for the annotation to work

colnames(subset_DF_extra) <- c("p","cor","split") # this is the header for the annotation

#plotName = "shaPRS - rG:0.5, n:full, no extra" # no plotname for final publication

# make pheatmap on the same colour scale:

Breaks <- seq(min(c(subsetResults, subsetResults_extra)), max(c(subsetResults, subsetResults_extra)), length = 100)

# make pheatmap on the same colour scale:

Breaks <- seq(min(subsetResults, subsetResults_extra), max(subsetResults, subsetResults_extra), length = 100)

Breaks

paste(outputLoc,"_subset.png", sep="" )

pheatmap(subsetResults, breaks = Breaks, main = plotName , filename=paste(outputLoc,"_subset.png", sep="" ),annotation_row=subset_DF, show_rownames = F, height=5, width=5 , cex=1 ,cluster_rows=F, cluster_cols=F)

pheatmap(subsetResults_extra, breaks = Breaks, main = plotName , filename=paste(outputLoc,"_subset_extra.png", sep="" ),annotation_row=subset_DF_extra, show_rownames = F, height=5, width=5 , cex=1 ,cluster_rows=F, cluster_cols=F)

args=vector()

args =c(args,"#causals_1000_rG_0.1_A50_B50_size_half")

args = c(args,"0.1")

args = c(args,"C:/softwares/Cluster/0shaPRS/debug/#causals_1000_rG_0.1_A50_B50_size_half")

args = c(args,"C:/softwares/Cluster/0shaPRS/debug/1000/10/0.1_1.0/A50_B50/size_half/")

args = c(args,"0.1")

args = c(args,"1.0")

args = c(args,"C:/softwares/Cluster/0shaPRS/debug/1000/10/0.55_0.1818182/A50_B50/size_half/")

args = c(args,"0.55")

args = c(args,"0.1818182")

args = c(args, "C:/softwares/Cluster/0shaPRS/debug/1000/10/1.0_0.1/A50_B50/size_half/")

args = c(args,"1.0")

args = c(args,"0.1")

limitsEnabled = F

# load input files for each method

combined= NULL

subpheno= NULL

shaPRS= NULL

SMTPred= NULL

xlabels=vector()

for(i in seq(from=4, to=length(args), by=4)){ # 4 as we also add the 'regular'

baseLoc=args[i]

print(paste0("baseLoc is: ", baseLoc))

current_combined = read.table(paste0(baseLoc,"combined")  ,header=F)

current_subpheno = read.table(paste0(baseLoc,"subpheno")  ,header=F)

current_shaPRS = read.table(paste0(baseLoc,"shaPRS_meta")  ,header=F)

current_SMTPred = read.table(paste0(baseLoc,"SMTPred")  ,header=F)

# replace NAs with col mean

current_combined[is.na(current_combined[,1]), 1] <- mean(current_combined[,1], na.rm = TRUE)

current_subpheno[is.na(current_subpheno[,1]), 1] <- mean(current_subpheno[,1], na.rm = TRUE)

current_shaPRS[is.na(current_shaPRS[,1]), 1] <- mean(current_shaPRS[,1], na.rm = TRUE)

current_SMTPred[is.na(current_SMTPred[,1]), 1] <- mean(current_SMTPred[,1], na.rm = TRUE)

if (is.null(combined)) {

combined= current_combined

subpheno= current_subpheno

shaPRS= current_shaPRS

SMTPred= current_SMTPred

} else {

combined= cbind( combined,current_combined )

subpheno= cbind( subpheno, current_subpheno)

shaPRS= cbind(shaPRS, current_shaPRS)

SMTPred= cbind( SMTPred, current_SMTPred)

}

p_current = round(as.numeric(args[(i+1)]),2)

shared_corr = round(as.numeric(args[(i+2)]),2)

print(paste0("p_current: ",p_current, " / shared_corr: ", shared_corr, " | baseLoc: ", baseLoc))

xlabels = c(xlabels, paste0("p:",p_current,"/r:",shared_corr) )

}

?plot

library(qvalue)

?qvalue_truncp

?qvalue::qvalue_truncp

?qvalue

inputDataLoc="C:/0Datasets/shaPRS/DEL/EUR_JAP_T2D_SE_meta"

blendFactorsLoc="C:/0Datasets/shaPRS/DEL/EUR_JAP_T2D_lFDR_meta_SNP_lFDR"

outputLoc="C:/0Datasets/shaPRS/DEL/QvalManhattan"

B12Loc = "C:/0Datasets/shaPRS/DEL/EUR_JAP_T2D_sumstats_meta"

plotTitle="Uga"

# 1. load data

inputData= read.table(inputDataLoc, header = T)

blendFactors= read.table(blendFactorsLoc, header = T)

B12= read.table(B12Loc, header = T)

inputData_blendFactors = merge(inputData,blendFactors, by ="SNP") # merge to make sure they are aligned

B12_blendFactors = merge(B12,blendFactors, by ="SNP") # merge to make sure they are aligned

# MANHATTAN PLOT

inputData_blendFactors$P=inputData_blendFactors$Qval # add the adjusted Q vals as 'P', as that is col I would be plotting next

lfdr_2_blending = 1-inputData_blendFactors$lFDR

base_colour1 = 0.20

base_colour2 = 0.20

manhattanBaseColours = c(rgb(0,base_colour2,0,1),rgb(0,0,base_colour1,1) )

allIndices = inputData_blendFactors$BP

# need to offset the SNP BP indices, by the previous number of indices in all previous chromosomes

inputData_blendFactors$Offsets = 0

for (i in 1:21) { # we always set the next offset, so we dont loop til last Chrom

message(i)

CHR_SNPs = inputData_blendFactors[inputData_blendFactors$CHR == i,]

maxBPCurrentChrom = max(CHR_SNPs$BP)

currentOffset = CHR_SNPs$Offsets[1]

nextOffset = currentOffset + maxBPCurrentChrom

inputData_blendFactors[inputData_blendFactors$CHR == (i+1),9] = nextOffset

}

hist(B12_blendFactors$lFDR, probability = T, col ="red", xlab = "lFDR", main ="")

plot(B12_blendFactors$lFDR, B12_blendFactors$b, col ="red", xlab = "lFDR", ylab = "SNP coef", main ="")

# get distribution

library(scales)

plot(B12_blendFactors$lFDR[1:5000], B12_blendFactors$b[1:5000], col =alpha("red", 0.4), xlab = "lFDR", ylab = "SNP coef", main ="")

plot(B12_blendFactors$lFDR[1:5000], B12_blendFactors$b[1:5000], col =alpha("red", 0.3), xlab = "lFDR", ylab = "SNP coef", main ="")

baseLoc="C:/0Datasets/ukbb/fis/"

phenoLoc=paste0(baseLoc, "pheno")

covarsLoc=paste0(baseLoc, "covariates")

pheno=read.table(phenoLoc  ,header=F)

View(pheno)

covars=read.table(covarsLoc  ,header=T)

View(covars)

men_index = which(covars$SEX =="M")

female_pheno = pheno[female_index]

female_index = which(covars$SEX =="F")

male_pheno = pheno[men_index]

men_index = which(covars$SEX =="M")

female_index = which(covars$SEX =="F")

female_pheno = pheno[female_index]

male_pheno = pheno[men_index]

men_index

female_pheno = pheno$V1[female_index]

male_pheno = pheno$V1[men_index]

mean(female_pheno)

mean(male_pheno) # 5.788613

t.test(male_pheno,female_pheno)

var(female_pheno)

var(male_pheno) #

percDifference = function (data1, data2) {

mean1 = mean(data1)

mean2 = mean(data2)

percDiff = round( (mean1 - mean2) / ( (mean1 + mean2)/2 ) * 100)

return(percDiff)

}

percDifference = function (data1, data2) {

percDiff = round( (data1 - data2) / ( (data1 + data2)/2 ) * 100)

return(percDiff)

}

percDifference( mean(male_pheno), mean(male_pheno) )

percDifference( mean(male_pheno), mean(female_pheno) )

female_townsend = covars$TOWNSEND[female_index]

male_townsend = covars$TOWNSEND[men_index]

mean(female_townsend)

mean(male_townsend) # -1.655731

var(female_townsend)

var(male_townsend)

t.test(male_townsend, female_townsend)

percDifference( var(male_pheno), var(female_pheno) ) # FIS: 4

baseLoc="C:/0Datasets/ukbb/height/"

phenoLoc=paste0(baseLoc, "pheno")

covarsLoc=paste0(baseLoc, "covariates")

pheno=read.table(phenoLoc  ,header=F)

covars=read.table(covarsLoc  ,header=T)

men_index = which(covars$SEX =="M")

female_index = which(covars$SEX =="F")

female_pheno = pheno$V1[female_index]

male_pheno = pheno$V1[men_index]

mean(female_pheno) # FIS: 5.788613

mean(male_pheno) # FIS: 5.995717

percDifference( mean(male_pheno), mean(female_pheno) ) # FIS: 4

t.test(male_pheno,female_pheno) # t = 22.523, df = 184654, p-value < 2.2e-16

var(female_pheno) # FIS: 3.778955

var(male_pheno) # FIS: 4.301227

percDifference( var(male_pheno), var(female_pheno) ) # FIS: 13

baseLoc="C:/0Datasets/ukbb/bmi/"

phenoLoc=paste0(baseLoc, "pheno")

covarsLoc=paste0(baseLoc, "covariates")

pheno=read.table(phenoLoc  ,header=F)

covars=read.table(covarsLoc  ,header=T)

men_index = which(covars$SEX =="M")

female_index = which(covars$SEX =="F")

female_pheno = pheno$V1[female_index]

male_pheno = pheno$V1[men_index]

mean(female_pheno) # FIS: 5.788613 # 162.6658

mean(male_pheno) # FIS: 5.995717 # 175.9021

percDifference( mean(male_pheno), mean(female_pheno) ) # FIS: 4, Height: 8

t.test(male_pheno,female_pheno) # FIS/Height: p-value < 2.2e-16

var(female_pheno) # FIS: 3.778955, Height: 38.63668

var(male_pheno) # FIS: 4.301227, 45.66119

percDifference( var(male_pheno), var(female_pheno) ) # FIS: 13, Height: 17

var(female_townsend) # 7.647166

var(male_townsend) # 7.98267 # variance is greater for men

mean(female_townsend) # -1.655731

# Compare male female townsends

female_townsend = covars$TOWNSEND[female_index]

male_townsend = covars$TOWNSEND[men_index]

mean(female_townsend) # -1.655731

mean(male_townsend) # -1.658059

t.test(male_townsend, female_townsend) #  p-value = 0.8553

var(female_townsend) # 7.647166

var(male_townsend) # 7.98267 # variance is greater for men

hist(male_pheno)

hist(female_pheno)

hist(male_pheno)

hist(female_pheno)

writeLines('PATH="${RTOOLS40_HOME}\\usr\\bin;${PATH}"', con = "~/.Renviron")

Sys.which("make")

library("devtools")

library(roxygen2)

# 1) set working Dir

#setwd("C:/Users/mk23/GoogleDrive_phd/PHD/!Publications/shaPRS/R_package_303021/shaPRS")

setwd("C:/Users/mk23/GoogleDrive_Cam/0Publications/shaPRS/R_package/shaPRS")

#install_github("jdstorey/qvalue")

#install_github("mkelcb/shaprs")

usethis::use_package("qvalue")

usethis::use_package("Matrix")

usethis::use_package("compiler")

document()

check(cran=TRUE)

document()

check(cran=TRUE)

document()

check(cran=TRUE)

document()

check(cran=TRUE)

# 1. load phenos

subphenoLoc='inst/extdata/phenoA_sumstats'

CombinedPhenoLoc='inst/extdata/Combined_sumstats'

blendFactorLoc='inst/extdata/myOutput_SNP_lFDR'

subpheno= read.table(subphenoLoc, header = T)

CombinedPheno= read.table(CombinedPhenoLoc, header = T)

blendingFactors= read.table(blendFactorLoc, header = F)

# 1. load phenos

subphenoLoc='inst/extdata/phenoA_sumstats'

subpheno_otherLoc='inst/extdata/phenoB_sumstats'

blendFactorLoc='inst/extdata/myOutput_SNP_lFDR'

subpheno= read.table(subphenoLoc, header = T)

subpheno_other= read.table(subpheno_otherLoc, header = T)

blendingFactors= read.table(blendFactorLoc, header = F)

View(subpheno)

View(subpheno_other)

View(subpheno)

View(blendingFactors)

blendingFactors= read.table(blendFactorLoc, header = T)

View(blendingFactors)

# 1.  Merge first the 3 tables together by RSid, so they are always aligned, x = subpheno  and    y = CombinedPheno ( ensure that when we check allele alignment we are comparing the same SNPs

subpheno_otherPheno = merge(subpheno,subpheno_other,by.x = "SNP",by.y = "SNP")

subpheno_otherPheno_blending = merge(subpheno_otherPheno,blendingFactors, by.x = "SNP", by.y = "SNP")

document()

check(cran=TRUE)

# Build package

build()

# Test install

install()

#install_github("jdstorey/qvalue")

install_github("mkelcb/shaprs")

library("devtools")

# Test install

install()

library("shaPRS")

# II) tests

?shaPRS_adjust

inputDataLoc <- system.file("extdata", "shapersToydata.txt", package = "shaPRS")

inputData= read.table(inputDataLoc, header = T)

results = shaPRS_adjust(inputData, thresholds=c(0.5,0.99))

# Test LD ref blend

?LDRefBlend

sumstatsData = readRDS(file = system.file("extdata", "sumstatsData_toy.rds", package = "shaPRS") )

# read SNP map files ( same toy data for the example)

pop1_map_rds = readRDS(file = system.file("extdata", "my_data.rds", package = "shaPRS") )

pop2_map_rds = readRDS(file = system.file("extdata", "my_data2.rds", package = "shaPRS") )

# use chrom 21 as an example

chromNum=21

# load the two chromosomes from each population ( same toy data for the example)

pop1LDmatrix = readRDS(file = system.file("extdata", "LDref.rds", package = "shaPRS") )

pop2LDmatrix = readRDS(file = system.file("extdata", "LDref2.rds", package = "shaPRS") )

# 2. grab the RSids from the map for the SNPS on this chrom,

# each LD mat has a potentiall different subset of SNPs

# this is guaranteed to be the same order as the pop1LDmatrix

pop1_chrom_SNPs = pop1_map_rds[ which(pop1_map_rds$chr == chromNum),]

# this is guaranteed to be the same order as the pop2LDmatrix

pop2_chrom_SNPs = pop2_map_rds[ which(pop2_map_rds$chr == chromNum),]

pop1_chrom_SNPs$pop1_id = 1:nrow(pop1_chrom_SNPs)

pop2_chrom_SNPs$pop2_id = 1:nrow(pop2_chrom_SNPs)

# intersect the 2 SNP lists so that we only use the ones common to both LD matrices by merging them

chrom_SNPs_df  <- merge(pop1_chrom_SNPs,pop2_chrom_SNPs, by = "rsid")

# align the two LD matrices

chrom_SNPs_df = alignStrands(chrom_SNPs_df, A1.x ="a1.x", A2.x ="a0.x", A1.y ="a1.y", A2.y ="a0.y")

# align the summary for phe A and B

sumstatsData = alignStrands(sumstatsData)

# subset sumstats data to the same chrom

sumstatsData = sumstatsData[which(sumstatsData$CHR == chromNum ),]

# merge sumstats with common LD map data

sumstatsData  <- merge(chrom_SNPs_df,sumstatsData, by.x="rsid", by.y = "SNP")

# remove duplicates

sumstatsData = sumstatsData[ !duplicated(sumstatsData$rsid) ,]

# use the effect alleles for the sumstats data with the effect allele of the LD mat

# as we are aligning the LD mats against each other, not against the summary stats

# we only use the lFDR /SE from the sumstats,

# which are directionless, so those dont need to be aligned

sumstatsData$A1.x =sumstatsData$a1.x

sumstatsData$A1.y =sumstatsData$a1.y

# make sure the sumstats is ordered the same way as the LD matrix:

sumstatsData = sumstatsData[order(sumstatsData$pop1_id), ]

# subset the LD matrices to the SNPs we actualy have

pop1LDmatrix = pop1LDmatrix[sumstatsData$pop1_id,sumstatsData$pop1_id]

pop2LDmatrix = pop2LDmatrix[sumstatsData$pop2_id,sumstatsData$pop2_id]

# generate the blended LD matrix

cormat = LDRefBlend(pop1LDmatrix,pop2LDmatrix, sumstatsData)

View(cormat)

uninstall()

# II) tests

?shaPRS_adjust

results$lFDRTable

# Test installing from remote

install_github("mkelcb/shaprs")

library("shaPRS")

# Test blend

?shaPRS_blend_overlap

# Test blend

?shaPRS_blend_overlap

# Test blend

?shaPRS_blend_overlap

subphenoLoc <- system.file("extdata", "phenoA_sumstats", package = "shaPRS")

subpheno_otherLoc <- system.file("extdata", "phenoB_sumstats", package = "shaPRS")

blendFactorLoc <- system.file("extdata", "myOutput_SNP_lFDR", package = "shaPRS")

subpheno= read.table(subphenoLoc, header = TRUE)

subpheno_other= read.table(subpheno_otherLoc, header = TRUE)

blendingFactors= read.table(blendFactorLoc, header = TRUE)

blendedSumstats = shaPRS_blend_overlap(subpheno, subpheno_other, blendingFactors)

View(blendedSumstats)

typeof(cormat)

map_rds_new = pop1_map_rds[which(pop1_map_rds$chr == chromNum),]

map_rds_new2 = map_rds_new[which(map_rds_new$rsid %in% sumstatsData$rsid),] # match the first to the second

View(map_rds_new2)

document()

check(cran=TRUE)

# Build package

build()

# Test installing from remote

install_github("mkelcb/shaprs")

uninstall()

# Test installing from remote

install_github("mkelcb/shaprs")

library("shaPRS")

# Test LD ref blend

?LDRefBlend

LDRefBlend

# Test LD ref blend

?LDRefBlend

uninstall()