Data loading
suppressPackageStartupMessages({
library(SummarizedExperiment)
library(SEtools)
library(edgeR)
library(DT)
library(pheatmap)
library(plotly)
library(dplyr)
library(sva)
library(reshape2)
library(overlapper)
library(ggrepel)
})
source("Functions/EDC_Functions.R")
source("Functions/CriFormatted.R")
load("Data/TotalSumExp.RData")
load("Data/geneLengths.RData")
load("Data/AllSEcorrected.RData",verbose=T)
## Loading objects:
## SEs
## DEAs
load("Data/brainspan.cortex.byAge.RData", verbose = T)
## Loading objects:
## byAge
load("Data/iPSC.RData", verbose = T)
## Loading objects:
## e
Brainspan correlation
Organoids <- assays(SEs$chronic.org[,which(SEs$chronic.org$EXPO=="CNT")])$counts
Fetal2D <- assays(SEs$chronic.fetal[,which(SEs$chronic.fetal$EXPO=="CNT")])$counts
iPSC <- e[,grep("3391B",colnames(e))]
Fetal <- byAge
i <- intersect(rownames(Organoids),intersect(rownames(iPSC),rownames(Fetal2D)))
tot <- cbind(Organoids[i,],iPSC[i,],Fetal2D[i,])
tot <- filterGenes(tot)
en <- donorm(tot)
#put the rows in the same order
i <- intersect(row.names(Fetal),row.names(en))
en <- en[i,]
Fetal <- Fetal[i,]
#create the matrix of pearson correlation coefficients and plot
m <- matrix(0,nrow=ncol(Fetal),ncol=ncol(en))
for(i in 1:ncol(en)){ m[,i] <- apply(Fetal,2,FUN=function(x){ cor(log(en[,i]+1),log(x+1)) }) }
row.names(m) <- colnames(Fetal)
colnames(m) <- colnames(en)
m2 <- m[grep("pcw",row.names(m)),]
colnames(m2) <- c(rep("Organoids",8), rep("iPSC",5), rep("FetalProgenitors",8))
byheatmap(m2,cluster_row=F)
DEGs heatmaps
org.chronic <- 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)]
MixN_org_chronic <- SEs$chronic.org[,which(SEs$chronic.org$EXPO=="CNT"|SEs$chronic.org$EXPO=="1X" | SEs$chronic.org$EXPO=="1000X")]
sehm(MixN_org_chronic, assayName = "corrected", org.chronic, do.scale = T, show_rownames = F, anno_columns = c("EXPO2"), main="Chronic org DEGs")
1X vs 1000X
FDRTh <- 0.05
LogFCTh <- 0.5
MixN1X <- DEAs$chronic.org
MixN1X$genes <- row.names(MixN1X)
MixN1X$logFC <- (MixN1X$logFC.EXPO1X)
MixN1000X <- DEAs$chronic.org
MixN1000X$genes <- row.names(MixN1000X)
MixN1000X$logFC <- (MixN1X$logFC.EXPO1000X)
#Comp <- compareResultsFC(MixN1X, MixN1000X, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, logCPMth = 0, FCceil=8, title='Comparing MixN 1X to 1000X', geneLabel=TRUE)
#Comp$Scatter
Comp <- compareResultsFCNew(MixN1X, MixN1000X, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, logCPMth=0, FCceil=2.5,
title='Comparing MixN 1X to 1000X', geneLabel=TRUE, topLab=-16)
## Warning in cor.test.default(AllRes$logFC_A, AllRes$logFC_B, method = corMethod):
## Cannot compute exact p-value with ties
Comp$Scatter
ggsave(Comp$Scatter, filename='Data/ChronicOrganoids/Scatter.pdf', width=10, height=10)
ggsave(Comp$Scatter, filename='Data/ChronicOrganoids/Scatter.png', width=10, height=10)
Neurodevelopmental Disorder (NDD) Genes
load("Data/ASD.RData", verbose = T)
## Loading objects:
## ASD
## SFARI
## SFARIgenes
## NeuropsychiatricDiseases
## PsychencodeNDD
org.chronic3 <- row.names(DEAs$chronic.org)[which(DEAs$chronic.org$FDR<=0.05)]
controlNeg <- row.names(DEAs$chronic.org)[-which(DEAs$chronic.org$FDR<=0.05)]
chronic.orgUp <- row.names(DEAs$chronic.org)[which(((DEAs$chronic.org$logFC.EXPO1X)>0 & (DEAs$chronic.org$logFC.EXPO1000X)>0) & DEAs$chronic.org$FDR<=0.05)]
chronic.orgDown <- row.names(DEAs$chronic.org)[which(((DEAs$chronic.org$logFC.EXPO1X)<0 & (DEAs$chronic.org$logFC.EXPO1000X)<0) & DEAs$chronic.org$FDR<=0.05)]
#MixN <- list(organoidsDEGs=org.chronic3, organoidsDEGsUp=chronic.orgUp, organoidsDEGsDown=chronic.orgDown,NonAffectedGenes=controlNeg)
MixN <- list(organoidsDEGs=org.chronic3, NonAffectedGenes=controlNeg)
m <- overlapper::multintersect(ll = MixN, ll2 = ASD, universe = row.names(DEAs$chronic.org),two.tailed = F)
org.chronicLong<-org.chronic3
save(org.chronic,org.chronicLong,file = "Data/DEGsOrganoidsChronic.RData")
dotplot.multintersect(m, sizeRange = c(0,15), th=0.05)
ggsave(dotplot.multintersect(m, th=0.05), filename='Data/ChronicOrganoids/DotplotOrgNDD.pdf', width=10, height=6)
P values of the overlaps
m$prob <- round(m$prob, digits=3)
DT::datatable(m$prob)
SFARI genes
m <- overlapper::multintersect(ll = MixN, ll2 = SFARIgenes, universe = row.names(DEAs$chronic.org),two.tailed = F)
dotplot.multintersect(m, th=0.05)
ggsave(dotplot.multintersect(m, th=0.05), filename='Data/ChronicOrganoids/DotplotOrgSFARI.pdf', width=10, height=6)
NPD genes
m <- overlapper::multintersect(ll = MixN, ll2 = NeuropsychiatricDiseases, universe = row.names(DEAs$chronic.org),two.tailed = F)
dotplot.multintersect(m, th=0.05)
ggsave(dotplot.multintersect(m, th=0.05), filename='Data/ChronicOrganoids/DotplotOrgNPD.pdf', width=10, height=6)
NDD psychencode
m <- overlapper::multintersect(ll = MixN, ll2 = PsychencodeNDD, universe = row.names(DEAs$chronic.org),two.tailed = F)
dotplot.multintersect(m, th=0.05)
ggsave(dotplot.multintersect(m, th=0.05), filename='Data/ChronicFetal/DotplotOrgNDD.pdf', width=10, height=6)
Gene expression dysregulation of the most important NDD genes
highest scores
geneStripPairEDCMix(SE = MixN_org_chronic,GeneSet =intersect(org.chronic, union(union(SFARIgenes$score1,SFARIgenes$score2),SFARIgenes$score3)), printExp = FALSE, SampleColors = "Default")
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.ymin` is deprecated. Use `fun.min` instead.
## Warning: `fun.ymax` is deprecated. Use `fun.max` instead.
all scores
geneStripPairEDCMix(SE = MixN_org_chronic,GeneSet =intersect(org.chronic, unlist(SFARIgenes)), printExp = FALSE, SampleColors = "Default")
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.ymin` is deprecated. Use `fun.min` instead.
## Warning: `fun.ymax` is deprecated. Use `fun.max` instead.
PVals <- DEAs$chronic.org[intersect(org.chronic, unlist(SFARIgenes)),]
PVals <- PVals[complete.cases(PVals), ]
PVals
## logFC.EXPO1X logFC.EXPO1000X logCPM F PValue
## SCN9A 0.16332099 -0.73394503 3.6238927 28.679243 4.486322e-08
## NEFL 0.17224825 -0.67496037 7.2905518 26.587394 1.004741e-07
## CACNA1E -0.22395439 -0.68438361 4.9734221 21.775175 7.506001e-07
## SCN2A -0.23500980 -0.75882170 4.4730437 21.537869 8.340935e-07
## PCDH9 -0.35329648 -0.51191575 6.2629341 18.479058 3.451823e-06
## CDH9 -0.63189859 -1.50276593 0.4498228 17.344314 6.030875e-06
## SPARCL1 0.28517770 -0.59496470 5.7544922 15.403606 1.636156e-05
## TH 0.31734183 -0.56675217 2.3380022 14.953642 2.079647e-05
## NEGR1 -0.04352262 -0.54219055 4.4452160 14.906921 2.132497e-05
## SLC6A1 -0.21067534 -0.57822638 3.0876260 13.077681 5.868259e-05
## NRXN3 -0.50939333 0.21190259 5.1255620 12.755676 7.057116e-05
## KCNMA1 -0.65819131 -0.33457106 2.9964461 12.306480 9.158739e-05
## GRM8 0.59891403 0.30232488 5.4549395 11.619416 1.375170e-04
## CACNB2 -0.25652132 -0.51771213 3.2092017 11.613580 1.379983e-04
## DPYD -0.41965092 -0.59127083 3.0251296 11.342770 1.624122e-04
## LMX1B 0.22470186 -0.62563738 4.3620879 11.132867 1.844648e-04
## KDM6B 0.28440194 0.54084969 5.1366501 10.480879 2.756171e-04
## WNT1 0.42357343 -0.73059040 1.7103686 10.118513 3.459453e-04
## FABP3 -0.54794627 -0.50839479 3.3966961 9.978136 3.780888e-04
## CNTN5 0.05897304 -0.60625916 1.6670097 8.563565 9.503601e-04
## DRD2 0.17067356 -0.65784667 1.0232263 8.552201 9.576166e-04
## KIRREL3 -0.06077077 -0.51471524 2.7742917 8.078676 1.318723e-03
## PCDH11X -0.32805473 -0.64062334 4.6909851 7.525175 1.931552e-03
## GRM7 0.05663132 -0.61625257 2.1733713 7.382859 2.133643e-03
## SLC4A10 0.10963873 0.57241339 5.2768430 7.242564 2.354874e-03
## ACE 0.50518870 -0.00501186 1.4608377 7.221730 2.389741e-03
## FABP7 -0.58477173 -0.21714537 4.9423324 6.660712 3.567096e-03
## THBS1 -0.55858202 -0.12754071 3.7180289 6.643538 3.611637e-03
## FEZF2 0.33938486 0.89684384 2.2713893 6.558978 3.839691e-03
## ADAMTS18 0.30426384 0.60810219 3.1172417 6.255539 4.791842e-03
## SAMD11 0.61769439 0.16814633 2.5613572 5.718011 7.145064e-03
## CACNA1I -0.56456711 -0.82961573 1.3985198 5.689791 7.298381e-03
## FOXP2 -0.57048773 0.03921951 5.9517196 5.616600 7.712461e-03
## FDR
## SCN9A 0.0001156523
## NEFL 0.0001170810
## CACNA1E 0.0003182873
## SCN2A 0.0003239857
## PCDH9 0.0006120982
## CDH9 0.0008920245
## SPARCL1 0.0014428240
## TH 0.0016448800
## NEGR1 0.0016487940
## SLC6A1 0.0029098719
## NRXN3 0.0031645244
## KCNMA1 0.0036541645
## GRM8 0.0047692248
## CACNB2 0.0047697124
## DPYD 0.0052581297
## LMX1B 0.0057095415
## KDM6B 0.0070921402
## WNT1 0.0081205050
## FABP3 0.0084843741
## CNTN5 0.0146529435
## DRD2 0.0147116784
## KIRREL3 0.0175908790
## PCDH11X 0.0222464372
## GRM7 0.0234873475
## SLC4A10 0.0246898511
## ACE 0.0249591803
## FABP7 0.0313881338
## THBS1 0.0316627215
## FEZF2 0.0327906209
## ADAMTS18 0.0371020919
## SAMD11 0.0473454810
## CACNA1I 0.0477600401
## FOXP2 0.0495748998
- Genes with conventional PVal > 0.05:
- Genes with conventional PVal < 0.05:
- Genes with conventional FDR < 0.05: SCN9A, NEFL, CACNA1E, SCN2A, PCDH9, CDH9, SPARCL1, TH, NEGR1, SLC6A1, NRXN3, KCNMA1, GRM8, CACNB2, DPYD, LMX1B, KDM6B, WNT1, FABP3, CNTN5, DRD2, KIRREL3, PCDH11X, GRM7, SLC4A10, ACE, FABP7, THBS1, FEZF2, ADAMTS18, SAMD11, CACNA1I, FOXP2
Gene expression dysregulation of the fetal NDD-DEGs in organoids dataset
load("Data/DEGsFetalChronic.RData", verbose = T)
## Loading objects:
## fetal.chronic
## fetal.chronicLong
geneStripPairEDCMix(SE = MixN_org_chronic,GeneSet = intersect(fetal.chronic, unlist(SFARIgenes)), printExp = FALSE,SampleColors = "Default")
## [1] "Expression values are not available for the following genes: CD38"
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.ymin` is deprecated. Use `fun.min` instead.
## Warning: `fun.ymax` is deprecated. Use `fun.max` instead.
PVals <- DEAs$chronic.org[intersect(fetal.chronic, unlist(SFARIgenes)),]
PVals <- PVals[complete.cases(PVals), ]
PVals
## logFC.EXPO1X logFC.EXPO1000X logCPM F PValue FDR
## RANBP17 -0.1626618 -0.2152368 5.277389 4.470447 0.018742025 0.08479478
## UPF3B 0.0328277 -0.2586830 3.815079 2.967897 0.064606056 0.17889598
## KCNJ2 0.4778245 0.1451159 6.134523 5.280253 0.009961453 0.05816433
- Genes with conventional PVal > 0.05: UPF3B
- Genes with conventional PVal < 0.05: RANBP17, KCNJ2
- Genes with conventional FDR < 0.05:
Gene expression dysregulation of the genes found in the acute expo
load("Data/DEGsFetalAcute.RData", verbose = T)
## Loading objects:
## fet.acute
geneStripPairEDCMix(SE = MixN_org_chronic,GeneSet = intersect(fet.acute,org.chronic), printExp = FALSE,SampleColors = "Default")
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.ymin` is deprecated. Use `fun.min` instead.
## Warning: `fun.ymax` is deprecated. Use `fun.max` instead.
PVals <- DEAs$chronic.org[intersect(fet.acute,org.chronic),]
PVals <- PVals[complete.cases(PVals), ]
PVals
## logFC.EXPO1X logFC.EXPO1000X logCPM F PValue FDR
## NEFL 0.1722482 -0.6749604 7.290552 26.58739 1.004741e-07 0.000117081
- Genes with conventional PVal > 0.05:
- Genes with conventional PVal < 0.05:
- Genes with conventional FDR < 0.05: NEFL
geneStripPairEDCMix(SE = MixN_org_chronic,GeneSet = c("CLSTN2", "EPHB2"), printExp = FALSE,SampleColors = "Default")
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.ymin` is deprecated. Use `fun.min` instead.
## Warning: `fun.ymax` is deprecated. Use `fun.max` instead.
PVals <- DEAs$chronic.org[c("CLSTN2", "EPHB2"),]
PVals <- PVals[complete.cases(PVals), ]
PVals
## logFC.EXPO1X logFC.EXPO1000X logCPM F PValue FDR
## CLSTN2 0.214876852 0.09619874 5.239390 1.447981 0.24891393 0.4130782
## EPHB2 0.006988673 0.16640719 7.820478 3.802717 0.03213182 0.1164109
- Genes with conventional PVal > 0.05: CLSTN2
- Genes with conventional PVal < 0.05: EPHB2
- Genes with conventional FDR < 0.05:
Comparison across systems
more stringent DEGs
load("Data/DEGsFetalChronic.RData", verbose = T)
## Loading objects:
## fetal.chronic
## fetal.chronicLong
load("Data/DEGsOrganoidsChronic.RData", verbose = T)
## Loading objects:
## org.chronic
## org.chronicLong
MixN <- list(OrganoidsDEGs=org.chronic, FetalAcuteDEGs=fet.acute, FetalChronicDEGs=fetal.chronic)
col <- c("low"="grey", "high"="yellow")
m<-multintersect(ll = MixN, universe = row.names(filterGenes(Total)),two.tailed = F)
dotplot.multintersect(m, th=0.05,colors = col)
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
less stringent DEGs
MixN <- list(OrganoidsDEGs=org.chronicLong, FetalAcuteDEGs=fet.acute, FetalChronicDEGs=fetal.chronicLong)
m<-multintersect(ll = MixN, universe = row.names(filterGenes(Total)),two.tailed = F)
dotplot.multintersect(m, th=0.05)
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
Data praparation
For details on data filtering, normalization, batch correction and differential expression analysis, see here
For details on the preparation of Neurodevelopmental disorder genes, see 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 overlapper_0.99.1
## [3] reshape2_1.4.4 sva_3.38.0
## [5] BiocParallel_1.24.1 genefilter_1.72.1
## [7] mgcv_1.8-34 nlme_3.1-152
## [9] dplyr_1.0.5 plotly_4.9.3
## [11] ggplot2_3.3.3 pheatmap_1.0.12
## [13] DT_0.17 edgeR_3.32.1
## [15] limma_3.46.0 SEtools_1.4.0
## [17] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [19] GenomicRanges_1.42.0 GenomeInfoDb_1.26.2
## [21] IRanges_2.24.1 S4Vectors_0.28.1
## [23] BiocGenerics_0.36.0 MatrixGenerics_1.2.1
## [25] matrixStats_0.58.0
##
## loaded via a namespace (and not attached):
## [1] circlize_0.4.12 plyr_1.8.6 lazyeval_0.2.2
## [4] shinydashboard_0.7.1 splines_4.0.3 crosstalk_1.1.1
## [7] digest_0.6.27 foreach_1.5.1 htmltools_0.5.1.1
## [10] fansi_0.4.2 magrittr_2.0.1 memoise_2.0.0
## [13] cluster_2.1.1 openxlsx_4.2.3 ComplexHeatmap_2.6.2
## [16] annotate_1.68.0 colorspace_2.0-0 blob_1.2.1
## [19] xfun_0.21 crayon_1.4.1 RCurl_1.98-1.2
## [22] jsonlite_1.7.2 survival_3.2-7 iterators_1.0.13
## [25] glue_1.4.2 registry_0.5-1 gtable_0.3.0
## [28] zlibbioc_1.36.0 XVector_0.30.0 UpSetR_1.4.0
## [31] GetoptLong_1.0.5 DelayedArray_0.16.1 V8_3.4.0
## [34] shape_1.4.5 scales_1.1.1 futile.options_1.0.1
## [37] DBI_1.1.1 randomcoloR_1.1.0.1 Rcpp_1.0.6
## [40] viridisLite_0.3.0 xtable_1.8-4 clue_0.3-58
## [43] bit_4.0.4 htmlwidgets_1.5.3 httr_1.4.2
## [46] RColorBrewer_1.1-2 ellipsis_0.3.1 farver_2.1.0
## [49] pkgconfig_2.0.3 XML_3.99-0.5 sass_0.3.1
## [52] locfit_1.5-9.4 utf8_1.2.1 labeling_0.4.2
## [55] tidyselect_1.1.0 rlang_0.4.10 later_1.1.0.1
## [58] AnnotationDbi_1.52.0 munsell_0.5.0 tools_4.0.3
## [61] cachem_1.0.4 generics_0.1.0 RSQLite_2.2.3
## [64] evaluate_0.14 stringr_1.4.0 fastmap_1.1.0
## [67] yaml_2.2.1 knitr_1.31 bit64_4.0.5
## [70] shinycssloaders_1.0.0 zip_2.1.1 purrr_0.3.4
## [73] mime_0.10 formatR_1.7 compiler_4.0.3
## [76] curl_4.3 png_0.1-7 tibble_3.1.0
## [79] bslib_0.2.4 stringi_1.5.3 highr_0.8
## [82] futile.logger_1.4.3 lattice_0.20-41 Matrix_1.3-2
## [85] ggsci_2.9 vctrs_0.3.7 pillar_1.5.1
## [88] lifecycle_1.0.0 jquerylib_0.1.3 GlobalOptions_0.1.2
## [91] cowplot_1.1.1 data.table_1.14.0 bitops_1.0-6
## [94] seriation_1.2-9 httpuv_1.5.5 R6_2.5.0
## [97] promises_1.2.0.1 TSP_1.1-10 gridExtra_2.3
## [100] codetools_0.2-18 lambda.r_1.2.4 assertthat_0.2.1
## [103] rjson_0.2.20 withr_2.4.1 GenomeInfoDbData_1.2.4
## [106] VennDiagram_1.6.20 grid_4.0.3 tidyr_1.1.3
## [109] rmarkdown_2.7 Cairo_1.5-12.2 Rtsne_0.15
## [112] shiny_1.6.0