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

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

p.cutoff = 0.01

data(lgrc.expr)
data(lgrc.expr.meta)
data(lgrc.genes)


###################################################
### code chunk number 2: lgrc_sdcd_expression.Rnw:54-68
###################################################
design.mtx = cbind(ctrl=1,
		copd=as.integer(grepl("COPD",colnames(expr))),
		age=expr.meta$age,
		pkyr=expr.meta$pkyrs)

good.idx = apply(design.mtx,1,function(x){!any(is.na(x))}) & (expr.meta$gender == "1-Male")
male.fit = lmFit(log(expr)[,good.idx], design.mtx[good.idx,])
male.fit = eBayes(male.fit)

good.idx = apply(design.mtx,1,function(x){!any(is.na(x))}) & (expr.meta$gender == "2-Female")
female.fit = lmFit(log(expr)[,good.idx], design.mtx[good.idx,])
female.fit = eBayes(female.fit)

male.female.copd.beta.diff.genes = sdcd(male.fit, female.fit, "copd", lgrc.genes, fdr.cutoff=0.25, file.prefix="male.female.copd", write.file=FALSE)


###################################################
### code chunk number 3: lgrc_sdcd_expression.Rnw:75-89
###################################################
design.mtx = cbind(ctrl=1,
		gender=expr.meta$gender,
		age=expr.meta$age,
		pkyr=expr.meta$pkyrs)

good.idx = apply(design.mtx,1,function(x){!any(is.na(x))}) & grepl("COPD",colnames(expr))
copd.fit = lmFit(log(expr)[,good.idx], design.mtx[good.idx,])
copd.fit = eBayes(copd.fit)

good.idx = apply(design.mtx,1,function(x){!any(is.na(x))}) & grepl("CTRL",colnames(expr))
ctrl.fit = lmFit(log(expr)[,good.idx], design.mtx[good.idx,])
ctrl.fit = eBayes(ctrl.fit)

copd.ctrl.gender.beta.diff.genes = sdcd(copd.fit, ctrl.fit, "gender", lgrc.genes, fdr.cutoff=0.25, file.prefix="copd.ctrl.gender", class.names=c("copd","ctrl"), write.file=FALSE)


###################################################
### code chunk number 4: lgrc_sdcd_expression.Rnw:96-102
###################################################
male.female.copd.beta.diff.genes.all = sdcd(male.fit, female.fit, "copd", lgrc.genes, fdr.cutoff=10, file.prefix="male.female.copd", write.file=FALSE)
copd.ctrl.gender.beta.diff.genes.all = sdcd(copd.fit, ctrl.fit, "gender", lgrc.genes, fdr.cutoff=10, file.prefix="copd.ctrl.gender", class.names=c("copd","ctrl"), write.file=FALSE)
all.beta.diff.genes = cbind(copd.ctrl.gender.beta.diff.genes.all, male.female.copd.beta.diff.genes.all)
rename.col = grep("beta.diff", names(all.beta.diff.genes))
names(all.beta.diff.genes)[rename.col[1:2]] = names(all.beta.diff.genes)[rename.col[1:2]] %+% ".copd.ctrl"
names(all.beta.diff.genes)[rename.col[3:4]] = names(all.beta.diff.genes)[rename.col[3:4]] %+% ".male.female"


###################################################
### code chunk number 5: lgrc_sdcd_expression.Rnw:105-107 (eval = FALSE)
###################################################
## sdcd.genes = merge(copd.ctrl.gender.beta.diff.genes, male.female.copd.beta.diff.genes, by=setdiff(intersect(names(copd.ctrl.gender.beta.diff.genes), names(male.female.copd.beta.diff.genes)),c("beta.diff","beta.diff.pooled.sd")))
## sdcd.genes = unique(sdcd.genes)


###################################################
### code chunk number 6: lgrc_sdcd_expression.Rnw:110-112
###################################################
data(lgrc.sdcd.genes)
print("There are " %+% nrow(sdcd.genes) %+% " SDCD genes")


###################################################
### code chunk number 7: lgrc_sdcd_expression.Rnw:117-125
###################################################
# FIGURE 1B
my.smart.plot(male.fit$coefficients[,"copd"], female.fit$coefficients[,"copd"], main="Coefficients of differential gene expression in males and females", xlab=expression(alpha['male']), ylab=expression(alpha['female']), colSet="Greys")
my.smart.plot(male.fit$coefficients[male.female.copd.beta.diff.genes$ensembl_gene_id,"copd"], female.fit$coefficients[male.female.copd.beta.diff.genes$ensembl_gene_id,"copd"], colSet="Blues", type="points")
my.smart.plot(male.fit$coefficients[copd.ctrl.gender.beta.diff.genes$ensembl_gene_id,"copd"], female.fit$coefficients[copd.ctrl.gender.beta.diff.genes$ensembl_gene_id,"copd"], colSet="Reds", type="points")
my.smart.plot(male.fit$coefficients[sdcd.genes$ensembl_gene_id,"copd"], female.fit$coefficients[sdcd.genes$ensembl_gene_id,"copd"], colSet="Purples", type="points")
abline(0,1,lty=2,col="gray")
abline(h=c(qnorm(0.025),qnorm(0.975)),v=c(qnorm(0.025),qnorm(0.975)),lty=3,col="gray")
smartlegend("right","bottom",c("stratified by gender","stratified by case-control status","both (SDCD)"),pch=20,col=c("blue","red","purple"))


###################################################
### code chunk number 8: lgrc_sdcd_expression.Rnw:128-143
###################################################
# FIGURE 1C
all.beta.diff.genes$copd.ctrl.beta.diff = all.beta.diff.genes$copd.beta - all.beta.diff.genes$ctrl.beta
this.pch = 20
my.smart.plot(all.beta.diff.genes$copd.ctrl.beta.diff, -log10(all.beta.diff.genes$copd.ctrl.p), main="Volcano plot for COPO-control differential expression", xlab=expression(beta['COPD'] - beta['control']), ylab=expression(-log(p['COPD,control'])), colSet="Greys", pch=this.pch)
my.smart.plot(all.beta.diff.genes[sdcd.genes$ensembl_gene_id,"copd.ctrl.beta.diff"], -log10(all.beta.diff.genes[sdcd.genes$ensembl_gene_id,"copd.ctrl.p"]), colSet="Purples", type="points", pch=this.pch)
smartlegend("right","top",c("SDCD Genes"),pch=this.pch,col=c("purple"))
CIpercent = 0.9
abline(v=quantile(all.beta.diff.genes$beta.diff.copd.ctrl, c((1-CIpercent)/2, (1+CIpercent)/2)), col="red", lty=2, lwd=1)

extreme.betas.idx = abs(sdcd.genes$beta.diff.x) > 0.25 | (abs(sdcd.genes$beta.diff.x) > 0.2 & sdcd.genes$copd.ctrl.p < 1e-5)
extreme.betas = cbind(sdcd.genes[extreme.betas.idx, c("hgnc_symbol","beta.diff.x","copd.ctrl.p","male.female.p.adj","copd.ctrl.p.adj","chromosome_name")], 
n.log.p=-log10(sdcd.genes[extreme.betas.idx, c("copd.ctrl.p")]))
print("Extreme beta_diff points are: ")
print(extreme.betas)
text(extreme.betas$beta.diff.x, extreme.betas$n.log.p, extreme.betas$hgnc_symbol, pos=1, cex=0.8)


###################################################
### code chunk number 9: lgrc_sdcd_expression.Rnw:146-154 (eval = FALSE)
###################################################
## # Figure S2
## all.beta.diff.genes$male.female.beta.diff = all.beta.diff.genes$male.beta - all.beta.diff.genes$female.beta
## this.pch = 20
## my.smart.plot(all.beta.diff.genes$male.female.beta.diff, -log10(all.beta.diff.genes$male.female.p), main="Volcano plot for male-female differential expression", xlab=expression(alpha['male'] - alpha['female']), ylab=expression(-log(p['male,female'])), colSet="Greys", pch=this.pch)
## my.smart.plot(all.beta.diff.genes[sdcd.genes$ensembl_gene_id,"male.female.beta.diff"], -log10(all.beta.diff.genes[sdcd.genes$ensembl_gene_id,"male.female.p"]), colSet="Purples", type="points", pch=this.pch)
## smartlegend("right","top",c("SDCD Genes"),pch=this.pch,col=c("purple"))
## CIpercent = 0.9
## abline(v=quantile(all.beta.diff.genes$beta.diff.male.female, c((1-CIpercent)/2, (1+CIpercent)/2)), col="red", lty=2, lwd=1)


###################################################
### code chunk number 10: lgrc_sdcd_expression.Rnw:159-160
###################################################
sessionInfo()