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