RELIEF: methylation study
Mark Ziemann
2020-08-18
- Introduction
- Analysis pipeline
- Obtaining array annotations
- Obtain reactome gene sets
- Curate metadata
- Import, normalise and filter
- plot mds
- extract probes on chrX
- Differential analysis including sex chromosomes
- Differential analysis excluding sex chromosomes
- Normalise for blood composition and repeat the analysis
- plot mds
- Differential analysis including sex chromosomes
- Differential analysis excluding sex chromosomes
Introduction
Analysis pipeline
In this analysis, I will use the following resources:
https://bioconductor.org/packages/release/bioc/vignettes/missMethyl/inst/doc/missMethyl.html
https://www.bioconductor.org/help/course-materials/2015/BioC2015/methylation450k.html
This script should be run with R 4.0
suppressPackageStartupMessages({
library("missMethyl")
library("limma")
library("minfi")
library("IlluminaHumanMethylationEPICanno.ilm10b2.hg19")
library("IlluminaHumanMethylationEPICmanifest")
library("beeswarm")
library("gplots")
library("topconfects")
library("FlowSorted.Blood.EPIC")
library("mitch")
})
data(Locations)
source("epic_850k_functions.R")Obtaining array annotations
anno <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
myann <- data.frame(anno[,c("UCSC_RefGene_Name","Regulatory_Feature_Group")])Obtain reactome gene sets
download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip",
destfile="ReactomePathways.gmt.zip")
unzip("ReactomePathways.gmt.zip",overwrite = TRUE)
genesets <- gmt_import("ReactomePathways.gmt")Curate metadata
tmp<-readLines("ILMLEPIC-15830_SampleSheet.csv")
sample.annotation<-read.csv(text=tmp[8:length(tmp)],stringsAsFactors=FALSE)
sample.annotation[grep("Pre",sample.annotation$Sample_Name),4]<-"NA"
sample.annotation[grep("T0",sample.annotation$Sample_Name),4]<-"Group1"
sample.annotation[grep("POD",sample.annotation$Sample_Name),4]<-"Group2"
sample.annotation$RG_number<-sapply(strsplit(sample.annotation$Sample_Name,"-"),"[",1)
# clinical metadata
samplesheet<-read.csv("samplesheet.tsv",sep="\t",stringsAsFactors=FALSE)
# merge these metadata together
sample.annotation2<-merge(sample.annotation,samplesheet,by="RG_number")
sample.annotation2$Basename<-sample.annotation2$Sentrix_ID
sample.annotation2$Basename<-paste(sample.annotation2$Sentrix_ID,sample.annotation2$Sentrix_Position,sep="_")Import, normalise and filter
basedir = "analysis/idat"
rgSet <- read.metharray.exp(basedir, targets = sample.annotation2)## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls# normalise for differences in cell type composition
mSet<-preprocessRaw(rgSet)
mSetSw <- SWAN(mSet,verbose=FALSE)## [SWAN] Preparing normalization subset
## EPIC
## [SWAN] Normalizing methylated channel
## [SWAN] Normalizing unmethylated channel# plots to show differences due to normalisation
par(mfrow=c(1,2), cex=1.25)
densityByProbeType(mSet[,1], main = "Raw")
densityByProbeType(mSetSw[,1], main = "SWAN")
detP <- detectionP(rgSet)
keep <- rowSums(detP < 0.01) == ncol(rgSet)
mSetSw <- mSetSw[keep,]
mset_reduced <- mSetSw
meth <- getMeth(mset_reduced)
unmeth <- getUnmeth(mset_reduced)
Mval <- log2((meth + 100)/(unmeth + 100))
beta <- getBeta(mset_reduced)
dim(Mval)## [1] 857902 96plot mds
plotMDS(Mval, labels=sample.annotation2$Sample_Name, col=as.integer(factor(sample.annotation2$CrpGroup)))
legend("topleft",legend=c("CRP high","CRP norm"),pch=16,cex=1.2,col=1:2)
mds<-cmdscale(dist(t(Mval)))
rownames(mds)<-sample.annotation2$Sample_Name
g1<-rownames(mds[which(mds[,1]>0),])
g1<-unique(sapply(strsplit(g1,"-"),"[",1))
g0<-rownames(mds[which(mds[,1]<0),])
g0<-unique(sapply(strsplit(g0,"-"),"[",1))
# change header
Mval2<-Mval
colnames(Mval2)<-sample.annotation2$Sample_Nameextract probes on chrX
cgx<-rownames(Locations[which(Locations$chr %in% "chrX"),])
cgy<-rownames(Locations[which(Locations$chr %in% "chrY"),])
mvx<-Mval2[which(rownames(Mval2) %in% cgx),]
mvy<-Mval2[which(rownames(Mval2) %in% cgy),]
plot(colMeans(mvy),colMeans(mvx))
female <- names(which(colMeans(mvy) < -0))
female <- unique(sapply(strsplit(female,"-"),"[",1))
male <- names(which(colMeans(mvy) > 0))
male <- unique(sapply(strsplit(male,"-"),"[",1))
intersect(female,male)## [1] "RG003" "RG012" "RG019" "RG066" "RG071" "RG820"# RG227 looks like a problem
sample.annotation2$malemeth <- as.numeric(colMeans(mvy) > 0)
excl <- c("RG210","RG227")
sample.annotation3 <- sample.annotation2[which(!sample.annotation2$RG_number %in% excl),]
Mval3 <- Mval2[,grep("RG227",colnames(Mval2),invert=TRUE)]
Mval3 <- Mval3[,grep("RG210",colnames(Mval3),invert=TRUE)]
write.table(sample.annotation3,file="samplesheet_mod.tsv",quote=FALSE,sep="\t")
save.image("missmethyl_epic_850k_analysis.Rdata")Differential analysis including sex chromosomes
T0 vs POD analysis
Here we compare pre-op samples to post-op samples.
name="t0vPod_in"
timepoint <- factor(as.numeric(grepl("POD",sample.annotation3$Sample_Name)))
age <- sample.annotation3$age
malemeth <- factor(sample.annotation3$malemeth)
design <- model.matrix(~age + malemeth + timepoint)
fit.reduced <- lmFit(Mval3,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))## (Intercept) age malemeth1 timepoint1
## Down 283086 57161 7747 42361
## NotSig 65917 785763 835269 800016
## Up 508899 14978 14886 15525dm <- topTable(fit.reduced,coef=4, number = Inf)
dma <- merge(myann,dm,by=0)
dma <- dma[order(dma$P.Value),]
head(dma,10)## Row.names UCSC_RefGene_Name
## 219838 cg06528771
## 439011 cg13373048
## 787677 cg25310867 GAB2;GAB2
## 229279 cg06804705 NCRNA00114;NCRNA00114;NCRNA00114
## 844305 cg27307975 LIPC
## 501755 cg15237707 KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15
## 651916 cg20361768 LINC00880
## 446034 cg13574263
## 500621 cg15200345
## 558377 cg17051905 CLN8
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 219838 -1.1165239 0.8441231 -8.352955 6.388786e-13
## 439011 -1.3120256 0.7408465 -7.941312 4.628994e-12
## 787677 -1.0161173 0.9628818 -7.736638 1.232016e-11
## 229279 -0.4807061 1.0314635 -7.530582 3.284968e-11
## 844305 -0.6051741 1.7051766 -7.477450 4.226437e-11
## 501755 -0.6731052 2.3662241 -7.181961 1.703888e-10
## 651916 -0.6671936 0.7881934 -7.174268 1.766544e-10
## 446034 -0.6050069 1.0790662 -7.074552 2.818693e-10
## 500621 -0.5595581 1.3379966 -7.054061 3.102098e-10
## 558377 -0.6125495 2.8867391 -6.724255 1.434478e-09
## adj.P.Val B
## 219838 5.480953e-07 18.35129
## 439011 1.985612e-06 16.55460
## 787677 3.523165e-06 15.66500
## 229279 7.045452e-06 14.77287
## 844305 7.251738e-06 14.54350
## 501755 2.165031e-05 13.27368
## 651916 2.165031e-05 13.24077
## 446034 2.956996e-05 12.81488
## 500621 2.956996e-05 12.72754
## 558377 1.212792e-04 11.33089dm_up <- rownames(subset(dm,adj.P.Val<0.05 & logFC>0))
dm_dn <- rownames(subset(dm,adj.P.Val<0.05 & logFC<0))
length(dm_up)## [1] 15525length(dm_dn)## [1] 42361confects <- limma_confects(fit.reduced, coef=4, fdr=0.05)
make_dm_plots(dm = dm ,name=name , mx=Mval3, beta=beta, groups=timepoint, confects=confects)
t0vPod_in <- list("dma"=dma, "dm_up"=dm_up, "dm_dn"=dm_dn, "confects"=confects)Low vs High CRP samples post-op
Here we compare the post-op samples of patients with low CRP to high CRP.
name="t1_crp_in"
sample.annotation3$timepoint <- factor(as.numeric(grepl("POD",sample.annotation3$Sample_Name)))
sample.annotation4 <- sample.annotation3[sample.annotation3$timepoint==1,]
Mval4 <- Mval3[,which(colnames(Mval3) %in% sample.annotation4$Sample_Name)]
crpgrp <- factor(sample.annotation4$CrpGroup,levels=c(0,1))
age <- sample.annotation4$age
malemeth <- factor(sample.annotation4$malemeth)
design <- model.matrix(~age + malemeth + crpgrp)
fit.reduced <- lmFit(Mval4,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))## (Intercept) age malemeth1 crpgrp1
## Down 272215 4 5396 474
## NotSig 99753 857892 849036 856736
## Up 485934 6 3470 692dm <- topTable(fit.reduced,coef=4, number = Inf)
dma <- merge(myann,dm,by=0)
dma <- dma[order(dma$P.Value),]
head(dma,10)## Row.names UCSC_RefGene_Name
## 501755 cg15237707 KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15
## 464536 cg14107762 TNFSF8;TNFSF8
## 193227 cg05705212 CXCR6;FYCO1
## 707434 cg22431883 CEP41;CEP41;CEP41;CEP41
## 705480 cg22355342 EXOC4
## 189285 cg05584657
## 81324 cg02364564
## 723115 cg23005797 C2orf48
## 398325 cg12048331 WDR82
## 595137 cg18289508
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 501755 -0.8588788 2.02131916 -6.716355 2.343642e-08
## 464536 0.9404818 0.06612697 6.271743 1.098431e-07
## 193227 0.9240981 1.71436691 6.102953 1.973110e-07
## 707434 0.7499884 1.40422320 6.005958 2.761474e-07
## 705480 0.7287797 -0.03419336 5.990044 2.917937e-07
## 189285 0.8060025 0.67974802 5.980152 3.019620e-07
## 81324 0.6561633 1.63938996 5.961280 3.223491e-07
## 723115 0.7809750 2.48776257 5.933411 3.549893e-07
## 398325 1.1017611 0.57664509 5.900080 3.983735e-07
## 595137 0.6270399 1.45079160 5.899955 3.985456e-07
## adj.P.Val B
## 501755 0.02010615 8.530667
## 464536 0.02831762 7.194074
## 193227 0.02831762 6.685308
## 707434 0.02831762 6.392885
## 705480 0.02831762 6.344914
## 189285 0.02831762 6.315094
## 81324 0.02831762 6.258209
## 723115 0.02831762 6.174209
## 398325 0.02831762 6.073763
## 595137 0.02831762 6.073387dm_up <- rownames(subset(dm,adj.P.Val<0.05 & logFC>0))
dm_dn <- rownames(subset(dm,adj.P.Val<0.05 & logFC<0))
length(dm_up)## [1] 692length(dm_dn)## [1] 474confects <- limma_confects(fit.reduced, coef=4, fdr=0.05)
make_dm_plots(dm = dm ,name=name , mx=Mval4, beta=beta, groups=crpgrp, confects=confects)
t0_crp_in <- list("dma"=dma, "dm_up"=dm_up, "dm_dn"=dm_dn, "confects"=confects)Low vs High CRP samples baseline
Here we compare the baseline samples of patients with low CRP to high CRP.
name="t0_crp_in"
sample.annotation4 <- sample.annotation3[sample.annotation3$timepoint==0,]
Mval4 <- Mval3[,which(colnames(Mval3) %in% sample.annotation4$Sample_Name)]
crpgrp <- factor(sample.annotation4$CrpGroup,levels=c(0,1))
age <- sample.annotation4$age
malemeth <- factor(sample.annotation4$malemeth)
design <- model.matrix(~age + malemeth + crpgrp)
fit.reduced <- lmFit(Mval4,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))## (Intercept) age malemeth1 crpgrp1
## Down 273416 8 6372 0
## NotSig 92300 857883 843385 857902
## Up 492186 11 8145 0dm <- topTable(fit.reduced,coef=4, number = Inf)
dma <- merge(myann,dm,by=0)
dma <- dma[order(dma$P.Value),]
head(dma,10)## Row.names UCSC_RefGene_Name Regulatory_Feature_Group logFC
## 319051 cg09538685 PPAP2A;PPAP2A;PPAP2A -0.2991391
## 68537 cg01987065 ESYT2 -0.2772638
## 478405 cg14507403 LMTK2 Unclassified 1.4098408
## 839138 cg27123414 0.3723418
## 30698 cg00873371 -0.2957734
## 658377 cg20604021 PTTG1IP 0.3478706
## 124384 cg03634967 SLC29A1;SLC29A1;SLC29A1 Unclassified 0.4677247
## 623821 cg19345678 -0.2817084
## 221316 cg06572904 -0.3885284
## 103641 cg03016486 MIR8081 -0.3364394
## AveExpr t P.Value adj.P.Val B
## 319051 2.585368 -5.081047 6.475872e-06 0.9999005 0.1615498
## 68537 3.193181 -4.724893 2.140865e-05 0.9999005 -0.3589884
## 478405 -1.157954 4.712643 2.229836e-05 0.9999005 -0.3769194
## 839138 1.918915 4.627384 2.957863e-05 0.9999005 -0.5017033
## 30698 1.999202 -4.603086 3.204988e-05 0.9999005 -0.5372562
## 658377 2.054446 4.575167 3.513993e-05 0.9999005 -0.5781006
## 124384 -4.387651 4.547590 3.847791e-05 0.9999005 -0.6184329
## 623821 2.152636 -4.540637 3.936745e-05 0.9999005 -0.6286009
## 221316 3.634324 -4.468882 4.980652e-05 0.9999005 -0.7334729
## 103641 2.078431 -4.448418 5.325028e-05 0.9999005 -0.7633589dm_up <- rownames(subset(dm,adj.P.Val<0.05 & logFC>0))
dm_dn <- rownames(subset(dm,adj.P.Val<0.05 & logFC<0))
length(dm_up)## [1] 0length(dm_dn)## [1] 0confects <- limma_confects(fit.reduced, coef=4, fdr=0.05)
make_dm_plots(dm = dm ,name=name , mx=Mval4, beta=beta, groups=crpgrp, confects=confects)
t0_crp_in <- list("dma"=dma, "dm_up"=dm_up, "dm_dn"=dm_dn, "confects"=confects)Differential analysis excluding sex chromosomes
T0 vs POD analysis
Here we compare pre-op samples to post-op samples, excluding sex chromosomes.
name="t0vPod_ex"
cgx<-rownames(Locations[which(Locations$chr %in% "chrX"),])
cgy<-rownames(Locations[which(Locations$chr %in% "chrY"),])
cgxy<-union(cgx,cgy)
Mval5 <- Mval3[which(!rownames(Mval3) %in% cgxy),]
timepoint <- factor(as.numeric(grepl("POD",sample.annotation3$Sample_Name)))
age <- sample.annotation3$age
malemeth <- factor(sample.annotation3$malemeth)
design <- model.matrix(~age + malemeth + timepoint)
fit.reduced <- lmFit(Mval5,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))## (Intercept) age malemeth1 timepoint1
## Down 279128 48357 477 42193
## NotSig 60253 779393 835128 781118
## Up 499569 11200 3345 15639dm <- topTable(fit.reduced,coef=4, number = Inf)
dma <- merge(myann,dm,by=0)
dma <- dma[order(dma$P.Value),]
head(dma,10)## Row.names UCSC_RefGene_Name
## 214611 cg06528771
## 429028 cg13373048
## 770258 cg25310867 GAB2;GAB2
## 223845 cg06804705 NCRNA00114;NCRNA00114;NCRNA00114
## 825671 cg27307975 LIPC
## 490426 cg15237707 KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15
## 637274 cg20361768 LINC00880
## 435911 cg13574263
## 489326 cg15200345
## 545758 cg17051905 CLN8
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 214611 -1.1165239 0.8441231 -8.361912 6.022724e-13
## 429028 -1.3120256 0.7408465 -7.950301 4.373389e-12
## 770258 -1.0161173 0.9628818 -7.744886 1.169481e-11
## 223845 -0.4807061 1.0314635 -7.534256 3.190418e-11
## 825671 -0.6051741 1.7051766 -7.483237 4.064925e-11
## 490426 -0.6731052 2.3662241 -7.188379 1.636421e-10
## 637274 -0.6671936 0.7881934 -7.180640 1.697033e-10
## 435911 -0.6050069 1.0790662 -7.080375 2.716167e-10
## 489326 -0.5595581 1.3379966 -7.059384 2.996574e-10
## 545758 -0.6125495 2.8867391 -6.730126 1.384436e-09
## adj.P.Val B
## 214611 5.052764e-07 18.41024
## 429028 1.834527e-06 16.61087
## 770258 3.270453e-06 15.71666
## 223845 6.691503e-06 14.80337
## 825671 6.820537e-06 14.58280
## 490426 2.033894e-05 13.31383
## 637274 2.033894e-05 13.28067
## 435911 2.793307e-05 12.85182
## 489326 2.793307e-05 12.76222
## 545758 1.161472e-04 11.36590dm_up <- rownames(subset(dm,adj.P.Val<0.05 & logFC>0))
dm_dn <- rownames(subset(dm,adj.P.Val<0.05 & logFC<0))
length(dm_up)## [1] 15639length(dm_dn)## [1] 42193confects <- limma_confects(fit.reduced, coef=4, fdr=0.05)
make_dm_plots(dm = dm ,name=name , mx=Mval5, beta=beta, groups=timepoint, confects=confects)
t0vPod_ex <- list("dma"=dma, "dm_up"=dm_up, "dm_dn"=dm_dn, "confects"=confects)Low vs High CRP samples post-op
Here we compare the post-op samples of patients with low CRP to high CRP excluding sex chromosomes.
name="t1_crp_ex"
sample.annotation4 <- sample.annotation3[sample.annotation3$timepoint==1,]
Mval6 <- Mval3[,which(colnames(Mval3) %in% sample.annotation4$Sample_Name)]
Mval6 <- Mval6[which(!rownames(Mval6) %in% cgxy),]
crpgrp <- factor(sample.annotation4$CrpGroup,levels=c(0,1))
age <- sample.annotation4$age
malemeth <- factor(sample.annotation4$malemeth)
design <- model.matrix(~age + malemeth + crpgrp)
fit.reduced <- lmFit(Mval6,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))## (Intercept) age malemeth1 crpgrp1
## Down 269066 2 15 584
## NotSig 91083 838942 838899 837596
## Up 478801 6 36 770dm <- topTable(fit.reduced,coef=4, number = Inf)
dma <- merge(myann,dm,by=0)
dma <- dma[order(dma$P.Value),]
head(dma,10)## Row.names UCSC_RefGene_Name
## 490426 cg15237707 KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15
## 453993 cg14107762 TNFSF8;TNFSF8
## 188665 cg05705212 CXCR6;FYCO1
## 691656 cg22431883 CEP41;CEP41;CEP41;CEP41
## 689746 cg22355342 EXOC4
## 184810 cg05584657
## 79352 cg02364564
## 707001 cg23005797 C2orf48
## 389265 cg12048331 WDR82
## 581675 cg18289508
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 490426 -0.8588788 2.02131916 -6.725911 2.217949e-08
## 453993 0.9404818 0.06612697 6.282334 1.039466e-07
## 188665 0.9240981 1.71436691 6.113342 1.870976e-07
## 691656 0.7499884 1.40422320 6.014218 2.639948e-07
## 689746 0.7287797 -0.03419336 5.997956 2.793246e-07
## 184810 0.8060025 0.67974802 5.989242 2.878999e-07
## 79352 0.6561633 1.63938996 5.967731 3.102096e-07
## 707001 0.7809750 2.48776257 5.942178 3.389595e-07
## 389265 1.1017611 0.57664509 5.911531 3.769579e-07
## 581675 0.6270399 1.45079160 5.905777 3.845511e-07
## adj.P.Val B
## 490426 0.01860749 8.584643
## 453993 0.02658155 7.246987
## 188665 0.02658155 6.736060
## 691656 0.02658155 6.436317
## 689746 0.02658155 6.387147
## 184810 0.02658155 6.360802
## 79352 0.02658155 6.295768
## 707001 0.02658155 6.218518
## 389265 0.02658155 6.125885
## 581675 0.02658155 6.108495dm_up <- rownames(subset(dm,adj.P.Val<0.05 & logFC>0))
dm_dn <- rownames(subset(dm,adj.P.Val<0.05 & logFC<0))
length(dm_up)## [1] 770length(dm_dn)## [1] 584confects <- limma_confects(fit.reduced, coef=4, fdr=0.05)
make_dm_plots(dm = dm ,name=name , mx=Mval6, beta=beta, groups=crpgrp, confects=confects)
t0_crp_ex <- list("dma"=dma, "dm_up"=dm_up, "dm_dn"=dm_dn, "confects"=confects)Low vs High CRP samples
Here we compare the baseliine samples of patients with low CRP to high CRP excluding sex chromosomes.
name="t0_crp_ex"
sample.annotation4 <- sample.annotation3[sample.annotation3$timepoint==0,]
Mval6 <- Mval3[,which(colnames(Mval3) %in% sample.annotation4$Sample_Name)]
Mval6 <- Mval6[which(!rownames(Mval6) %in% cgxy),]
crpgrp <- factor(sample.annotation4$CrpGroup,levels=c(0,1))
age <- sample.annotation4$age
malemeth <- factor(sample.annotation4$malemeth)
design <- model.matrix(~age + malemeth + crpgrp)
fit.reduced <- lmFit(Mval6,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))## (Intercept) age malemeth1 crpgrp1
## Down 270096 8 156 0
## NotSig 84117 838931 838126 838950
## Up 484737 11 668 0dm <- topTable(fit.reduced,coef=4, number = Inf)
dma <- merge(myann,dm,by=0)
dma <- dma[order(dma$P.Value),]
head(dma,10)## Row.names UCSC_RefGene_Name Regulatory_Feature_Group logFC
## 311708 cg09538685 PPAP2A;PPAP2A;PPAP2A -0.2991391
## 467564 cg14507403 LMTK2 Unclassified 1.4098408
## 66884 cg01987065 ESYT2 -0.2772638
## 820601 cg27123414 0.3723418
## 29973 cg00873371 -0.2957734
## 643605 cg20604021 PTTG1IP 0.3478706
## 121467 cg03634967 SLC29A1;SLC29A1;SLC29A1 Unclassified 0.4677247
## 609757 cg19345678 -0.2817084
## 216056 cg06572904 -0.3885284
## 101150 cg03016486 MIR8081 -0.3364394
## AveExpr t P.Value adj.P.Val B
## 311708 2.585368 -5.070976 6.633986e-06 0.9999842 0.1773347
## 467564 -1.157954 4.722942 2.137394e-05 0.9999842 -0.3355410
## 66884 3.193181 -4.715395 2.191779e-05 0.9999842 -0.3466769
## 820601 1.918915 4.627503 2.934207e-05 0.9999842 -0.4763340
## 29973 1.999202 -4.597148 3.244015e-05 0.9999842 -0.5210991
## 643605 2.054446 4.574012 3.501462e-05 0.9999842 -0.5552103
## 121467 -4.387651 4.551807 3.767325e-05 0.9999842 -0.5879413
## 609757 2.152636 -4.533585 4.000188e-05 0.9999842 -0.6147953
## 216056 3.634324 -4.470485 4.920603e-05 0.9999842 -0.7077284
## 101150 2.078431 -4.447172 5.310609e-05 0.9999842 -0.7420386dm_up <- rownames(subset(dm,adj.P.Val<0.05 & logFC>0))
dm_dn <- rownames(subset(dm,adj.P.Val<0.05 & logFC<0))
length(dm_up)## [1] 0length(dm_dn)## [1] 0confects <- limma_confects(fit.reduced, coef=4, fdr=0.05)
make_dm_plots(dm = dm ,name=name , mx=Mval6, beta=beta, groups=crpgrp, confects=confects)
t0_crp_ex <- list("dma"=dma, "dm_up"=dm_up, "dm_dn"=dm_dn, "confects"=confects)Normalise for blood composition and repeat the analysis
basedir = "analysis/idat"
rgSet <- read.metharray.exp(basedir, targets = sample.annotation2)## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nuls
## Warning in readChar(con, nchars = n): truncating string with embedded nulscells <- estimateCellCounts2(rgSet, referencePlatform= "IlluminaHumanMethylationEPIC", returnAll = TRUE)## snapshotDate(): 2020-04-27## see ?FlowSorted.Blood.EPIC and browseVignettes('FlowSorted.Blood.EPIC') for documentation## loading from cache## [estimateCellCounts2] Combining user data with reference (flow sorted) data.## Warning in asMethod(object): NAs introduced by coercion## [estimateCellCounts2] Processing user and reference data together.## [estimateCellCounts2] Picking probes for composition estimation.## [estimateCellCounts2] Estimating composition.mset <- cells$normalizedData
densityByProbeType(mSet[,1], main = "Raw")
densityByProbeType(mset[,1], main = "estimateCellCounts2")
mset_reduced <- mset[which(rownames(mset) %in% names(keep[keep==TRUE])),]
meth <- getMeth(mset_reduced)
unmeth <- getUnmeth(mset_reduced)
Mval <- log2((meth + 100)/(unmeth + 100))
beta <- getBeta(mset_reduced)
dim(Mval)## [1] 857902 96plot mds
plotMDS(Mval, labels=sample.annotation2$Sample_Name, col=as.integer(factor(sample.annotation2$CrpGroup)))
legend("topleft",legend=c("CRP high","CRP norm"),pch=16,cex=1.2,col=1:2)
mds<-cmdscale(dist(t(Mval)))
rownames(mds)<-sample.annotation2$Sample_Name
g1<-rownames(mds[which(mds[,1]>0),])
g1<-unique(sapply(strsplit(g1,"-"),"[",1))
g0<-rownames(mds[which(mds[,1]<0),])
g0<-unique(sapply(strsplit(g0,"-"),"[",1))
# change header
Mval2<-Mval
colnames(Mval2)<-sample.annotation2$Sample_NameDifferential analysis including sex chromosomes
T0 vs POD analysis
Here we compare pre-op samples to post-op samples.
name="t0vPod_in_blood"
timepoint <- factor(as.numeric(grepl("POD",sample.annotation3$Sample_Name)))
age <- sample.annotation3$age
malemeth <- factor(sample.annotation3$malemeth)
design <- model.matrix(~age + malemeth + timepoint)
fit.reduced <- lmFit(Mval3,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))## (Intercept) age malemeth1 timepoint1
## Down 283086 57161 7747 42361
## NotSig 65917 785763 835269 800016
## Up 508899 14978 14886 15525dm <- topTable(fit.reduced,coef=4, number = Inf)
dma <- merge(myann,dm,by=0)
dma <- dma[order(dma$P.Value),]
head(dma,10)## Row.names UCSC_RefGene_Name
## 219838 cg06528771
## 439011 cg13373048
## 787677 cg25310867 GAB2;GAB2
## 229279 cg06804705 NCRNA00114;NCRNA00114;NCRNA00114
## 844305 cg27307975 LIPC
## 501755 cg15237707 KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15
## 651916 cg20361768 LINC00880
## 446034 cg13574263
## 500621 cg15200345
## 558377 cg17051905 CLN8
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 219838 -1.1165239 0.8441231 -8.352955 6.388786e-13
## 439011 -1.3120256 0.7408465 -7.941312 4.628994e-12
## 787677 -1.0161173 0.9628818 -7.736638 1.232016e-11
## 229279 -0.4807061 1.0314635 -7.530582 3.284968e-11
## 844305 -0.6051741 1.7051766 -7.477450 4.226437e-11
## 501755 -0.6731052 2.3662241 -7.181961 1.703888e-10
## 651916 -0.6671936 0.7881934 -7.174268 1.766544e-10
## 446034 -0.6050069 1.0790662 -7.074552 2.818693e-10
## 500621 -0.5595581 1.3379966 -7.054061 3.102098e-10
## 558377 -0.6125495 2.8867391 -6.724255 1.434478e-09
## adj.P.Val B
## 219838 5.480953e-07 18.35129
## 439011 1.985612e-06 16.55460
## 787677 3.523165e-06 15.66500
## 229279 7.045452e-06 14.77287
## 844305 7.251738e-06 14.54350
## 501755 2.165031e-05 13.27368
## 651916 2.165031e-05 13.24077
## 446034 2.956996e-05 12.81488
## 500621 2.956996e-05 12.72754
## 558377 1.212792e-04 11.33089dm_up <- rownames(subset(dm,adj.P.Val<0.05 & logFC>0))
dm_dn <- rownames(subset(dm,adj.P.Val<0.05 & logFC<0))
length(dm_up)## [1] 15525length(dm_dn)## [1] 42361confects <- limma_confects(fit.reduced, coef=4, fdr=0.05)
make_dm_plots(dm = dm ,name=name , mx=Mval3, beta=beta, groups=timepoint, confects=confects)
t0vPod_in_blood <- list("dma"=dma, "dm_up"=dm_up, "dm_dn"=dm_dn, "confects"=confects)Low vs High CRP samples post-op
Here we compare the post-op samples of patients with low CRP to high CRP.
name="t1_crp_in_blood"
sample.annotation4 <- sample.annotation3[sample.annotation3$timepoint==1,]
Mval4 <- Mval3[,which(colnames(Mval3) %in% sample.annotation4$Sample_Name)]
crpgrp <- factor(sample.annotation4$CrpGroup,levels=c(0,1))
age <- sample.annotation4$age
malemeth <- factor(sample.annotation4$malemeth)
design <- model.matrix(~age + malemeth + crpgrp)
fit.reduced <- lmFit(Mval4,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))## (Intercept) age malemeth1 crpgrp1
## Down 272215 4 5396 474
## NotSig 99753 857892 849036 856736
## Up 485934 6 3470 692dm <- topTable(fit.reduced,coef=4, number = Inf)
dma <- merge(myann,dm,by=0)
dma <- dma[order(dma$P.Value),]
head(dma,10)## Row.names UCSC_RefGene_Name
## 501755 cg15237707 KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15
## 464536 cg14107762 TNFSF8;TNFSF8
## 193227 cg05705212 CXCR6;FYCO1
## 707434 cg22431883 CEP41;CEP41;CEP41;CEP41
## 705480 cg22355342 EXOC4
## 189285 cg05584657
## 81324 cg02364564
## 723115 cg23005797 C2orf48
## 398325 cg12048331 WDR82
## 595137 cg18289508
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 501755 -0.8588788 2.02131916 -6.716355 2.343642e-08
## 464536 0.9404818 0.06612697 6.271743 1.098431e-07
## 193227 0.9240981 1.71436691 6.102953 1.973110e-07
## 707434 0.7499884 1.40422320 6.005958 2.761474e-07
## 705480 0.7287797 -0.03419336 5.990044 2.917937e-07
## 189285 0.8060025 0.67974802 5.980152 3.019620e-07
## 81324 0.6561633 1.63938996 5.961280 3.223491e-07
## 723115 0.7809750 2.48776257 5.933411 3.549893e-07
## 398325 1.1017611 0.57664509 5.900080 3.983735e-07
## 595137 0.6270399 1.45079160 5.899955 3.985456e-07
## adj.P.Val B
## 501755 0.02010615 8.530667
## 464536 0.02831762 7.194074
## 193227 0.02831762 6.685308
## 707434 0.02831762 6.392885
## 705480 0.02831762 6.344914
## 189285 0.02831762 6.315094
## 81324 0.02831762 6.258209
## 723115 0.02831762 6.174209
## 398325 0.02831762 6.073763
## 595137 0.02831762 6.073387dm_up <- rownames(subset(dm,adj.P.Val<0.05 & logFC>0))
dm_dn <- rownames(subset(dm,adj.P.Val<0.05 & logFC<0))
length(dm_up)## [1] 692length(dm_dn)## [1] 474confects <- limma_confects(fit.reduced, coef=4, fdr=0.05)
make_dm_plots(dm = dm ,name=name , mx=Mval4, beta=beta, groups=crpgrp, confects=confects)
t0_crp_in_blood <- list("dma"=dma, "dm_up"=dm_up, "dm_dn"=dm_dn, "confects"=confects)Low vs High CRP samples baseline
Here we compare the baseline samples of patients with low CRP to high CRP.
name="t0_crp_in_blood"
sample.annotation4 <- sample.annotation3[sample.annotation3$timepoint==0,]
Mval4 <- Mval3[,which(colnames(Mval3) %in% sample.annotation4$Sample_Name)]
crpgrp <- factor(sample.annotation4$CrpGroup,levels=c(0,1))
age <- sample.annotation4$age
malemeth <- factor(sample.annotation4$malemeth)
design <- model.matrix(~age + malemeth + crpgrp)
fit.reduced <- lmFit(Mval4,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))## (Intercept) age malemeth1 crpgrp1
## Down 273416 8 6372 0
## NotSig 92300 857883 843385 857902
## Up 492186 11 8145 0dm <- topTable(fit.reduced,coef=4, number = Inf)
dma <- merge(myann,dm,by=0)
dma <- dma[order(dma$P.Value),]
head(dma,10)## Row.names UCSC_RefGene_Name Regulatory_Feature_Group logFC
## 319051 cg09538685 PPAP2A;PPAP2A;PPAP2A -0.2991391
## 68537 cg01987065 ESYT2 -0.2772638
## 478405 cg14507403 LMTK2 Unclassified 1.4098408
## 839138 cg27123414 0.3723418
## 30698 cg00873371 -0.2957734
## 658377 cg20604021 PTTG1IP 0.3478706
## 124384 cg03634967 SLC29A1;SLC29A1;SLC29A1 Unclassified 0.4677247
## 623821 cg19345678 -0.2817084
## 221316 cg06572904 -0.3885284
## 103641 cg03016486 MIR8081 -0.3364394
## AveExpr t P.Value adj.P.Val B
## 319051 2.585368 -5.081047 6.475872e-06 0.9999005 0.1615498
## 68537 3.193181 -4.724893 2.140865e-05 0.9999005 -0.3589884
## 478405 -1.157954 4.712643 2.229836e-05 0.9999005 -0.3769194
## 839138 1.918915 4.627384 2.957863e-05 0.9999005 -0.5017033
## 30698 1.999202 -4.603086 3.204988e-05 0.9999005 -0.5372562
## 658377 2.054446 4.575167 3.513993e-05 0.9999005 -0.5781006
## 124384 -4.387651 4.547590 3.847791e-05 0.9999005 -0.6184329
## 623821 2.152636 -4.540637 3.936745e-05 0.9999005 -0.6286009
## 221316 3.634324 -4.468882 4.980652e-05 0.9999005 -0.7334729
## 103641 2.078431 -4.448418 5.325028e-05 0.9999005 -0.7633589dm_up <- rownames(subset(dm,adj.P.Val<0.05 & logFC>0))
dm_dn <- rownames(subset(dm,adj.P.Val<0.05 & logFC<0))
length(dm_up)## [1] 0length(dm_dn)## [1] 0confects <- limma_confects(fit.reduced, coef=4, fdr=0.05)
make_dm_plots(dm = dm ,name=name , mx=Mval4, beta=beta, groups=crpgrp, confects=confects)
t0_crp_in_blood <- list("dma"=dma, "dm_up"=dm_up, "dm_dn"=dm_dn, "confects"=confects)Differential analysis excluding sex chromosomes
T0 vs POD analysis
Here we compare pre-op samples to post-op samples, excluding X and Y chromosomes.
name="t0vPod_ex_blood"
cgx<-rownames(Locations[which(Locations$chr %in% "chrX"),])
cgy<-rownames(Locations[which(Locations$chr %in% "chrY"),])
cgxy<-union(cgx,cgy)
Mval5 <- Mval3[which(!rownames(Mval3) %in% cgxy),]
timepoint <- factor(as.numeric(grepl("POD",sample.annotation3$Sample_Name)))
age <- sample.annotation3$age
malemeth <- factor(sample.annotation3$malemeth)
design <- model.matrix(~age + malemeth + timepoint)
fit.reduced <- lmFit(Mval5,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))## (Intercept) age malemeth1 timepoint1
## Down 279128 48357 477 42193
## NotSig 60253 779393 835128 781118
## Up 499569 11200 3345 15639dm <- topTable(fit.reduced,coef=4, number = Inf)
dma <- merge(myann,dm,by=0)
dma <- dma[order(dma$P.Value),]
head(dma,10)## Row.names UCSC_RefGene_Name
## 214611 cg06528771
## 429028 cg13373048
## 770258 cg25310867 GAB2;GAB2
## 223845 cg06804705 NCRNA00114;NCRNA00114;NCRNA00114
## 825671 cg27307975 LIPC
## 490426 cg15237707 KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15
## 637274 cg20361768 LINC00880
## 435911 cg13574263
## 489326 cg15200345
## 545758 cg17051905 CLN8
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 214611 -1.1165239 0.8441231 -8.361912 6.022724e-13
## 429028 -1.3120256 0.7408465 -7.950301 4.373389e-12
## 770258 -1.0161173 0.9628818 -7.744886 1.169481e-11
## 223845 -0.4807061 1.0314635 -7.534256 3.190418e-11
## 825671 -0.6051741 1.7051766 -7.483237 4.064925e-11
## 490426 -0.6731052 2.3662241 -7.188379 1.636421e-10
## 637274 -0.6671936 0.7881934 -7.180640 1.697033e-10
## 435911 -0.6050069 1.0790662 -7.080375 2.716167e-10
## 489326 -0.5595581 1.3379966 -7.059384 2.996574e-10
## 545758 -0.6125495 2.8867391 -6.730126 1.384436e-09
## adj.P.Val B
## 214611 5.052764e-07 18.41024
## 429028 1.834527e-06 16.61087
## 770258 3.270453e-06 15.71666
## 223845 6.691503e-06 14.80337
## 825671 6.820537e-06 14.58280
## 490426 2.033894e-05 13.31383
## 637274 2.033894e-05 13.28067
## 435911 2.793307e-05 12.85182
## 489326 2.793307e-05 12.76222
## 545758 1.161472e-04 11.36590dm_up <- rownames(subset(dm,adj.P.Val<0.05 & logFC>0))
dm_dn <- rownames(subset(dm,adj.P.Val<0.05 & logFC<0))
length(dm_up)## [1] 15639length(dm_dn)## [1] 42193confects <- limma_confects(fit.reduced, coef=4, fdr=0.05)
make_dm_plots(dm = dm ,name=name , mx=Mval5, beta=beta, groups=timepoint, confects=confects)
t0vPod_ex_blood <- list("dma"=dma, "dm_up"=dm_up, "dm_dn"=dm_dn, "confects"=confects)
# pathway
dmap <- dma[grep("Promoter_Associated",dma$Regulatory_Feature_Group),]
dmap[which(dmap$UCSC_RefGene_Name==""),2] <- "NA"
dmap$genename <- sapply(strsplit(dmap$UCSC_RefGene_Name,";"),"[[",1)
dmap2 <- dmap[,c("genename","t")]
rank <- aggregate(. ~ genename,dmap2,mean)
rownames(rank) <- rank$genename
rank$genename=NULL
write.table(rank,file="t0vPod_meth.tsv",sep="\t",quote=FALSE)
capture.output(
res <- mitch_calc(x = rank,genesets = genesets, priority = "significance",resrows=20)
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)## Note: When prioritising by significance (ie: small
## p-values), large effect sizes might be missed.head(res$enrichment_result,20)## set
## 510 Metabolism of RNA
## 855 SRP-dependent cotranslational protein targeting to membrane
## 417 Influenza Infection
## 418 Influenza Viral RNA Transcription and Replication
## 306 Formation of a pool of free 40S subunits
## 658 Peptide chain elongation
## 274 Eukaryotic Translation Elongation
## 341 GTP hydrolysis and joining of the 60S ribosomal subunit
## 473 L13a-mediated translational silencing of Ceruloplasmin expression
## 276 Eukaryotic Translation Termination
## 132 Cellular responses to external stimuli
## 110 Cap-dependent Translation Initiation
## 275 Eukaryotic Translation Initiation
## 133 Cellular responses to stress
## 522 Metabolism of proteins
## 501 Major pathway of rRNA processing in the nucleolus and cytosol
## 1096 Viral mRNA Translation
## 140 Chromatin modifying enzymes
## 141 Chromatin organization
## 608 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)
## setSize pANOVA s.dist p.adjustANOVA
## 510 542 4.429953e-07 0.12826184 0.0005023567
## 855 85 5.019081e-06 0.28678256 0.0028213095
## 417 129 7.463782e-06 0.22895905 0.0028213095
## 418 111 1.273834e-05 0.24027785 0.0030608593
## 306 76 1.349585e-05 0.28909155 0.0030608593
## 658 67 1.938307e-05 0.30211897 0.0032967624
## 274 70 2.207766e-05 0.29360051 0.0032967624
## 341 86 2.325758e-05 0.26434218 0.0032967624
## 473 85 4.020172e-05 0.25804686 0.0044061647
## 276 70 4.028612e-05 0.28413148 0.0044061647
## 132 400 4.604405e-05 0.11972562 0.0044061647
## 110 93 5.051159e-05 0.24359888 0.0044061647
## 275 93 5.051159e-05 0.24359888 0.0044061647
## 133 395 5.689335e-05 0.11899103 0.0046083614
## 522 1288 6.792172e-05 0.06803682 0.0049488382
## 501 149 7.371317e-05 0.18868564 0.0049488382
## 1096 68 7.418893e-05 0.27819386 0.0049488382
## 140 155 8.821978e-05 0.18303769 0.0052653279
## 141 155 8.821978e-05 0.18303769 0.0052653279
## 608 72 1.184672e-04 0.26268661 0.0067170879unlink("t0vPod_ex_blood.html")
capture.output(
mitch_plots(res,outfile="t0vPod_ex_blood.pdf")
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)Low vs High CRP samples post-op
Here we compare the post-op samples of patients with low CRP to high CRP excluding sex chromosomes.
name="t1_crp_ex_blood"
sample.annotation4 <- sample.annotation3[sample.annotation3$timepoint==1,]
Mval6 <- Mval3[,which(colnames(Mval3) %in% sample.annotation4$Sample_Name)]
Mval6 <- Mval6[which(!rownames(Mval6) %in% cgxy),]
crpgrp <- factor(sample.annotation4$CrpGroup,levels=c(0,1))
age <- sample.annotation4$age
malemeth <- factor(sample.annotation4$malemeth)
design <- model.matrix(~age + malemeth + crpgrp)
fit.reduced <- lmFit(Mval6,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))## (Intercept) age malemeth1 crpgrp1
## Down 269066 2 15 584
## NotSig 91083 838942 838899 837596
## Up 478801 6 36 770dm <- topTable(fit.reduced,coef=4, number = Inf)
dma <- merge(myann,dm,by=0)
dma <- dma[order(dma$P.Value),]
head(dma,10)## Row.names UCSC_RefGene_Name
## 490426 cg15237707 KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15
## 453993 cg14107762 TNFSF8;TNFSF8
## 188665 cg05705212 CXCR6;FYCO1
## 691656 cg22431883 CEP41;CEP41;CEP41;CEP41
## 689746 cg22355342 EXOC4
## 184810 cg05584657
## 79352 cg02364564
## 707001 cg23005797 C2orf48
## 389265 cg12048331 WDR82
## 581675 cg18289508
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 490426 -0.8588788 2.02131916 -6.725911 2.217949e-08
## 453993 0.9404818 0.06612697 6.282334 1.039466e-07
## 188665 0.9240981 1.71436691 6.113342 1.870976e-07
## 691656 0.7499884 1.40422320 6.014218 2.639948e-07
## 689746 0.7287797 -0.03419336 5.997956 2.793246e-07
## 184810 0.8060025 0.67974802 5.989242 2.878999e-07
## 79352 0.6561633 1.63938996 5.967731 3.102096e-07
## 707001 0.7809750 2.48776257 5.942178 3.389595e-07
## 389265 1.1017611 0.57664509 5.911531 3.769579e-07
## 581675 0.6270399 1.45079160 5.905777 3.845511e-07
## adj.P.Val B
## 490426 0.01860749 8.584643
## 453993 0.02658155 7.246987
## 188665 0.02658155 6.736060
## 691656 0.02658155 6.436317
## 689746 0.02658155 6.387147
## 184810 0.02658155 6.360802
## 79352 0.02658155 6.295768
## 707001 0.02658155 6.218518
## 389265 0.02658155 6.125885
## 581675 0.02658155 6.108495dm_up <- rownames(subset(dm,adj.P.Val<0.05 & logFC>0))
dm_dn <- rownames(subset(dm,adj.P.Val<0.05 & logFC<0))
length(dm_up)## [1] 770length(dm_dn)## [1] 584confects <- limma_confects(fit.reduced, coef=4, fdr=0.05)
make_dm_plots(dm = dm ,name=name , mx=Mval6, beta=beta, groups=crpgrp, confects=confects)
t1_crp_ex_blood <- list("dma"=dma, "dm_up"=dm_up, "dm_dn"=dm_dn, "confects"=confects)
# pathway
dmap <- dma[grep("Promoter_Associated",dma$Regulatory_Feature_Group),]
dmap[which(dmap$UCSC_RefGene_Name==""),2] <- "NA"
dmap$genename <- sapply(strsplit(dmap$UCSC_RefGene_Name,";"),"[[",1)
dmap2 <- dmap[,c("genename","t")]
rank <- aggregate(. ~ genename,dmap2,mean)
rownames(rank) <- rank$genename
rank$genename=NULL
write.table(rank,file="t1_crp_meth.tsv",sep="\t",quote=FALSE)
capture.output(
res <- mitch_calc(x = rank,genesets = genesets, priority = "significance",resrows=20)
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)## Note: When prioritising by significance (ie: small
## p-values), large effect sizes might be missed.head(res$enrichment_result,20)## set
## 944 Signaling by WNT
## 983 TCF dependent signaling in response to WNT
## 350 Gene expression (Transcription)
## 510 Metabolism of RNA
## 417 Influenza Infection
## 348 Gene Silencing by RNA
## 418 Influenza Viral RNA Transcription and Replication
## 415 Infectious disease
## 809 Regulation of expression of SLITs and ROBOs
## 133 Cellular responses to stress
## 132 Cellular responses to external stimuli
## 644 PIP3 activates AKT signaling
## 916 Signaling by NOTCH
## 746 RNA Polymerase II Transcription
## 607 Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC)
## 609 Nonsense-Mediated Decay (NMD)
## 383 HIV Infection
## 81 Beta-catenin independent WNT signaling
## 453 Interleukin-4 and Interleukin-13 signaling
## 935 Signaling by ROBO receptors
## setSize pANOVA s.dist p.adjustANOVA
## 944 178 1.553938e-05 0.18840805 0.01180133
## 983 119 3.049523e-05 0.22178944 0.01180133
## 350 1025 3.122046e-05 0.07870903 0.01180133
## 510 542 4.420602e-05 0.10376576 0.01253241
## 417 129 6.801923e-05 0.20358268 0.01542676
## 348 59 1.106539e-04 0.29127638 0.02078066
## 418 111 1.282757e-04 0.21085594 0.02078066
## 415 490 1.810120e-04 0.09980980 0.02397368
## 809 122 1.970751e-04 0.19560490 0.02397368
## 133 395 2.309236e-04 0.10886399 0.02397368
## 132 400 2.325489e-04 0.10815360 0.02397368
## 644 182 3.067759e-04 0.15569665 0.02899033
## 916 133 5.032262e-04 0.17516984 0.03967837
## 746 913 5.417276e-04 0.06891174 0.03967837
## 607 86 5.598359e-04 0.21559442 0.03967837
## 609 86 5.598359e-04 0.21559442 0.03967837
## 383 188 6.514794e-04 0.14473489 0.04345751
## 81 107 7.494407e-04 0.18900610 0.04721476
## 453 48 8.646267e-04 0.27815315 0.05160456
## 935 153 1.080346e-03 0.15358865 0.05477481unlink("t1_crp_ex_blood.html")
capture.output(
mitch_plots(res,outfile="t1_crp_ex_blood.pdf")
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)Low vs High CRP samples
Here we compare the baseline samples of patients with low CRP to high CRP excluding sex chromosomes.
name="t0_crp_ex_blood"
sample.annotation4 <- sample.annotation3[sample.annotation3$timepoint==0,]
Mval6 <- Mval3[,which(colnames(Mval3) %in% sample.annotation4$Sample_Name)]
Mval6 <- Mval6[which(!rownames(Mval6) %in% cgxy),]
crpgrp <- factor(sample.annotation4$CrpGroup,levels=c(0,1))
age <- sample.annotation4$age
malemeth <- factor(sample.annotation4$malemeth)
design <- model.matrix(~age + malemeth + crpgrp)
fit.reduced <- lmFit(Mval6,design)
fit.reduced <- eBayes(fit.reduced)
summary(decideTests(fit.reduced))## (Intercept) age malemeth1 crpgrp1
## Down 270096 8 156 0
## NotSig 84117 838931 838126 838950
## Up 484737 11 668 0dm <- topTable(fit.reduced,coef=4, number = Inf)
dma <- merge(myann,dm,by=0)
dma <- dma[order(dma$P.Value),]
head(dma,10)## Row.names UCSC_RefGene_Name Regulatory_Feature_Group logFC
## 311708 cg09538685 PPAP2A;PPAP2A;PPAP2A -0.2991391
## 467564 cg14507403 LMTK2 Unclassified 1.4098408
## 66884 cg01987065 ESYT2 -0.2772638
## 820601 cg27123414 0.3723418
## 29973 cg00873371 -0.2957734
## 643605 cg20604021 PTTG1IP 0.3478706
## 121467 cg03634967 SLC29A1;SLC29A1;SLC29A1 Unclassified 0.4677247
## 609757 cg19345678 -0.2817084
## 216056 cg06572904 -0.3885284
## 101150 cg03016486 MIR8081 -0.3364394
## AveExpr t P.Value adj.P.Val B
## 311708 2.585368 -5.070976 6.633986e-06 0.9999842 0.1773347
## 467564 -1.157954 4.722942 2.137394e-05 0.9999842 -0.3355410
## 66884 3.193181 -4.715395 2.191779e-05 0.9999842 -0.3466769
## 820601 1.918915 4.627503 2.934207e-05 0.9999842 -0.4763340
## 29973 1.999202 -4.597148 3.244015e-05 0.9999842 -0.5210991
## 643605 2.054446 4.574012 3.501462e-05 0.9999842 -0.5552103
## 121467 -4.387651 4.551807 3.767325e-05 0.9999842 -0.5879413
## 609757 2.152636 -4.533585 4.000188e-05 0.9999842 -0.6147953
## 216056 3.634324 -4.470485 4.920603e-05 0.9999842 -0.7077284
## 101150 2.078431 -4.447172 5.310609e-05 0.9999842 -0.7420386dm_up <- rownames(subset(dm,adj.P.Val<0.05 & logFC>0))
dm_dn <- rownames(subset(dm,adj.P.Val<0.05 & logFC<0))
length(dm_up)## [1] 0length(dm_dn)## [1] 0confects <- limma_confects(fit.reduced, coef=4, fdr=0.05)
make_dm_plots(dm = dm ,name=name , mx=Mval6, beta=beta, groups=crpgrp, confects=confects)
t0_crp_ex_blood <- list("dma"=dma, "dm_up"=dm_up, "dm_dn"=dm_dn, "confects"=confects)
# pathway
dmap <- dma[grep("Promoter_Associated",dma$Regulatory_Feature_Group),]
dmap[which(dmap$UCSC_RefGene_Name==""),2] <- "NA"
dmap$genename <- sapply(strsplit(dmap$UCSC_RefGene_Name,";"),"[[",1)
dmap2 <- dmap[,c("genename","t")]
rank <- aggregate(. ~ genename,dmap2,mean)
rownames(rank) <- rank$genename
rank$genename=NULL
write.table(rank,file="t0_crp_meth.tsv",sep="\t",quote=FALSE)
capture.output(
res <- mitch_calc(x = rank,genesets = genesets, priority = "significance",resrows=20)
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)## Note: When prioritising by significance (ie: small
## p-values), large effect sizes might be missed.head(res$enrichment_result,20)## set
## 866 SUMOylation of intracellular receptors
## 221 Diseases associated with glycosaminoglycan metabolism
## 1000 TP53 Regulates Transcription of Genes Involved in G2 Cell Cycle Arrest
## 396 Heparan sulfate/heparin (HS-GAG) metabolism
## 217 Disassembly of the destruction complex and recruitment of AXIN to the membrane
## 516 Metabolism of lipids
## 983 TCF dependent signaling in response to WNT
## 944 Signaling by WNT
## 367 Glycosaminoglycan metabolism
## 1054 Translation
## 945 Signaling by WNT in cancer
## 541 Mitochondrial tRNA aminoacylation
## 909 Signaling by Hippo
## 698 Processing of Capped Intronless Pre-mRNA
## 1045 Transcriptional activity of SMAD2/SMAD3:SMAD4 heterotrimer
## 851 SHC1 events in ERBB2 signaling
## 497 MET promotes cell motility
## 685 Post-translational modification: synthesis of GPI-anchored proteins
## 991 TNFR1-induced NFkappaB signaling pathway
## 536 Mitochondrial Fatty Acid Beta-Oxidation
## setSize pANOVA s.dist p.adjustANOVA
## 866 19 0.004390609 0.37760179 0.9697381
## 221 16 0.007884791 0.38375242 0.9697381
## 1000 16 0.008280384 -0.38136340 0.9697381
## 396 21 0.008312757 0.33278787 0.9697381
## 217 26 0.009719241 0.29309293 0.9697381
## 516 435 0.014966067 0.06868066 0.9697381
## 983 119 0.015249645 0.12909329 0.9697381
## 944 178 0.020228712 0.10128121 0.9697381
## 367 63 0.023364893 0.16538849 0.9697381
## 1054 234 0.025267619 -0.08532349 0.9697381
## 945 23 0.028366270 0.26417632 0.9697381
## 541 18 0.029666924 -0.29614948 0.9697381
## 909 12 0.031553816 0.35852850 0.9697381
## 698 21 0.033241797 0.26849859 0.9697381
## 1045 37 0.034022893 -0.20153295 0.9697381
## 851 11 0.034323558 0.36857438 0.9697381
## 497 13 0.035171546 0.33748755 0.9697381
## 685 31 0.039475564 0.21382383 0.9697381
## 991 20 0.044740014 -0.25934802 0.9697381
## 536 25 0.045470489 0.23123138 0.9697381unlink("t0_crp_ex_blood.html")
capture.output(
mitch_plots(res,outfile="t0_crp_ex_blood.pdf")
, file = "/dev/null", append = FALSE,
type = c("output", "message"), split = FALSE)