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