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