Data loading
suppressPackageStartupMessages({
library(SummarizedExperiment)
library(SEtools)
library(edgeR)
library(DT)
library(pheatmap)
library(plotly)
library(dplyr)
library(sva)
library(ggrepel)
})
source("Functions/EDC_Functions.R")
source("Functions/CriFormatted.R")
load("Data/AllSEcorrected.RData",verbose=T)
## Loading objects:
## SEs
## DEAs
load("Data/HormonalGenes/HormonalGenes.RData",verbose = T)
## Loading objects:
## Thyroid
## Androgen
## Estrogen
## Corticoid
## Progesterone
## PPAR
## Retinoic
load("Data/DEAfetSingleCompunds.RData", verbose=TRUE)
## Loading objects:
## res.fetT3
## res.fetBPA
load("Data/DEAorgSingleCompunds.RData", verbose=TRUE)
## Loading objects:
## res.orgT3
## res.orgBPA
MixN vs single compounds Fetal
Fetal MixN DEGs are specific for MixN treated samples
fetal <- SEs$chronic.fetal[,which(SEs$chronic.fetal$EXPO=="CNT" | SEs$chronic.fetal$EXPO=="1X" | SEs$chronic.fetal$EXPO=="1000X" | SEs$chronic.fetal$EXPO=="T3"| SEs$chronic.fetal$EXPO=="BPA0.04X")]
degs.fe <- row.names(DEAs$chronic.fetal)[which((abs(DEAs$chronic.fetal$logFC.EXPO1X)>0.5 | abs(DEAs$chronic.fetal$logFC.EXPO1000X)>0.5) & DEAs$chronic.fetal$FDR<=0.05 & DEAs$chronic.fetal$logCPM>0)]
sehm(fetal, assayName = "corrected", degs.fe, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="MixNDEGs")
Fetal T3 treated samples respond to T3 known targets and the DEGs are specific for T3 samples
sehm(fetal[,which(SEs$chronic.fetal$EXPO=="CNT" | SEs$chronic.fetal$EXPO=="T3")], assayName = "corrected", Thyroid, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="T3 target genes")
degs.fetT3 <- row.names(res.fetT3)[which(abs(res.fetT3$logFC)> 0.5 & res.fetT3$FDR<=0.05 & res.fetT3$logCPM > 0)]
sehm(fetal, assayName = "corrected", degs.fetT3, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="T3DEGs")
T3 <- res.fetT3
T3$genes <- row.names(T3)
MixN <- DEAs$chronic.fetal
MixN$genes <- row.names(MixN)
MixN$logFC <- (MixN$logFC.EXPO1X + MixN$logFC.EXPO1000X)/2
Comp <- compareResultsFCNew(T3, MixN, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, logCPMth=0, FCceil=5,
title='Fetal: Comparing T3 to MixN (mean FC)', geneLabel=TRUE, topLab=-30)
Comp$Scatter
#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterFetal_MixN_T3.pdf', width=10, height=10)
#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterFetal_MixN_T3.png', width=10, height=10)
# Check overlap between MixN and Thyroid
table(Comp$AllRes$Status)
##
## NotSignificant SignificantA SignificantA&B SignificantB
## 12488 374 17 44
## Fisher Test Significant
SigBoth <- dplyr::filter(Comp$Significant, Comp$Significant$Status=='SignificantA&B')
Contingency <- matrix(c(dim(dplyr::filter(SigBoth, log2FCA < 0 & log2FCB < 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA >= 0 & log2FCB < 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA < 0 & log2FCB >= 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA >= 0 & log2FCB >= 0 ))[1]),
nrow=2, dimnames=list(c('DownA', 'UpA'), c('DownB', 'UpB')))
Contingency
## DownB UpB
## DownA 11 3
## UpA 1 2
Fetal BPA treated samples show specific DEGs
degs.fetBPA <- row.names(res.fetBPA)[which(abs(res.fetBPA$logFC)> 0.5 & res.fetBPA$FDR<=0.05 & res.fetBPA$logCPM > 0)]
sehm(fetal, assayName = "corrected", degs.fetBPA, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="BPADEGs")
BPA <- res.fetBPA
BPA$genes <- row.names(BPA)
MixN <- DEAs$chronic.fetal
MixN$genes <- row.names(MixN)
MixN$logFC <- MixN$logFC.EXPO1X
Comp <- compareResultsFCNew(BPA, MixN, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, logCPMth=0, FCceil=4,
title='Fetal: Comparing BPA to MixN (1X FC)', geneLabel=TRUE, topLab=-30)
Comp$Scatter
#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterFetal_MixN_BPA.pdf', width=10, height=10)
#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterFetal_MixN_BPA.png', width=10, height=10)
MixN vs single compounds Organoids
Organoids MixN DEGs are specific for MixN treated samples
Organoids <- SEs$chronic.org[,which(SEs$chronic.org$EXPO=="CNT"|SEs$chronic.org$EXPO=="1X" | SEs$chronic.org$EXPO=="1000X" | SEs$chronic.org$EXPO=="T3"| SEs$chronic.org$EXPO=="BPA0.04X")]
degs.org <- row.names(DEAs$chronic.org)[which((abs(DEAs$chronic.org$logFC.EXPO1X)>0.5 | abs(DEAs$chronic.org$logFC.EXPO1000X)>0.5) & DEAs$chronic.org$FDR<=0.05 & DEAs$chronic.org$logCPM>0)]
sehm(Organoids, assayName = "corrected", degs.org, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="MixNDEGs")
Organoids T3 treated samples respond to T3 known targets and the DEGs are specific for T3 samples
sehm(Organoids[,which(SEs$chronic.org$EXPO=="CNT" | SEs$chronic.org$EXPO=="T3")], assayName = "corrected", Thyroid, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="T3 target genes")
degs.orgT3 <- row.names(res.orgT3)[which(abs(res.orgT3$logFC)> 0.5 & res.orgT3$FDR<=0.05 & res.orgT3$logCPM > 0)]
sehm(Organoids, assayName = "corrected", degs.orgT3, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="T3DEGs")
T3 <- res.orgT3
T3$genes <- row.names(T3)
MixN <- DEAs$chronic.org
MixN$genes <- row.names(MixN)
MixN$logFC <- (MixN$logFC.EXPO1X + MixN$logFC.EXPO1000X)/2
Comp <- compareResultsFCNew(T3, MixN, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, logCPMth=0, FCceil=6,
title='Comparing T3 to MixN (mean FC)', geneLabel=TRUE, topLab=-40)
Comp$Scatter
#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterOrg_MixN_T3.pdf', width=10, height=10)
#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterOrg_MixN_T3.png', width=10, height=10)
# Check overlap between MixN and Thyroid
table(Comp$AllRes$Status)
##
## NotSignificant SignificantA SignificantA&B SignificantB
## 12331 2212 60 286
## Contingency table on shared genes
SigBoth <- dplyr::filter(Comp$Significant, Comp$Significant$Status=='SignificantA&B')
ContingencyII <- matrix(c(dim(dplyr::filter(SigBoth, log2FCA < 0 & log2FCB < 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA >= 0 & log2FCB < 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA < 0 & log2FCB >= 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA >= 0 & log2FCB >= 0 ))[1]),
nrow=2, dimnames=list(c('DownA', 'UpA'), c('DownB', 'UpB')))
ContingencyII
## DownB UpB
## DownA 14 36
## UpA 10 0
prop.table(ContingencyII)
## DownB UpB
## DownA 0.2333333 0.6
## UpA 0.1666667 0.0
dplyr::filter(SigBoth, log2FCA < 0 & log2FCB > 0 ) %>% dplyr::pull(genes)
## [1] "SP5" "COL9A2" "KRT8" "DDIT4L" "IGSF21" "NEUROG1"
## [7] "ID1" "HIST2H3A" "HIST2H3C" "CDC20B" "METRN" "ESM1"
## [13] "SPATA18" "RPE65" "SLC35F2" "TRIL" "HIST2H3D" "TRAPPC5"
## [19] "TMEM121" "S1PR3" "SALL3" "FEZF2" "PTHLH" "NDUFS7"
## [25] "TPH1" "SOD3" "ZNF503" "TGFB2" "PKMYT1" "DMRT3"
## [31] "COL9A3" "C4orf48" "CCDC85B" "APOE" "COL18A1" "C19orf24"
Organoids BPA treated samples show specific DEGs
degs.orgBPA <- row.names(res.orgBPA)[which(abs(res.orgBPA$logFC)> 0.5 & res.orgBPA$FDR<=0.05 & res.orgBPA$logCPM > 0)]
sehm(Organoids, assayName = "corrected", degs.orgBPA, do.scale = T, show_rownames = F, anno_columns = c("EXPO2","MixtureBatch"), main="BPADEGs")
BPA <- res.orgBPA
BPA$genes <- row.names(BPA)
MixN <- DEAs$chronic.org
MixN$genes <- row.names(MixN)
MixN$logFC <- MixN$logFC.EXPO1X
#Comp <- compareResultsFC(BPA, MixN, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, FCceil=8,logCPMth = 0, title='Comparing BPA to MixN 1X', geneLabel=TRUE)
#Comp$Scatter
Comp <- compareResultsFCNew(BPA, MixN, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, logCPMth=0, FCceil=4,
title='Organoids: Comparing BPA to MixN (1X FC)', geneLabel=TRUE, topLab=-30)
Comp$Scatter
#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterOrg_MixN_BPA.pdf', width=10, height=10)
#ggsave(Comp$Scatter, filename='Data/SingleCompounds/ScatterOrg_MixN_BPA.png', width=10, height=10)
# Check overlap between MixN and BPA
table(Comp$AllRes$Status)
##
## NotSignificant SignificantA SignificantA&B SignificantB
## 13770 821 50 248
## Contingency table on shared genes
SigBoth <- dplyr::filter(Comp$Significant, Comp$Significant$Status=='SignificantA&B')
ContingencyII <- matrix(c(dim(dplyr::filter(SigBoth, log2FCA < 0 & log2FCB < 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA >= 0 & log2FCB < 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA < 0 & log2FCB >= 0 ))[1],
dim(dplyr::filter(SigBoth, log2FCA >= 0 & log2FCB >= 0 ))[1]),
nrow=2, dimnames=list(c('DownA', 'UpA'), c('DownB', 'UpB')))
ContingencyII
## DownB UpB
## DownA 32 0
## UpA 5 13
Data praparation
For details on data filtering, normalization, batch correction and differential expression analysis, see here and here
Authors
Nicolò Caporale: nicolo.caporale@ieo.it, nicolo.caporale@unimi.it
Cristina Cheroni: cristina.cheroni@ieo.it
Pierre-Luc Germain: pierre-luc.germain@ieo.it
Giuseppe Testa: giuseppe.testa@ieo.it
Lab: http://www.testalab.eu/
‘Date: December 01, 2021’
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ggrepel_0.9.1 sva_3.38.0
## [3] BiocParallel_1.24.1 genefilter_1.72.1
## [5] mgcv_1.8-34 nlme_3.1-152
## [7] dplyr_1.0.5 plotly_4.9.3
## [9] ggplot2_3.3.3 pheatmap_1.0.12
## [11] DT_0.17 edgeR_3.32.1
## [13] limma_3.46.0 SEtools_1.4.0
## [15] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [17] GenomicRanges_1.42.0 GenomeInfoDb_1.26.2
## [19] IRanges_2.24.1 S4Vectors_0.28.1
## [21] BiocGenerics_0.36.0 MatrixGenerics_1.2.1
## [23] matrixStats_0.58.0
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.15 colorspace_2.0-0 rjson_0.2.20
## [4] ellipsis_0.3.1 circlize_0.4.12 XVector_0.30.0
## [7] GlobalOptions_0.1.2 clue_0.3-58 farver_2.1.0
## [10] bit64_4.0.5 AnnotationDbi_1.52.0 fansi_0.4.2
## [13] codetools_0.2-18 splines_4.0.3 cachem_1.0.4
## [16] knitr_1.31 jsonlite_1.7.2 Cairo_1.5-12.2
## [19] annotate_1.68.0 cluster_2.1.1 png_0.1-7
## [22] compiler_4.0.3 httr_1.4.2 assertthat_0.2.1
## [25] Matrix_1.3-2 fastmap_1.1.0 lazyeval_0.2.2
## [28] htmltools_0.5.1.1 tools_4.0.3 gtable_0.3.0
## [31] glue_1.4.2 GenomeInfoDbData_1.2.4 V8_3.4.0
## [34] Rcpp_1.0.6 jquerylib_0.1.3 vctrs_0.3.7
## [37] iterators_1.0.13 xfun_0.21 stringr_1.4.0
## [40] openxlsx_4.2.3 lifecycle_1.0.0 XML_3.99-0.5
## [43] zlibbioc_1.36.0 scales_1.1.1 TSP_1.1-10
## [46] RColorBrewer_1.1-2 ComplexHeatmap_2.6.2 yaml_2.2.1
## [49] curl_4.3 memoise_2.0.0 sass_0.3.1
## [52] stringi_1.5.3 RSQLite_2.2.3 highr_0.8
## [55] randomcoloR_1.1.0.1 foreach_1.5.1 seriation_1.2-9
## [58] zip_2.1.1 shape_1.4.5 rlang_0.4.10
## [61] pkgconfig_2.0.3 bitops_1.0-6 evaluate_0.14
## [64] lattice_0.20-41 purrr_0.3.4 labeling_0.4.2
## [67] htmlwidgets_1.5.3 bit_4.0.4 tidyselect_1.1.0
## [70] magrittr_2.0.1 R6_2.5.0 generics_0.1.0
## [73] DelayedArray_0.16.1 DBI_1.1.1 pillar_1.5.1
## [76] withr_2.4.1 survival_3.2-7 RCurl_1.98-1.2
## [79] tibble_3.1.0 crayon_1.4.1 utf8_1.2.1
## [82] rmarkdown_2.7 GetoptLong_1.0.5 locfit_1.5-9.4
## [85] grid_4.0.3 data.table_1.14.0 blob_1.2.1
## [88] digest_0.6.27 xtable_1.8-4 tidyr_1.1.3
## [91] munsell_0.5.0 registry_0.5-1 viridisLite_0.3.0
## [94] bslib_0.2.4