--- a
+++ b/.Rhistory
@@ -0,0 +1,512 @@
+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()