DNA meth and IVF
The GMEA team and ART epigenetics team
2024-01-04
Introduction
Here we are looking at putting together the easiest possible workflow for the package.
We will analyze methylation data for IVF conceived infants compared to naturally conceived.
Data from:
Estill MS, Bolnick JM, Waterland RA, Bolnick AD, Diamond MP, Krawetz SA. Assisted reproductive technology alters deoxyribonucleic acid methylation profiles in bloodspots of newborn infants. Fertil Steril. 2016 Sep 1;106(3):629-639.e10. doi: 10.1016/j.fertnstert.2016.05.006. Epub 2016 Jun 8. PMID: 27288894.
Novakovic B, Lewis S, Halliday J, Kennedy J, Burgner DP, Czajko A, Kim B, Sexton-Oates A, Juonala M, Hammarberg K, Amor DJ, Doyle LW, Ranganathan S, Welsh L, Cheung M, McBain J, McLachlan R, Saffery R. Assisted reproductive technologies are associated with limited epigenetic variation at birth that largely resolves by adulthood. Nat Commun. 2019 Sep 2;10(1):3922. doi: 10.1038/s41467-019-11929-9. PMID: 31477727; PMCID: PMC6718382.
Requirements
suppressPackageStartupMessages({
library("limma")
library("eulerr")
library("IlluminaHumanMethylation450kmanifest")
library("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
library("HGNChelper")
library("tictoc")
library("mitch")
library("gplots")
library("kableExtra")
library("beeswarm")
library("missMethyl")
library("gridExtra")
library("png")
})
CORES=16Load data
Reactome pathways were downloaded on the 14th Sept 2023 from MsigDB.
We’ll also score the probes by the signed -log pvalue
# estill 2016 450k
art1 <- readRDS("estill_nat_vs_fh.rds")
art1 <- art1[[1]]
rownames(art1) <- art1$Row.names
str(art1)## 'data.frame': 418833 obs. of 11 variables:
## $ Row.names : 'AsIs' chr "cg20757019" "cg08260549" "cg16905586" "cg01931797" ...
## $ UCSC_RefGene_Name : chr "" "XPO7;XPO7" "GMPPB;GMPPB;GMPPB;GMPPB" "KCTD21;USP35" ...
## $ Regulatory_Feature_Group: chr "Unclassified" "Promoter_Associated" "Promoter_Associated" "Promoter_Associated" ...
## $ Islands_Name : chr "chr2:95872891-95873548" "chr8:21776842-21777942" "chr3:49760863-49761895" "chr11:77899564-77899990" ...
## $ Relation_to_Island : chr "Island" "Island" "Island" "Island" ...
## $ logFC : num -1.224 -0.993 -0.742 -1.018 -0.893 ...
## $ AveExpr : num -4.12 -4.16 -3.59 -4.48 -4.32 ...
## $ t : num -22.3 -19 -19 -18.5 -18.4 ...
## $ P.Value : num 8.34e-37 5.31e-32 6.93e-32 3.92e-31 5.01e-31 ...
## $ adj.P.Val : num 3.49e-31 9.68e-27 9.68e-27 3.82e-26 3.82e-26 ...
## $ B : num 71.8 61.5 61.2 59.6 59.3 ...# novakovic 2019 EPIC
art2 <- readRDS("novakovic_nat_vs_fh.rds")
art2 <-art2[[1]]
rownames(art2) <- art2$Row.names
str(art2)## 'data.frame': 793844 obs. of 11 variables:
## $ Row.names : 'AsIs' chr "cg25867694" "cg16794867" "cg00352360" "cg18789663" ...
## $ UCSC_RefGene_Name : chr "SELM" "AADACL2;MIR548H2" "RPL23;SNORA21" "PLD5" ...
## $ Regulatory_Feature_Group: chr "Promoter_Associated" "" "" "" ...
## $ Islands_Name : chr "chr22:31503236-31503871" "" "chr17:37009471-37010118" "chr1:242687922-242688682" ...
## $ Relation_to_Island : chr "Island" "OpenSea" "N_Shore" "Island" ...
## $ logFC : num -0.32 -0.32 -0.385 -0.335 -0.243 ...
## $ AveExpr : num -3.154 -0.905 -2.665 -2.748 -2.064 ...
## $ t : num -6.74 -6.64 -6.36 -6.34 -6.22 ...
## $ P.Value : num 4.00e-10 6.58e-10 2.75e-09 3.15e-09 5.65e-09 ...
## $ adj.P.Val : num 0.000261 0.000261 0.000625 0.000625 0.000633 ...
## $ B : num 12.26 11.81 10.53 10.41 9.89 ...gs_symbols <- gmt_import("c2.cp.reactome.v2023.1.Hs.symbols.gmt")Curate the annotation
Use all probes on the chip. Start with 450k chip.
Update gene symbols.
tic()
anno <- getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
myann <- data.frame(anno[,c("UCSC_RefGene_Name","UCSC_RefGene_Group","Islands_Name","Relation_to_Island")])
gp <- myann[,"UCSC_RefGene_Name",drop=FALSE]
gp2 <- strsplit(gp$UCSC_RefGene_Name,";")
names(gp2) <- rownames(gp)
gp2 <- lapply(gp2,unique)
gt1 <- stack(gp2)
colnames(gt1) <- c("gene","probe")
gt1$probe <- as.character(gt1$probe)
dim(gt1)## [1] 407090 2str(gt1)## 'data.frame': 407090 obs. of 2 variables:
## $ gene : chr "TSPY4" "FAM197Y2" "TTTY14" "TMSB4Y" ...
## $ probe: chr "cg00050873" "cg00050873" "cg00212031" "cg00214611" ...length(unique(gt1$gene))## [1] 21231toc() #6.3s## 4.559 sec elapsedtic()
#new.hgnc.table <- getCurrentHumanMap()
new.hgnc.table <- readRDS("new.hgnc.table.rds")
fix <- checkGeneSymbols(gt1$gene,map=new.hgnc.table)## Warning in checkGeneSymbols(gt1$gene, map = new.hgnc.table): Human gene symbols
## should be all upper-case except for the 'orf' in open reading frames. The case
## of some letters was corrected.## Warning in checkGeneSymbols(gt1$gene, map = new.hgnc.table): x contains
## non-approved gene symbolsfix2 <- fix[which(fix$x != fix$Suggested.Symbol),]
length(unique(fix2$x))## [1] 2788gt1$gene <- fix$Suggested.Symbol
toc()## 35.226 sec elapsedNow EPIC chip.
tic()
anno <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
myann <- data.frame(anno[,c("UCSC_RefGene_Name","UCSC_RefGene_Group","Islands_Name","Relation_to_Island")])
gp <- myann[,"UCSC_RefGene_Name",drop=FALSE]
gp2 <- strsplit(gp$UCSC_RefGene_Name,";")
names(gp2) <- rownames(gp)
gp2 <- lapply(gp2,unique)
gt2 <- stack(gp2)
colnames(gt2) <- c("gene","probe")
gt2$probe <- as.character(gt2$probe)
dim(gt2)## [1] 684970 2str(gt2)## 'data.frame': 684970 obs. of 2 variables:
## $ gene : chr "YTHDF1" "EIF2S3" "PKN3" "CCDC57" ...
## $ probe: chr "cg18478105" "cg09835024" "cg14361672" "cg01763666" ...length(unique(gt2$gene))## [1] 27364toc() #10.6s## 6.469 sec elapsedtic()
#new.hgnc.table <- getCurrentHumanMap()
fix <- checkGeneSymbols(gt2$gene,map=new.hgnc.table)## Warning in checkGeneSymbols(gt2$gene, map = new.hgnc.table): Human gene symbols
## should be all upper-case except for the 'orf' in open reading frames. The case
## of some letters was corrected.## Warning in checkGeneSymbols(gt2$gene, map = new.hgnc.table): x contains
## non-approved gene symbolsfix2 <- fix[which(fix$x != fix$Suggested.Symbol),]
length(unique(fix2$x))## [1] 3253gt2$gene <- fix$Suggested.Symbol
toc()## 44.852 sec elapsedImport for mitch
m1 <- mitch_import(x=art1,DEtype="limma",geneTable=gt1)## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 418833## Note: no. genes in output = 19240## Warning in mitch_import(x = art1, DEtype = "limma", geneTable = gt1): Warning: less than half of the input genes are also in the
## outputdim(m1)## [1] 19240 1head(m1)## x
## A1BG -2.7341135
## A1BG-AS1 -2.9086004
## A1CF -1.4455477
## A2M 0.1755092
## A2ML1 1.1420743
## A4GALT -1.2155001m2 <- mitch_import(x=art2,DEtype="limma",geneTable=gt2)## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 793844## Note: no. genes in output = 22588## Warning in mitch_import(x = art2, DEtype = "limma", geneTable = gt2): Warning: less than half of the input genes are also in the
## outputdim(m2)## [1] 22588 1head(m2)## x
## A1BG -0.7186568
## A1BG-AS1 -0.1813097
## A1CF -0.5404743
## A2M -0.4595169
## A2M-AS1 1.6811677
## A2ML1 -0.4082271m <- merge(m1,m2,by=0)
rownames(m) <- m[,1]
m[,1]=NULL
colnames(m) <- c("Estill","Novakovic")
dim(m)## [1] 19234 2head(m)## Estill Novakovic
## A1BG -2.7341135 -0.7186568
## A1BG-AS1 -2.9086004 -0.1813097
## A1CF -1.4455477 -0.5404743
## A2M 0.1755092 -0.4595169
## A2ML1 1.1420743 -0.4082271
## A4GALT -1.2155001 -0.5418459Run mitch
tic()
mres <- mitch_calc(x=m,genesets=gs_symbols,minsetsize=5,cores=16, priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.toc() #11.45s## 50.848 sec elapsedmtable <- mres$enrichment_result
sig <- subset(mtable,p.adjustMANOVA<0.05)
head(sig,20)## set setSize
## 145 REACTOME_AMINO_ACID_CONJUGATION 8
## 1144 REACTOME_MRNA_EDITING_C_TO_U_CONVERSION 8
## 657 REACTOME_ADRENOCEPTORS 9
## 1183 REACTOME_MRNA_EDITING 10
## 110 REACTOME_BETA_DEFENSINS 32
## 646 REACTOME_PD_1_SIGNALING 21
## 651 REACTOME_FOLDING_OF_ACTIN_BY_CCT_TRIC 10
## 111 REACTOME_DEFENSINS 41
## 180 REACTOME_FORMATION_OF_ATP_BY_CHEMIOSMOTIC_COUPLING 16
## 1060 REACTOME_ANTIMICROBIAL_PEPTIDES 81
## 1479 REACTOME_FOXO_MEDIATED_TRANSCRIPTION_OF_CELL_CYCLE_GENES 15
## 628 REACTOME_OLFACTORY_SIGNALING_PATHWAY 313
## 394 REACTOME_GENERATION_OF_SECOND_MESSENGER_MOLECULES 30
## 1304 REACTOME_CRISTAE_FORMATION 27
## 214 REACTOME_INITIAL_TRIGGERING_OF_COMPLEMENT 21
## 1158 REACTOME_RNA_POLYMERASE_III_CHAIN_ELONGATION 18
## 660 REACTOME_SIGNAL_REGULATORY_PROTEIN_FAMILY_INTERACTIONS 16
## 1082 REACTOME_FORMATION_OF_THE_CORNIFIED_ENVELOPE 123
## 212 REACTOME_COMPLEMENT_CASCADE 54
## 912 REACTOME_DECTIN_2_FAMILY 23
## pMANOVA s.Estill s.Novakovic p.Estill p.Novakovic s.dist
## 145 1.296663e-03 0.71026215 0.3976776 5.029933e-04 5.144144e-02 0.8140146
## 1144 4.144411e-03 0.47281026 0.5886430 2.056827e-02 3.935646e-03 0.7550166
## 657 4.987003e-03 -0.42261523 -0.5557867 2.813103e-02 3.884375e-03 0.6982138
## 1183 4.489133e-03 0.41268206 0.5274865 2.383564e-02 3.870387e-03 0.6697376
## 110 1.112075e-07 0.44740782 0.4682943 1.182208e-05 4.529767e-06 0.6476676
## 646 6.850562e-05 0.47535027 0.3931044 1.623552e-04 1.816773e-03 0.6168378
## 651 5.213322e-03 -0.14824698 -0.5920828 4.169493e-01 1.185440e-03 0.6103599
## 111 2.563631e-08 0.34054463 0.4843115 1.611358e-04 7.996533e-08 0.5920543
## 180 1.629547e-03 -0.34700866 -0.4599334 1.625360e-02 1.445543e-03 0.5761543
## 1060 1.200282e-14 0.28955075 0.4851620 6.636331e-06 4.283235e-14 0.5649972
## 1479 3.949713e-03 -0.36231854 -0.4205387 1.511603e-02 4.801356e-03 0.5550924
## 628 1.528217e-48 0.46289575 0.2647252 3.580471e-45 8.146077e-16 0.5332466
## 394 3.848302e-05 0.38586753 0.3677116 2.539894e-04 4.902626e-04 0.5330155
## 1304 2.039429e-04 -0.29481921 -0.4148275 8.011615e-03 1.905250e-04 0.5089207
## 214 5.372769e-04 0.08043909 0.4869248 5.234224e-01 1.118959e-04 0.4935242
## 1158 4.690039e-03 -0.20394752 -0.4355167 1.341470e-01 1.378569e-03 0.4809048
## 660 4.700215e-03 0.07628916 0.4705094 5.972950e-01 1.119371e-03 0.4766541
## 1082 9.207837e-16 0.22802983 0.4156781 1.259184e-05 1.647017e-15 0.4741159
## 212 4.491245e-07 0.23598077 0.4024949 2.704396e-03 3.102100e-07 0.4665716
## 912 1.516021e-03 0.16882764 0.4297561 1.610665e-01 3.596542e-04 0.4617283
## SD p.adjustMANOVA
## 145 0.22103063 1.717032e-02
## 1144 0.08190611 3.866547e-02
## 657 0.09416648 4.402505e-02
## 1183 0.08117898 4.050086e-02
## 110 0.01476898 8.695368e-06
## 646 0.05815663 1.648963e-03
## 651 0.31383933 4.501974e-02
## 111 0.10165851 2.550148e-06
## 180 0.07984984 1.938924e-02
## 1060 0.13831808 2.463579e-12
## 1479 0.04116787 3.705959e-02
## 628 0.14012773 2.509332e-45
## 394 0.01283819 1.089468e-03
## 1304 0.08485870 3.986597e-03
## 214 0.28742879 8.649105e-03
## 1158 0.16374414 4.171758e-02
## 660 0.27875582 4.171758e-02
## 1082 0.13268737 2.159896e-13
## 212 0.11774326 2.458208e-05
## 912 0.18450425 1.911850e-02Rank-rank
mrank <- mres$ranked_profile
palette <- colorRampPalette(c("white", "yellow", "orange", "red", "darkred", "black"))
k <- MASS::kde2d(mrank[,1], mrank[,2])
X_AXIS = "Estill 450k"
Y_AXIS = "Novakovic EPIC"
filled.contour(k, color.palette = palette, plot.title = {
abline(v = 0, h = 0, lty = 2, lwd = 2, col = "blue")
title(main = "rank-rank plot", xlab = X_AXIS,
ylab = Y_AXIS)
})
pdf("fig7a.pdf")
palette <- colorRampPalette(c("white", "yellow", "orange", "red", "darkred", "black"))
k <- MASS::kde2d(mrank[,1], mrank[,2])
X_AXIS = "Estill 450k"
Y_AXIS = "Novakovic EPIC"
filled.contour(k, color.palette = palette, plot.title = {
abline(v = 0, h = 0, lty = 2, lwd = 2, col = "blue")
title(main = "rank-rank plot", xlab = X_AXIS,
ylab = Y_AXIS)
})
dev.off()## pdf
## 2PW Scatter
head(mtable)## set setSize
## 145 REACTOME_AMINO_ACID_CONJUGATION 8
## 1144 REACTOME_MRNA_EDITING_C_TO_U_CONVERSION 8
## 1221 REACTOME_PTK6_REGULATES_PROTEINS_INVOLVED_IN_RNA_PROCESSING 5
## 657 REACTOME_ADRENOCEPTORS 9
## 256 REACTOME_ACTIVATION_OF_C3_AND_C5 5
## 1183 REACTOME_MRNA_EDITING 10
## pMANOVA s.Estill s.Novakovic p.Estill p.Novakovic s.dist
## 145 0.001296663 0.7102621 0.3976776 0.0005029933 0.051441438 0.8140146
## 1144 0.004144411 0.4728103 0.5886430 0.0205682664 0.003935646 0.7550166
## 1221 0.024853952 -0.6950855 -0.0837381 0.0071071085 0.745746126 0.7001114
## 657 0.004987003 -0.4226152 -0.5557867 0.0281310336 0.003884375 0.6982138
## 256 0.016055339 -0.1023766 0.6842894 0.6917904614 0.008049976 0.6919053
## 1183 0.004489133 0.4126821 0.5274865 0.0238356393 0.003870387 0.6697376
## SD p.adjustMANOVA
## 145 0.22103063 0.01717032
## 1144 0.08190611 0.03866547
## 1221 0.43228792 0.12282708
## 657 0.09416648 0.04402505
## 256 0.55625684 0.09727995
## 1183 0.08117898 0.04050086sig <- subset(mtable,p.adjustMANOVA<0.05)
plot(mtable$s.Estill,mtable$s.Novakovic,pch=19,col="gray",
xlab="Estill 450k",ylab="Novakovic EPIC",main="pathway enrichment score")
points(sig$s.Estill,sig$s.Novakovic,pch=19,col="red")
abline(v = 0, h = 0, lty = 2, lwd = 2, col = "blue")
pdf("fig7b.pdf")
plot(mtable$s.Estill,mtable$s.Novakovic,pch=19,col="gray",
xlab="Estill 450k",ylab="Novakovic EPIC",main="pathway enrichment score")
points(sig$s.Estill,sig$s.Novakovic,pch=19,col="red")
abline(v = 0, h = 0, lty = 2, lwd = 2, col = "blue")
dev.off()## pdf
## 2Heatmap
mt <- head(mtable,30)
rownames(mt) <- mt$set
mt <- mt[,c("s.Estill","s.Novakovic")]
rownames(mt) <- gsub("REACTOME_","",rownames(mt))
rownames(mt) <- gsub("_"," ",rownames(mt))
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2(as.matrix(mt),scale="none",trace="none",margins=c(6,25),col=colfunc(25),cexRow=0.6,cexCol=1)
pdf("fig7c.pdf")
heatmap.2(as.matrix(mt),scale="none",trace="none",margins=c(6,25),col=colfunc(25),cexRow=0.6,cexCol=1)
dev.off()## pdf
## 2Example 2d
tic()
mres <- mitch_calc(x=m,genesets=gs_symbols,minsetsize=5,cores=16, priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.toc() #11.45s## 53.388 sec elapsedmtable <- mres$enrichment_result
sig <- subset(mtable,p.adjustMANOVA<0.05)
head(sig,20)## set setSize
## 145 REACTOME_AMINO_ACID_CONJUGATION 8
## 1144 REACTOME_MRNA_EDITING_C_TO_U_CONVERSION 8
## 657 REACTOME_ADRENOCEPTORS 9
## 1183 REACTOME_MRNA_EDITING 10
## 110 REACTOME_BETA_DEFENSINS 32
## 646 REACTOME_PD_1_SIGNALING 21
## 651 REACTOME_FOLDING_OF_ACTIN_BY_CCT_TRIC 10
## 111 REACTOME_DEFENSINS 41
## 180 REACTOME_FORMATION_OF_ATP_BY_CHEMIOSMOTIC_COUPLING 16
## 1060 REACTOME_ANTIMICROBIAL_PEPTIDES 81
## 1479 REACTOME_FOXO_MEDIATED_TRANSCRIPTION_OF_CELL_CYCLE_GENES 15
## 628 REACTOME_OLFACTORY_SIGNALING_PATHWAY 313
## 394 REACTOME_GENERATION_OF_SECOND_MESSENGER_MOLECULES 30
## 1304 REACTOME_CRISTAE_FORMATION 27
## 214 REACTOME_INITIAL_TRIGGERING_OF_COMPLEMENT 21
## 1158 REACTOME_RNA_POLYMERASE_III_CHAIN_ELONGATION 18
## 660 REACTOME_SIGNAL_REGULATORY_PROTEIN_FAMILY_INTERACTIONS 16
## 1082 REACTOME_FORMATION_OF_THE_CORNIFIED_ENVELOPE 123
## 212 REACTOME_COMPLEMENT_CASCADE 54
## 912 REACTOME_DECTIN_2_FAMILY 23
## pMANOVA s.Estill s.Novakovic p.Estill p.Novakovic s.dist
## 145 1.296663e-03 0.71026215 0.3976776 5.029933e-04 5.144144e-02 0.8140146
## 1144 4.144411e-03 0.47281026 0.5886430 2.056827e-02 3.935646e-03 0.7550166
## 657 4.987003e-03 -0.42261523 -0.5557867 2.813103e-02 3.884375e-03 0.6982138
## 1183 4.489133e-03 0.41268206 0.5274865 2.383564e-02 3.870387e-03 0.6697376
## 110 1.112075e-07 0.44740782 0.4682943 1.182208e-05 4.529767e-06 0.6476676
## 646 6.850562e-05 0.47535027 0.3931044 1.623552e-04 1.816773e-03 0.6168378
## 651 5.213322e-03 -0.14824698 -0.5920828 4.169493e-01 1.185440e-03 0.6103599
## 111 2.563631e-08 0.34054463 0.4843115 1.611358e-04 7.996533e-08 0.5920543
## 180 1.629547e-03 -0.34700866 -0.4599334 1.625360e-02 1.445543e-03 0.5761543
## 1060 1.200282e-14 0.28955075 0.4851620 6.636331e-06 4.283235e-14 0.5649972
## 1479 3.949713e-03 -0.36231854 -0.4205387 1.511603e-02 4.801356e-03 0.5550924
## 628 1.528217e-48 0.46289575 0.2647252 3.580471e-45 8.146077e-16 0.5332466
## 394 3.848302e-05 0.38586753 0.3677116 2.539894e-04 4.902626e-04 0.5330155
## 1304 2.039429e-04 -0.29481921 -0.4148275 8.011615e-03 1.905250e-04 0.5089207
## 214 5.372769e-04 0.08043909 0.4869248 5.234224e-01 1.118959e-04 0.4935242
## 1158 4.690039e-03 -0.20394752 -0.4355167 1.341470e-01 1.378569e-03 0.4809048
## 660 4.700215e-03 0.07628916 0.4705094 5.972950e-01 1.119371e-03 0.4766541
## 1082 9.207837e-16 0.22802983 0.4156781 1.259184e-05 1.647017e-15 0.4741159
## 212 4.491245e-07 0.23598077 0.4024949 2.704396e-03 3.102100e-07 0.4665716
## 912 1.516021e-03 0.16882764 0.4297561 1.610665e-01 3.596542e-04 0.4617283
## SD p.adjustMANOVA
## 145 0.22103063 1.717032e-02
## 1144 0.08190611 3.866547e-02
## 657 0.09416648 4.402505e-02
## 1183 0.08117898 4.050086e-02
## 110 0.01476898 8.695368e-06
## 646 0.05815663 1.648963e-03
## 651 0.31383933 4.501974e-02
## 111 0.10165851 2.550148e-06
## 180 0.07984984 1.938924e-02
## 1060 0.13831808 2.463579e-12
## 1479 0.04116787 3.705959e-02
## 628 0.14012773 2.509332e-45
## 394 0.01283819 1.089468e-03
## 1304 0.08485870 3.986597e-03
## 214 0.28742879 8.649105e-03
## 1158 0.16374414 4.171758e-02
## 660 0.27875582 4.171758e-02
## 1082 0.13268737 2.159896e-13
## 212 0.11774326 2.458208e-05
## 912 0.18450425 1.911850e-02det <- mres$detailed_sets$REACTOME_ADRENOCEPTORS
k <- MASS::kde2d(det[,1], det[,2])
X_AXIS = "Estill 450k"
Y_AXIS = "Novakovic EPIC"
XMAX=max(mrank[,1])
XMIN=min(mrank[,1])
YMAX=max(mrank[,2])
YMIN=min(mrank[,2])
filled.contour(k, color.palette = palette, xlim=c(XMIN,XMAX), ylim=c(YMIN,YMAX),
plot.title = {
abline(v = 0, h = 0, lty = 2, lwd = 2, col = "blue")
title(main = "Adrenoreceptors", xlab = X_AXIS,
ylab = Y_AXIS)
})
pdf("fig7d.pdf")
filled.contour(k, color.palette = palette, xlim=c(XMIN,XMAX), ylim=c(YMIN,YMAX),
plot.title = {
abline(v = 0, h = 0, lty = 2, lwd = 2, col = "blue")
title(main = "Adrenoreceptors", xlab = X_AXIS,
ylab = Y_AXIS)
})
dev.off()## pdf
## 2det <- cbind(det, apply(det,1,median))
det <- det[order(det[,3]),]
head(det,20) %>% kbl(caption="Top hypomethylated genes of the selected pathway") %>% kable_paper("hover", full_width = F)| Estill | Novakovic | ||
|---|---|---|---|
| ADRA2C | -15978 | -13696 | -14837.0 |
| ADRB2 | -13214 | -16090 | -14652.0 |
| ADRA2B | -11579 | -16095 | -13837.0 |
| ADRA1B | -13787 | -13520 | -13653.5 |
| ADRA2A | -9457 | -15725 | -12591.0 |
| ADRB1 | -12232 | -11869 | -12050.5 |
| ADRA1D | -11483 | -10528 | -11005.5 |
| ADRB3 | -7071 | -13364 | -10217.5 |
| ADRA1A | -3766 | -6410 | -5088.0 |
Make a html report and some charts.
mitch_report(res=mres,outfile="art_mitchreport.html",overwrite=TRUE)## Note: overwriting existing report## Dataset saved as " /tmp/RtmpXg8hBZ/art_mitchreport.rds ".##
##
## processing file: mitch.Rmd## 1/34
## 2/34 [checklibraries]
## 3/34
## 4/34 [peek]
## 5/34
## 6/34 [metrics]
## 7/34
## 8/34 [scatterplot]## 9/34
## 10/34 [contourplot]## 11/34
## 12/34 [input_geneset_metrics1]
## 13/34
## 14/34 [input_geneset_metrics2]## 15/34
## 16/34 [input_geneset_metrics3]## 17/34
## 18/34 [echart1d]
## 19/34 [echart2d]
## 20/34
## 21/34 [heatmap]## 22/34
## 23/34 [effectsize]
## 24/34
## 25/34 [results_table]
## 26/34
## 27/34 [results_table_complete]
## 28/34
## 29/34 [detailed_geneset_reports1d]
## 30/34
## 31/34 [detailed_geneset_reports2d]## 32/34
## 33/34 [session_info]
## 34/34## output file: /gmea/figs/mitch.knit.md## /usr/local/bin/pandoc +RTS -K512m -RTS /gmea/figs/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpXg8hBZ/mitch_report.html --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/latex-div.lua --embed-resources --standalone --variable bs3=TRUE --section-divs --template /usr/local/lib/R/site-library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --mathjax --variable 'mathjax-url=https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --include-in-header /tmp/RtmpXg8hBZ/rmarkdown-str2ba0e4627db8d0.html##
## Output created: /tmp/RtmpXg8hBZ/mitch_report.html## [1] TRUEmitch_plots(res=mres,outfile="art_mitchcharts.pdf")## pdf
## 2Session information
save.image("art_example.Rdata")
sessionInfo()## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] pkgload_1.3.3
## [2] GGally_2.2.0
## [3] ggplot2_3.4.4
## [4] reshape2_1.4.4
## [5] gtools_3.9.5
## [6] tibble_3.2.1
## [7] dplyr_1.1.4
## [8] echarts4r_0.4.5
## [9] png_0.1-8
## [10] gridExtra_2.3
## [11] missMethyl_1.36.0
## [12] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1
## [13] beeswarm_0.4.0
## [14] kableExtra_1.3.4
## [15] gplots_3.1.3
## [16] mitch_1.15.0
## [17] tictoc_1.2
## [18] HGNChelper_0.8.1
## [19] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [20] IlluminaHumanMethylation450kmanifest_0.4.0
## [21] minfi_1.48.0
## [22] bumphunter_1.44.0
## [23] locfit_1.5-9.8
## [24] iterators_1.0.14
## [25] foreach_1.5.2
## [26] Biostrings_2.70.1
## [27] XVector_0.42.0
## [28] SummarizedExperiment_1.32.0
## [29] Biobase_2.62.0
## [30] MatrixGenerics_1.14.0
## [31] matrixStats_1.2.0
## [32] GenomicRanges_1.54.1
## [33] GenomeInfoDb_1.38.2
## [34] IRanges_2.36.0
## [35] S4Vectors_0.40.2
## [36] BiocGenerics_0.48.1
## [37] eulerr_7.0.0
## [38] limma_3.58.1
##
## loaded via a namespace (and not attached):
## [1] splines_4.3.2 later_1.3.2
## [3] BiocIO_1.12.0 bitops_1.0-7
## [5] filelock_1.0.3 preprocessCore_1.64.0
## [7] XML_3.99-0.16 lifecycle_1.0.4
## [9] lattice_0.22-5 MASS_7.3-60
## [11] base64_2.0.1 scrime_1.3.5
## [13] magrittr_2.0.3 sass_0.4.8
## [15] rmarkdown_2.25 jquerylib_0.1.4
## [17] yaml_2.3.8 httpuv_1.6.13
## [19] doRNG_1.8.6 askpass_1.2.0
## [21] DBI_1.1.3 RColorBrewer_1.1-3
## [23] abind_1.4-5 zlibbioc_1.48.0
## [25] rvest_1.0.3 quadprog_1.5-8
## [27] purrr_1.0.2 RCurl_1.98-1.13
## [29] rappdirs_0.3.3 GenomeInfoDbData_1.2.11
## [31] genefilter_1.84.0 annotate_1.80.0
## [33] svglite_2.1.3 DelayedMatrixStats_1.24.0
## [35] codetools_0.2-19 DelayedArray_0.28.0
## [37] xml2_1.3.6 tidyselect_1.2.0
## [39] farver_2.1.1 beanplot_1.3.1
## [41] BiocFileCache_2.10.1 webshot_0.5.5
## [43] illuminaio_0.44.0 GenomicAlignments_1.38.0
## [45] jsonlite_1.8.8 multtest_2.58.0
## [47] ellipsis_0.3.2 survival_3.5-7
## [49] systemfonts_1.0.5 tools_4.3.2
## [51] progress_1.2.3 Rcpp_1.0.11
## [53] glue_1.6.2 SparseArray_1.2.2
## [55] xfun_0.41 HDF5Array_1.30.0
## [57] withr_2.5.2 fastmap_1.1.1
## [59] rhdf5filters_1.14.1 fansi_1.0.6
## [61] openssl_2.1.1 caTools_1.18.2
## [63] digest_0.6.33 R6_2.5.1
## [65] mime_0.12 colorspace_2.1-0
## [67] biomaRt_2.58.0 RSQLite_2.3.4
## [69] utf8_1.2.4 tidyr_1.3.0
## [71] generics_0.1.3 data.table_1.14.10
## [73] rtracklayer_1.62.0 prettyunits_1.2.0
## [75] httr_1.4.7 htmlwidgets_1.6.4
## [77] S4Arrays_1.2.0 ggstats_0.5.1
## [79] pkgconfig_2.0.3 gtable_0.3.4
## [81] blob_1.2.4 siggenes_1.76.0
## [83] htmltools_0.5.7 scales_1.3.0
## [85] rstudioapi_0.15.0 knitr_1.45
## [87] tzdb_0.4.0 rjson_0.2.21
## [89] nlme_3.1-163 curl_5.2.0
## [91] org.Hs.eg.db_3.18.0 cachem_1.0.8
## [93] rhdf5_2.46.1 stringr_1.5.1
## [95] KernSmooth_2.23-22 AnnotationDbi_1.64.1
## [97] restfulr_0.0.15 GEOquery_2.70.0
## [99] pillar_1.9.0 grid_4.3.2
## [101] reshape_0.8.9 vctrs_0.6.5
## [103] promises_1.2.1 dbplyr_2.4.0
## [105] xtable_1.8-4 evaluate_0.23
## [107] readr_2.1.4 GenomicFeatures_1.54.1
## [109] cli_3.6.2 compiler_4.3.2
## [111] Rsamtools_2.18.0 rlang_1.1.2
## [113] crayon_1.5.2 rngtools_1.5.2
## [115] labeling_0.4.3 nor1mix_1.3-2
## [117] mclust_6.0.1 plyr_1.8.9
## [119] stringi_1.8.3 viridisLite_0.4.2
## [121] BiocParallel_1.36.0 assertthat_0.2.1
## [123] topconfects_1.18.0 munsell_0.5.0
## [125] Matrix_1.6-1.1 hms_1.1.3
## [127] sparseMatrixStats_1.14.0 bit64_4.0.5
## [129] Rhdf5lib_1.24.1 KEGGREST_1.42.0
## [131] statmod_1.5.0 shiny_1.8.0
## [133] highr_0.10 memoise_2.0.1
## [135] bslib_0.6.1 bit_4.0.5