### R code from vignette source 'vignettes/COPDSexualDimorphism/inst/doc/lgrc_sdcd_eQTL.Rnw'

###################################################
### code chunk number 1: lgrc_sdcd_eQTL.Rnw:33-35
###################################################
library(COPDSexualDimorphism)
`%+%` <- function(x,y) paste(x,y,sep="")


###################################################
### code chunk number 2: lgrc_sdcd_eQTL.Rnw:58-68
###################################################
data(lgrc.eqtl)
dim(eqtl)

print("There are " %+% length(unique(eqtl$SNP)) %+% " cis SNPs of SDCD genes.")

fdr.cutoff = 0.05
eqtl$FDR_male = p.adjust(eqtl$P_male, "BH")
eqtl$FDR_female = p.adjust(eqtl$P_female, "BH")
print(sum(eqtl$FDR_male < fdr.cutoff, na.rm=T) %+% " male, " %+% sum(eqtl$FDR_female < fdr.cutoff, na.rm=T) %+% " female, " %+% sum(eqtl$FDR_male < fdr.cutoff & eqtl$FDR_female < fdr.cutoff, na.rm=T) %+% " both.")
fisher.test(eqtl$FDR_male < fdr.cutoff, eqtl$FDR_female < fdr.cutoff)


###################################################
### code chunk number 3: lgrc_sdcd_eQTL.Rnw:73-76
###################################################
discord.ref.allele = which(eqtl$A1_male != eqtl$A1_female)
eqtl$STAT_female[discord.ref.allele] = -eqtl$STAT_female[discord.ref.allele]
eqtl$BETA_female[discord.ref.allele] = -eqtl$BETA_female[discord.ref.allele]


###################################################
### code chunk number 4: lgrc_sdcd_eQTL.Rnw:81-104
###################################################
# package the info as limma fit object to pass to sdcd.core
eqtl.male = list(
  coefficients = data.frame(copd=eqtl$BETA_male),
  stdev.unscaled = data.frame(copd=eqtl$BETA_male/eqtl$STAT_male),
  sigma = 1,
  df.residual = eqtl$NMISS_male - 4, 
  df.prior = eqtl$NMISS_male - 4
  )
eqtl.female = list(
  coefficients = data.frame(copd=eqtl$BETA_female),
  stdev.unscaled = data.frame(copd=eqtl$BETA_female/eqtl$STAT_female),
  sigma = 1,
  df.residual = eqtl$NMISS_female - 4, 
  df.prior = eqtl$NMISS_female - 4
  )
# The SDCD analysis
eqtl.sdcd = sdcd.core(eqtl.male, eqtl.female, "copd")
eqtl = cbind(eqtl, eqtl.sdcd)

all.eqtl = eqtl
eqtl = subset(eqtl, beta.diff.pval.adj < fdr.cutoff & !is.na(beta.diff.pval.adj))

print("Male-female difference: " %+% nrow(eqtl) %+% " eQTL are significant at level " %+% fdr.cutoff %+% ", covering " %+% length(unique(eqtl$Ensembl_Gene)) %+% " genes.")


###################################################
### code chunk number 5: lgrc_sdcd_eQTL.Rnw:111-112
###################################################
sessionInfo()