RELIEF: methylation study
Mark Ziemann
2020-09-14
- 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" "RG867"# 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 283900 57501 7937 41135
## NotSig 66029 784148 835091 802310
## Up 507973 16253 14874 14457dm <- 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
## 651916 cg20361768 LINC00880
## 501755 cg15237707 KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15
## 446034 cg13574263
## 500621 cg15200345
## 558377 cg17051905 CLN8
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 219838 -1.1206545 0.8453538 -8.337029 6.904892e-13
## 439011 -1.3119275 0.7366940 -7.917059 5.203197e-12
## 787677 -1.0103604 0.9571802 -7.692180 1.524001e-11
## 229279 -0.4887401 1.0189108 -7.594301 2.428540e-11
## 844305 -0.6092721 1.7155302 -7.463128 4.525974e-11
## 651916 -0.6704890 0.7854368 -7.154501 1.939281e-10
## 501755 -0.6686642 2.3644606 -7.114250 2.341954e-10
## 446034 -0.6021258 1.0695885 -7.068230 2.904778e-10
## 500621 -0.5526212 1.3349773 -7.016926 3.691473e-10
## 558377 -0.6099956 2.8805567 -6.703254 1.580946e-09
## adj.P.Val B
## 219838 5.923720e-07 18.26786
## 439011 2.231916e-06 16.43705
## 787677 4.358146e-06 15.46114
## 229279 5.208623e-06 15.03766
## 844305 7.765684e-06 14.47158
## 651916 2.772855e-05 13.14730
## 501755 2.870239e-05 12.97548
## 446034 3.115019e-05 12.77932
## 500621 3.518802e-05 12.56100
## 558377 1.356296e-04 11.23533dm_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] 14457length(dm_dn)## [1] 41135confects <- 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 272985 3 5375 959
## NotSig 100214 857891 849399 855966
## Up 484703 8 3128 977dm <- 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
## 812829 cg26222198
## 193227 cg05705212 CXCR6;FYCO1
## 707434 cg22431883 CEP41;CEP41;CEP41;CEP41
## 189285 cg05584657
## 705480 cg22355342 EXOC4
## 649171 cg20260859
## 81324 cg02364564
## 723115 cg23005797 C2orf48
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 501755 -0.8535777 2.01817454 -6.681149 2.651792e-08
## 464536 0.9532704 0.05788172 6.277654 1.077234e-07
## 812829 -0.8294957 0.30623331 -6.135492 1.764273e-07
## 193227 0.9317736 1.71482186 6.125368 1.827319e-07
## 707434 0.7668647 1.41053482 6.099598 1.998121e-07
## 189285 0.8270036 0.67585418 6.077018 2.160807e-07
## 705480 0.7431985 -0.04410430 6.014149 2.686710e-07
## 649171 -0.4789004 1.49831629 -5.968798 3.143540e-07
## 81324 0.6676420 1.64936763 5.956814 3.276680e-07
## 723115 0.7846637 2.48167450 5.949308 3.362911e-07
## adj.P.Val B
## 501755 0.02274978 8.441729
## 464536 0.02391174 7.225805
## 812829 0.02391174 6.796416
## 193227 0.02391174 6.765830
## 707434 0.02391174 6.687974
## 189285 0.02391174 6.619757
## 705480 0.02391174 6.429823
## 649171 0.02391174 6.292828
## 81324 0.02391174 6.256630
## 723115 0.02391174 6.233959dm_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] 977length(dm_dn)## [1] 959confects <- 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 273942 8 6372 0
## NotSig 92858 857884 843330 857902
## Up 491102 10 8200 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.3009477
## 68537 cg01987065 ESYT2 -0.2764491
## 478405 cg14507403 LMTK2 Unclassified 1.4130530
## 839138 cg27123414 0.3749537
## 30698 cg00873371 -0.2930230
## 623821 cg19345678 -0.2836077
## 658377 cg20604021 PTTG1IP 0.3481162
## 124384 cg03634967 SLC29A1;SLC29A1;SLC29A1 Unclassified 0.4674893
## 103641 cg03016486 MIR8081 -0.3403144
## 221316 cg06572904 -0.3848850
## AveExpr t P.Value adj.P.Val B
## 319051 2.586628 -5.095932 6.160339e-06 0.9995568 0.1759891
## 68537 3.182933 -4.719718 2.178844e-05 0.9995568 -0.3728468
## 478405 -1.169378 4.706028 2.280233e-05 0.9995568 -0.3928500
## 839138 1.910324 4.658406 2.670331e-05 0.9995568 -0.4624318
## 30698 2.001459 -4.620528 3.026708e-05 0.9995568 -0.5177671
## 623821 2.166406 -4.581200 3.446038e-05 0.9995568 -0.5752089
## 658377 2.047815 4.579755 3.462486e-05 0.9995568 -0.5773195
## 124384 -4.387072 4.549232 3.828385e-05 0.9995568 -0.6218859
## 103641 2.069271 -4.500378 4.494227e-05 0.9995568 -0.6931843
## 221316 3.625364 -4.475905 4.869125e-05 0.9995568 -0.7288825dm_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 279902 48648 534 41026
## NotSig 60524 778176 834995 783358
## Up 498524 12126 3421 14566dm <- 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
## 637274 cg20361768 LINC00880
## 490426 cg15237707 KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15
## 435911 cg13574263
## 489326 cg15200345
## 545758 cg17051905 CLN8
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 214611 -1.1206545 0.8453538 -8.346078 6.505713e-13
## 429028 -1.3119275 0.7366940 -7.926126 4.913922e-12
## 770258 -1.0103604 0.9571802 -7.700465 1.446527e-11
## 223845 -0.4887401 1.0189108 -7.598077 2.356389e-11
## 825671 -0.6092721 1.7155302 -7.468994 4.351087e-11
## 637274 -0.6704890 0.7854368 -7.160943 1.862323e-10
## 490426 -0.6686642 2.3644606 -7.120670 2.249720e-10
## 435911 -0.6021258 1.0695885 -7.074059 2.798788e-10
## 489326 -0.5526212 1.3349773 -7.022193 3.567066e-10
## 545758 -0.6099956 2.8805567 -6.709143 1.525669e-09
## adj.P.Val B
## 214611 5.457968e-07 18.32744
## 429028 2.061267e-06 16.49379
## 770258 4.045213e-06 15.51294
## 223845 4.942232e-06 15.06929
## 825671 7.300689e-06 14.51140
## 637274 2.603994e-05 13.18761
## 490426 2.696289e-05 13.01545
## 435911 2.935054e-05 12.81647
## 489326 3.325100e-05 12.59544
## 545758 1.279960e-04 11.27048dm_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] 14566length(dm_dn)## [1] 41026confects <- 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 269813 2 12 1053
## NotSig 91500 838940 838906 836850
## Up 477637 8 32 1047dm <- 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
## 184810 cg05584657
## 689746 cg22355342 EXOC4
## 634597 cg20260859
## 79352 cg02364564
## 707001 cg23005797 C2orf48
## 389265 cg12048331 WDR82
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 490426 -0.8535777 2.01817454 -6.690660 2.510338e-08
## 453993 0.9532704 0.05788172 6.288405 1.018670e-07
## 188665 0.9317736 1.71482186 6.135874 1.731537e-07
## 691656 0.7668647 1.41053482 6.108084 1.907110e-07
## 184810 0.8270036 0.67585418 6.086382 2.056462e-07
## 689746 0.7431985 -0.04410430 6.022301 2.569014e-07
## 634597 -0.4789004 1.49831629 -5.968246 3.099047e-07
## 79352 0.6676420 1.64936763 5.963530 3.150165e-07
## 707001 0.7846637 2.48167450 5.958140 3.209608e-07
## 389265 1.1185390 0.56815119 5.941602 3.399058e-07
## adj.P.Val B
## 490426 0.02106048 8.494001
## 453993 0.02453424 7.278234
## 188665 0.02453424 6.816187
## 691656 0.02453424 6.731987
## 184810 0.02453424 6.666236
## 689746 0.02453424 6.472086
## 634597 0.02453424 6.308340
## 79352 0.02453424 6.294054
## 707001 0.02453424 6.277729
## 389265 0.02453424 6.227642dm_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] 1047length(dm_dn)## [1] 1053confects <- 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 270597 8 156 0
## NotSig 84679 838932 838097 838950
## Up 483674 10 697 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.3009477
## 467564 cg14507403 LMTK2 Unclassified 1.4130530
## 66884 cg01987065 ESYT2 -0.2764491
## 820601 cg27123414 0.3749537
## 29973 cg00873371 -0.2930230
## 643605 cg20604021 PTTG1IP 0.3481162
## 609757 cg19345678 -0.2836077
## 121467 cg03634967 SLC29A1;SLC29A1;SLC29A1 Unclassified 0.4674893
## 101150 cg03016486 MIR8081 -0.3403144
## 216056 cg06572904 -0.3848850
## AveExpr t P.Value adj.P.Val B
## 311708 2.586628 -5.086085 6.306186e-06 0.9999373 0.1923354
## 467564 -1.169378 4.716211 2.186842e-05 0.9999373 -0.3517234
## 66884 3.182933 -4.710268 2.230521e-05 0.9999373 -0.3604772
## 820601 1.910324 4.658539 2.648666e-05 0.9999373 -0.4366605
## 29973 2.001459 -4.614201 3.067558e-05 0.9999373 -0.5019445
## 643605 2.047815 4.578612 3.450197e-05 0.9999373 -0.5543318
## 609757 2.166406 -4.574094 3.501998e-05 0.9999373 -0.5609807
## 121467 -4.387072 4.553399 3.749191e-05 0.9999373 -0.5914343
## 101150 2.069271 -4.499135 4.481335e-05 0.9999373 -0.6712471
## 216056 3.625364 -4.477305 4.813820e-05 0.9999373 -0.7033372dm_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 283900 57501 7937 41135
## NotSig 66029 784148 835091 802310
## Up 507973 16253 14874 14457dm <- 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
## 651916 cg20361768 LINC00880
## 501755 cg15237707 KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15
## 446034 cg13574263
## 500621 cg15200345
## 558377 cg17051905 CLN8
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 219838 -1.1206545 0.8453538 -8.337029 6.904892e-13
## 439011 -1.3119275 0.7366940 -7.917059 5.203197e-12
## 787677 -1.0103604 0.9571802 -7.692180 1.524001e-11
## 229279 -0.4887401 1.0189108 -7.594301 2.428540e-11
## 844305 -0.6092721 1.7155302 -7.463128 4.525974e-11
## 651916 -0.6704890 0.7854368 -7.154501 1.939281e-10
## 501755 -0.6686642 2.3644606 -7.114250 2.341954e-10
## 446034 -0.6021258 1.0695885 -7.068230 2.904778e-10
## 500621 -0.5526212 1.3349773 -7.016926 3.691473e-10
## 558377 -0.6099956 2.8805567 -6.703254 1.580946e-09
## adj.P.Val B
## 219838 5.923720e-07 18.26786
## 439011 2.231916e-06 16.43705
## 787677 4.358146e-06 15.46114
## 229279 5.208623e-06 15.03766
## 844305 7.765684e-06 14.47158
## 651916 2.772855e-05 13.14730
## 501755 2.870239e-05 12.97548
## 446034 3.115019e-05 12.77932
## 500621 3.518802e-05 12.56100
## 558377 1.356296e-04 11.23533dm_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] 14457length(dm_dn)## [1] 41135confects <- 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 272985 3 5375 959
## NotSig 100214 857891 849399 855966
## Up 484703 8 3128 977dm <- 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
## 812829 cg26222198
## 193227 cg05705212 CXCR6;FYCO1
## 707434 cg22431883 CEP41;CEP41;CEP41;CEP41
## 189285 cg05584657
## 705480 cg22355342 EXOC4
## 649171 cg20260859
## 81324 cg02364564
## 723115 cg23005797 C2orf48
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 501755 -0.8535777 2.01817454 -6.681149 2.651792e-08
## 464536 0.9532704 0.05788172 6.277654 1.077234e-07
## 812829 -0.8294957 0.30623331 -6.135492 1.764273e-07
## 193227 0.9317736 1.71482186 6.125368 1.827319e-07
## 707434 0.7668647 1.41053482 6.099598 1.998121e-07
## 189285 0.8270036 0.67585418 6.077018 2.160807e-07
## 705480 0.7431985 -0.04410430 6.014149 2.686710e-07
## 649171 -0.4789004 1.49831629 -5.968798 3.143540e-07
## 81324 0.6676420 1.64936763 5.956814 3.276680e-07
## 723115 0.7846637 2.48167450 5.949308 3.362911e-07
## adj.P.Val B
## 501755 0.02274978 8.441729
## 464536 0.02391174 7.225805
## 812829 0.02391174 6.796416
## 193227 0.02391174 6.765830
## 707434 0.02391174 6.687974
## 189285 0.02391174 6.619757
## 705480 0.02391174 6.429823
## 649171 0.02391174 6.292828
## 81324 0.02391174 6.256630
## 723115 0.02391174 6.233959dm_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] 977length(dm_dn)## [1] 959confects <- 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 273942 8 6372 0
## NotSig 92858 857884 843330 857902
## Up 491102 10 8200 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.3009477
## 68537 cg01987065 ESYT2 -0.2764491
## 478405 cg14507403 LMTK2 Unclassified 1.4130530
## 839138 cg27123414 0.3749537
## 30698 cg00873371 -0.2930230
## 623821 cg19345678 -0.2836077
## 658377 cg20604021 PTTG1IP 0.3481162
## 124384 cg03634967 SLC29A1;SLC29A1;SLC29A1 Unclassified 0.4674893
## 103641 cg03016486 MIR8081 -0.3403144
## 221316 cg06572904 -0.3848850
## AveExpr t P.Value adj.P.Val B
## 319051 2.586628 -5.095932 6.160339e-06 0.9995568 0.1759891
## 68537 3.182933 -4.719718 2.178844e-05 0.9995568 -0.3728468
## 478405 -1.169378 4.706028 2.280233e-05 0.9995568 -0.3928500
## 839138 1.910324 4.658406 2.670331e-05 0.9995568 -0.4624318
## 30698 2.001459 -4.620528 3.026708e-05 0.9995568 -0.5177671
## 623821 2.166406 -4.581200 3.446038e-05 0.9995568 -0.5752089
## 658377 2.047815 4.579755 3.462486e-05 0.9995568 -0.5773195
## 124384 -4.387072 4.549232 3.828385e-05 0.9995568 -0.6218859
## 103641 2.069271 -4.500378 4.494227e-05 0.9995568 -0.6931843
## 221316 3.625364 -4.475905 4.869125e-05 0.9995568 -0.7288825dm_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 279902 48648 534 41026
## NotSig 60524 778176 834995 783358
## Up 498524 12126 3421 14566dm <- 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
## 637274 cg20361768 LINC00880
## 490426 cg15237707 KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15;KCNJ15
## 435911 cg13574263
## 489326 cg15200345
## 545758 cg17051905 CLN8
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 214611 -1.1206545 0.8453538 -8.346078 6.505713e-13
## 429028 -1.3119275 0.7366940 -7.926126 4.913922e-12
## 770258 -1.0103604 0.9571802 -7.700465 1.446527e-11
## 223845 -0.4887401 1.0189108 -7.598077 2.356389e-11
## 825671 -0.6092721 1.7155302 -7.468994 4.351087e-11
## 637274 -0.6704890 0.7854368 -7.160943 1.862323e-10
## 490426 -0.6686642 2.3644606 -7.120670 2.249720e-10
## 435911 -0.6021258 1.0695885 -7.074059 2.798788e-10
## 489326 -0.5526212 1.3349773 -7.022193 3.567066e-10
## 545758 -0.6099956 2.8805567 -6.709143 1.525669e-09
## adj.P.Val B
## 214611 5.457968e-07 18.32744
## 429028 2.061267e-06 16.49379
## 770258 4.045213e-06 15.51294
## 223845 4.942232e-06 15.06929
## 825671 7.300689e-06 14.51140
## 637274 2.603994e-05 13.18761
## 490426 2.696289e-05 13.01545
## 435911 2.935054e-05 12.81647
## 489326 3.325100e-05 12.59544
## 545758 1.279960e-04 11.27048dm_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] 14566length(dm_dn)## [1] 41026confects <- 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
## 306 Formation of a pool of free 40S subunits
## 418 Influenza Viral RNA Transcription and Replication
## 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
## 140 Chromatin modifying enzymes
## 141 Chromatin organization
## 1096 Viral mRNA Translation
## 608 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)
## setSize pANOVA s.dist p.adjustANOVA
## 510 542 6.102430e-07 0.12670035 0.0006920155
## 855 85 6.762840e-06 0.28282937 0.0038345300
## 417 129 1.185932e-05 0.22386005 0.0042580150
## 306 76 1.653991e-05 0.28611839 0.0042580150
## 418 111 1.877432e-05 0.23557130 0.0042580150
## 658 67 2.459173e-05 0.29834768 0.0045682471
## 274 70 2.945412e-05 0.28909735 0.0045682471
## 341 86 3.222749e-05 0.25972565 0.0045682471
## 473 85 5.247361e-05 0.25415413 0.0056337860
## 276 70 5.475349e-05 0.27918837 0.0056337860
## 132 400 5.916332e-05 0.11800137 0.0056337860
## 110 93 7.038642e-05 0.23889809 0.0056337860
## 275 93 7.038642e-05 0.23889809 0.0056337860
## 133 395 7.151396e-05 0.11739154 0.0056337860
## 522 1288 7.452098e-05 0.06765994 0.0056337860
## 501 149 8.741466e-05 0.18674178 0.0058805839
## 140 155 9.793846e-05 0.18185955 0.0058805839
## 141 155 9.793846e-05 0.18185955 0.0058805839
## 1096 68 9.852830e-05 0.27340554 0.0058805839
## 608 72 1.619077e-04 0.25742019 0.0091801642unlink("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 269813 2 12 1053
## NotSig 91500 838940 838906 836850
## Up 477637 8 32 1047dm <- 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
## 184810 cg05584657
## 689746 cg22355342 EXOC4
## 634597 cg20260859
## 79352 cg02364564
## 707001 cg23005797 C2orf48
## 389265 cg12048331 WDR82
## Regulatory_Feature_Group logFC AveExpr t P.Value
## 490426 -0.8535777 2.01817454 -6.690660 2.510338e-08
## 453993 0.9532704 0.05788172 6.288405 1.018670e-07
## 188665 0.9317736 1.71482186 6.135874 1.731537e-07
## 691656 0.7668647 1.41053482 6.108084 1.907110e-07
## 184810 0.8270036 0.67585418 6.086382 2.056462e-07
## 689746 0.7431985 -0.04410430 6.022301 2.569014e-07
## 634597 -0.4789004 1.49831629 -5.968246 3.099047e-07
## 79352 0.6676420 1.64936763 5.963530 3.150165e-07
## 707001 0.7846637 2.48167450 5.958140 3.209608e-07
## 389265 1.1185390 0.56815119 5.941602 3.399058e-07
## adj.P.Val B
## 490426 0.02106048 8.494001
## 453993 0.02453424 7.278234
## 188665 0.02453424 6.816187
## 691656 0.02453424 6.731987
## 184810 0.02453424 6.666236
## 689746 0.02453424 6.472086
## 634597 0.02453424 6.308340
## 79352 0.02453424 6.294054
## 707001 0.02453424 6.277729
## 389265 0.02453424 6.227642dm_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] 1047length(dm_dn)## [1] 1053confects <- 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
## 350 Gene expression (Transcription)
## 983 TCF dependent signaling in response to WNT
## 510 Metabolism of RNA
## 417 Influenza Infection
## 348 Gene Silencing by RNA
## 418 Influenza Viral RNA Transcription and Replication
## 809 Regulation of expression of SLITs and ROBOs
## 415 Infectious disease
## 644 PIP3 activates AKT signaling
## 133 Cellular responses to stress
## 132 Cellular responses to external stimuli
## 746 RNA Polymerase II Transcription
## 916 Signaling by NOTCH
## 81 Beta-catenin independent WNT signaling
## 453 Interleukin-4 and Interleukin-13 signaling
## 607 Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC)
## 609 Nonsense-Mediated Decay (NMD)
## 383 HIV Infection
## 352 Generic Transcription Pathway
## setSize pANOVA s.dist p.adjustANOVA
## 944 178 1.760953e-05 0.18720234 0.01251486
## 350 1025 2.736340e-05 0.07927544 0.01251486
## 983 119 3.310809e-05 0.22079158 0.01251486
## 510 542 6.959681e-05 0.10105787 0.01973069
## 417 129 1.266716e-04 0.19590729 0.02868111
## 348 59 1.517519e-04 0.28542037 0.02868111
## 418 111 2.154997e-04 0.20372331 0.03484119
## 809 122 2.884195e-04 0.19049798 0.03484119
## 415 490 3.092984e-04 0.09616634 0.03484119
## 644 182 3.202473e-04 0.15521536 0.03484119
## 133 395 3.586641e-04 0.10550261 0.03484119
## 132 400 3.686898e-04 0.10465269 0.03484119
## 746 913 4.321998e-04 0.07011376 0.03770112
## 916 133 6.094102e-04 0.17257101 0.04936223
## 81 107 7.879484e-04 0.18823109 0.05317974
## 453 48 8.567161e-04 0.27836660 0.05317974
## 607 86 8.594108e-04 0.20825948 0.05317974
## 609 86 8.594108e-04 0.20825948 0.05317974
## 383 188 9.087755e-04 0.14083615 0.05317974
## 352 816 9.379143e-04 0.06939730 0.05317974unlink("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 270597 8 156 0
## NotSig 84679 838932 838097 838950
## Up 483674 10 697 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.3009477
## 467564 cg14507403 LMTK2 Unclassified 1.4130530
## 66884 cg01987065 ESYT2 -0.2764491
## 820601 cg27123414 0.3749537
## 29973 cg00873371 -0.2930230
## 643605 cg20604021 PTTG1IP 0.3481162
## 609757 cg19345678 -0.2836077
## 121467 cg03634967 SLC29A1;SLC29A1;SLC29A1 Unclassified 0.4674893
## 101150 cg03016486 MIR8081 -0.3403144
## 216056 cg06572904 -0.3848850
## AveExpr t P.Value adj.P.Val B
## 311708 2.586628 -5.086085 6.306186e-06 0.9999373 0.1923354
## 467564 -1.169378 4.716211 2.186842e-05 0.9999373 -0.3517234
## 66884 3.182933 -4.710268 2.230521e-05 0.9999373 -0.3604772
## 820601 1.910324 4.658539 2.648666e-05 0.9999373 -0.4366605
## 29973 2.001459 -4.614201 3.067558e-05 0.9999373 -0.5019445
## 643605 2.047815 4.578612 3.450197e-05 0.9999373 -0.5543318
## 609757 2.166406 -4.574094 3.501998e-05 0.9999373 -0.5609807
## 121467 -4.387072 4.553399 3.749191e-05 0.9999373 -0.5914343
## 101150 2.069271 -4.499135 4.481335e-05 0.9999373 -0.6712471
## 216056 3.625364 -4.477305 4.813820e-05 0.9999373 -0.7033372dm_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
## 698 Processing of Capped Intronless Pre-mRNA
## 541 Mitochondrial tRNA aminoacylation
## 909 Signaling by Hippo
## 497 MET promotes cell motility
## 1045 Transcriptional activity of SMAD2/SMAD3:SMAD4 heterotrimer
## 851 SHC1 events in ERBB2 signaling
## 685 Post-translational modification: synthesis of GPI-anchored proteins
## 991 TNFR1-induced NFkappaB signaling pathway
## 555 Muscle contraction
## setSize pANOVA s.dist p.adjustANOVA
## 866 19 0.004119427 0.38028182 0.9693997
## 221 16 0.007204634 0.38812313 0.9693997
## 1000 16 0.007915160 -0.38356526 0.9693997
## 396 21 0.007971496 0.33457531 0.9693997
## 217 26 0.009623496 0.29347944 0.9693997
## 516 435 0.013703213 0.06957665 0.9693997
## 983 119 0.016633217 0.12740828 0.9693997
## 944 178 0.020554471 0.10101919 0.9693997
## 367 63 0.022273914 0.16671931 0.9693997
## 1054 234 0.025308406 -0.08529970 0.9693997
## 945 23 0.029420297 0.26244439 0.9693997
## 698 21 0.029779335 0.27402876 0.9693997
## 541 18 0.030840380 -0.29405489 0.9693997
## 909 12 0.032470036 0.35662088 0.9693997
## 497 13 0.034228359 0.33924858 0.9693997
## 1045 37 0.034926234 -0.20052657 0.9693997
## 851 11 0.036547134 0.36414057 0.9693997
## 685 31 0.039649069 0.21363607 0.9693997
## 991 20 0.043591036 -0.26075771 0.9693997
## 555 69 0.044572993 -0.14004019 0.9693997unlink("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)