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

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

p.cutoff = 0.01

data(lgrc.methp)
data(lgrc.meta)
data(lgrc.sdcd.genes)

sampleID = names(methp)[grepl("^LT",names(methp),perl=TRUE) & (names(methp) %in% meta$tissueid)]
only.methp = as.matrix(methp[,sampleID])
row.names(only.methp) = methp$name
colnames(only.methp) = sampleID


###################################################
### code chunk number 2: lgrc_sdcd_methp.Rnw:57-79
###################################################
sdcd.genes = subset(sdcd.genes, chromosome_name != "HSCHR6_MHC_QBL")
sdcd.genes.bed = GRanges(seqnames=Rle("chr" %+% sdcd.genes$chromosome_name),
                         ranges=IRanges(start=sdcd.genes$start_position, 
                                       end=sdcd.genes$end_position, 
                                       names=sdcd.genes$ensembl_gene_id),
                         strand=Rle(strand(sdcd.genes$strand)))

methp = subset(methp, !(chr %in% c("chrX", "chrY")))
methp.bed = GRanges(seqnames=Rle(methp$chr),
                    ranges=IRanges(start=methp$start,
                                  end=methp$end,
                                  names=methp$name))

window = 2e4
sum(countOverlaps(sdcd.genes.bed, resize(methp.bed, window, fix="center")) != 0) # 397 SDCD genes have VMR w/in 10kb
sum(countOverlaps(sdcd.genes.bed, resize(methp.bed, 2e6)) != 0)
sum(countOverlaps(resize(methp.bed, window, fix="center"), sdcd.genes.bed) != 0) # 892 VMR have SDCD genes w/in 10kb

sdcd.genes.vmr.bed = subsetByOverlaps(resize(methp.bed, window, fix="center"), sdcd.genes.bed)
# 892 VMRs have SDCD genes

sdcd.genes.vmr = names(sdcd.genes.vmr.bed)


###################################################
### code chunk number 3: lgrc_sdcd_methp.Rnw:86-98
###################################################
design = cbind(ctrl=1,
	  	gender=as.integer(meta[sampleID,"GENDER"] == "1-Male"),
		age=meta[sampleID,"age"],
		pkyr=meta[sampleID,"pkyrs"])

good.idx = apply(design,1,function(x){!any(is.na(x))}) & meta[sampleID,"diagmaj"] == "2-COPD/Emphysema"
copd.fit = lmFit(logit(only.methp)[sdcd.genes.vmr,good.idx], design[good.idx,])
copd.fit = eBayes(copd.fit)

good.idx = apply(design,1,function(x){!any(is.na(x))}) & meta[sampleID,"diagmaj"] == "3-Control"
ctrl.fit = lmFit(logit(only.methp)[sdcd.genes.vmr,good.idx], design[good.idx,])
ctrl.fit = eBayes(ctrl.fit)


###################################################
### code chunk number 4: lgrc_sdcd_methp.Rnw:103-106
###################################################
copd.ctrl.gender.beta.diff.genes = sdcd.vmr(copd.fit, ctrl.fit, "gender", sdcd.genes, annotate=TRUE, annotate.with="genes", fdr.cutoff=0.05, file.prefix="copd.ctrl.gender", class.names=c("copd","ctrl"))

copd.ctrl.gender.beta.diff.vmr = copd.ctrl.gender.beta.diff.genes$vmr


###################################################
### code chunk number 5: lgrc_sdcd_methp.Rnw:113-134 (eval = FALSE)
###################################################
## vmr.sdcd.gene = sapply(as.character(copd.ctrl.gender.beta.diff.genes$genesymbol), function(g) { 
## 		this.vmr.genes = unlist(strsplit(g,","))
## 		this.vmr.sdcd = this.vmr.genes[which(this.vmr.genes %in% sdcd.genes$hgnc_symbol)]
## 		if (length(this.vmr.sdcd) > 1) {
## 			print(g %+% " has more than one SDCD")
## 			better.sdcd = as.character(sdcd.genes[sdcd.genes$hgnc_symbol %in% this.vmr.sdcd,"hgnc_symbol"][which.min(sdcd.genes[sdcd.genes$hgnc_symbol %in% this.vmr.sdcd,"copd.ctrl.p.adj"])])
## 			print("Keeping " %+% better.sdcd %+% " because of copd.ctrl.p.adj")
## 			this.vmr.sdcd = better.sdcd
## 		}
## 		if (length(this.vmr.sdcd) == 0) this.vmr.sdcd = NA
## 		return(this.vmr.sdcd)
## 	} )
## interesting.vmrs = copd.ctrl.gender.beta.diff.genes$vmr[vmr.sdcd.gene %in% sdcd.genes$hgnc_symbol]
## interesting.vmrs.genes = vmr.sdcd.gene[vmr.sdcd.gene %in% sdcd.genes$hgnc_symbol]
## names(interesting.vmrs.genes) = interesting.vmrs
## copd.bool = meta[sampleID,"diagmaj"] == "2-COPD/Emphysema"
## male.bool = meta[sampleID,"GENDER"] == "1-Male"
## for (ivmr in interesting.vmrs) {
## 	this.gene = interesting.vmrs.genes[ivmr]
## 	do.sdcd.boxplot(ivmr, only.methp, copd.bool, male.bool, symbol=this.gene, filename=this.gene %+% "." %+% ivmr %+% ".pdf")
## }


###################################################
### code chunk number 6: lgrc_sdcd_methp.Rnw:139-140
###################################################
sessionInfo()