suppressPackageStartupMessages({
library(SummarizedExperiment)
library(edgeR)
library(DT)
library(pheatmap)
library(plotly)
library(dplyr)
library(sva)
})
load("../Data/geneLengths.RData")
source("../Functions/EDC_Functions.R")
load("../Data/TotalSumExp.RData",verbose=T)
## Loading objects:
## Total
Chronic expo
# subset to chronic data
se <- Total[,which(Total$EXPOduration %in% c("Chronic 15 days","Chronic50days"))]
# refactor EXPO to have DMSO as the base factor, and the other levels in a decent order
se$EXPO <- factor(as.character(se$EXPO),levels=c("DMSO","CNT","1X","1000X","BPA0.04X","T3","VITC"))
se$EXPO2 <- se$EXPO
# we will consider CNT and DMSO together; comment out following to use DMSO as baseline
levels(se$EXPO) <- c("CNT","CNT","1X","1000X","BPA0.04X","T3","VITC")
# we treat the batch as factor and use 2018 as the base level because it has the most samples
se$MixtureBatch <- relevel(as.factor(se$MixtureBatch),"2018")
se <- se[,order(se$type, se$EXPO, se$MixtureBatch)]
# we filter and normalize the two datasets separately:
org <- filterGenes(se[,which(se$EXPOduration=="Chronic50days")])
## 9986 rows out of 26300 were removed.
dds.org <- calcNormFactors(DGEList(assay(org)))
assays(org)$logcpm <- log(cpm(dds.org)+0.1)
fe <- filterGenes(se[,which(se$EXPOduration=="Chronic 15 days")])
## 12283 rows out of 26300 were removed.
dds.fe <- calcNormFactors(DGEList(assay(fe)))
assays(fe)$logcpm <- log(cpm(dds.fe)+0.1)
Organoids chronic
# performs the SVA analysis
sv <- svacor(org, ~MixtureBatch+EXPO, ~MixtureBatch,regressOutNull = T)
## Using variance-stabilizing transformation
## converting counts to integer mode
## Number of significant surrogate variables is: 6
## Iteration (out of 5 ):1 2 3 4 5
# we put the corrected data in the total SE object
assays(org)$corrected <- sv$cor
# we include the SVs in the object's metadata
colData(org) <- cbind(colData(org),sv$sv)
#DEA T3
fit <- glmQLFit(estimateDisp(dds.org,sv$mm),sv$mm)
res.orgT3 <- as.data.frame(topTags(glmQLFTest(fit,c("EXPOT3")),n=Inf))
degs.orgT3 <- row.names(res.orgT3)[which(abs(res.orgT3$logFC)> 0.5 & res.orgT3$FDR<=0.05 & res.orgT3$logCPM>0)]
res.orgBPA <- as.data.frame(topTags(glmQLFTest(fit,c("EXPOBPA0.04X")),n=Inf))
degs.orgBPA <- row.names(res.orgBPA)[which(abs(res.orgBPA$logFC)> 0.5 & res.orgBPA$FDR<=0.05 & res.orgBPA$logCPM>0)]
#save(res.orgT3,res.orgBPA, file = "../Data/DEAorgSingleCompunds.RData")
Fetal chronic
# performs the SVA analysis
sv <- svacor(fe, ~MixtureBatch+EXPO, ~MixtureBatch,regressOutNull = T)
## Using variance-stabilizing transformation
## converting counts to integer mode
## Number of significant surrogate variables is: 6
## Iteration (out of 5 ):1 2 3 4 5
# we put the corrected data in the total SE object
assays(fe)$corrected <- sv$cor
# we include the SVs in the object's metadata
colData(fe) <- cbind(colData(fe),sv$sv)
#DEA T3
fit <- glmQLFit(estimateDisp(dds.fe,sv$mm),sv$mm)
res.fetT3 <- as.data.frame(topTags(glmQLFTest(fit,c("EXPOT3")),n=Inf))
degs.fetT3 <- row.names(res.fetT3)[which(abs(res.fetT3$logFC)> 0.5 & res.fetT3$FDR<=0.05 & res.fetT3$logCPM>0)]
res.fetBPA <- as.data.frame(topTags(glmQLFTest(fit,c("EXPOBPA0.04X")),n=Inf))
degs.fetBPA <- row.names(res.fetBPA)[which(abs(res.fetBPA$logFC)> 0.5 & res.fetBPA$FDR<=0.05 & res.fetBPA$logCPM>0)]
#save(res.fetT3,res.fetBPA, file = "../Data/DEAfetSingleCompunds.RData")
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 13, 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] DESeq2_1.30.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 SummarizedExperiment_1.20.0
## [15] Biobase_2.50.0 GenomicRanges_1.42.0
## [17] GenomeInfoDb_1.26.2 IRanges_2.24.1
## [19] S4Vectors_0.28.1 BiocGenerics_0.36.0
## [21] MatrixGenerics_1.2.1 matrixStats_0.58.0
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.2 sass_0.3.1 tidyr_1.1.3
## [4] bit64_4.0.5 jsonlite_1.7.2 viridisLite_0.3.0
## [7] splines_4.0.3 bslib_0.2.4 assertthat_0.2.1
## [10] blob_1.2.1 GenomeInfoDbData_1.2.4 yaml_2.2.1
## [13] RSQLite_2.2.3 pillar_1.5.1 lattice_0.20-41
## [16] glue_1.4.2 digest_0.6.27 RColorBrewer_1.1-2
## [19] XVector_0.30.0 colorspace_2.0-0 htmltools_0.5.1.1
## [22] Matrix_1.3-2 XML_3.99-0.5 pkgconfig_2.0.3
## [25] zlibbioc_1.36.0 purrr_0.3.4 xtable_1.8-4
## [28] scales_1.1.1 tibble_3.1.0 annotate_1.68.0
## [31] generics_0.1.0 ellipsis_0.3.1 cachem_1.0.4
## [34] withr_2.4.1 lazyeval_0.2.2 survival_3.2-7
## [37] magrittr_2.0.1 crayon_1.4.1 memoise_2.0.0
## [40] evaluate_0.14 fansi_0.4.2 tools_4.0.3
## [43] data.table_1.14.0 lifecycle_1.0.0 stringr_1.4.0
## [46] munsell_0.5.0 locfit_1.5-9.4 DelayedArray_0.16.1
## [49] AnnotationDbi_1.52.0 compiler_4.0.3 jquerylib_0.1.3
## [52] rlang_0.4.10 grid_4.0.3 RCurl_1.98-1.2
## [55] htmlwidgets_1.5.3 bitops_1.0-6 rmarkdown_2.7
## [58] gtable_0.3.0 DBI_1.1.1 R6_2.5.0
## [61] knitr_1.31 fastmap_1.1.0 bit_4.0.4
## [64] utf8_1.2.1 stringi_1.5.3 Rcpp_1.0.6
## [67] geneplotter_1.68.0 vctrs_0.3.7 tidyselect_1.1.0
## [70] xfun_0.21