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()