Libraries and files loading
library(readxl)
library(metap)
library(writexl)
library(org.Hs.eg.db)
library(enrichplot)
library(clusterProfiler)
library(ggplot2)
hs=get("org.Hs.eg.db")
MEDIAsave="/home/mclaudia/MEDIA Project/Codici_finali/"
MEDIAgraph="/home/mclaudia/MEDIA Project/Codici_finali/Paper_figures/"
load(paste0(MEDIAsave,"DE_DESeq2_SCell.RData"))
singleCell <- df.res_sc
singleCell$gene_name <-rownames(singleCell)
load(paste0(MEDIAsave,"genide-bulk-unpaired_DE.RData"))
unpaired <-OAvsnonOA
load(paste0(MEDIAsave,"results_DE_EGApaired.RData"))
paired <-df.res
rm(df.res,OAvsnonOA,df.res_sc)
Retain up and down regulated genes for each signature
up_sc<-singleCell[singleCell$log2FoldChange >= 1.5 & singleCell$padj <=0.05,c(2,4,5,6)] #156
dw_sc <-singleCell[singleCell$log2FoldChange <= -1.5 & singleCell$padj <=0.05,c(2,4,5,6)] #79
up_paired <-paired[paired$log2FoldChange >= 0.8 & paired$padj <=0.05,c(3,6:8)] #190
dw_paired <-paired[paired$log2FoldChange <= -0.8 & paired$padj <=0.05,c(3,6:8)] #28
up_unpaired <-unpaired[unpaired$log2FoldChange >= 1.5 & unpaired$padj <=0.05,c(3,6:8)] #272
dw_unpaired <-unpaired[unpaired$log2FoldChange <= -1.5 & unpaired$padj <=0.05,c(3,6:8)] #170
DE_sc <-rbind(up_sc,dw_sc)
DE_paired <-rbind(up_paired,dw_paired)
DE_unpaired <-rbind(up_unpaired,dw_unpaired)
DE_sc <-DE_sc[order(DE_sc$log2FoldChange,decreasing = TRUE),]
DE_paired <-DE_paired[order(DE_paired$log2FoldChange,decreasing = TRUE),]
DE_unpaired <-DE_unpaired[order(DE_unpaired$log2FoldChange,decreasing = TRUE),]
Save files
write_xlsx(list(DE_sc,DE_paired,DE_unpaired),path=paste0(MEDIAsave,"up_dw_genes_eachSignature.xlsx"))
Pre-processing of the datasets: combine all the results, pValues and
log(FC)
paired <-paired[,c(3,6,8)]
unpaired <-unpaired[,c(3,6,8)]
singleCell <-singleCell[,c(2,4)]
singleCell <-cbind(gene_name=rownames(singleCell),singleCell)
ddsc <-singleCell$gene_name[grep("\\.",singleCell$gene_name)]
#check for coding duplicates not previously filtered
ss=ddsc[-grep("^AP0.+|^RP.+-|^CT.-.+|^AC.+|^AL.|^GS1-.+|^XX.+|^LL.+",ddsc)]
tmmp1 <-singleCell[singleCell$gene_name %in% ss,]
tm1 <- gsub("\\..*","",tmmp1$gene_name)
table(singleCell$gene_name[singleCell$gene_name %in% tm1])
##
## CRYBG3 SRSF10
## 1 1
#CRYBG3 1 SRSF10 1, remove two duplicates and correct the gene name for 16 coding genes
ss2=singleCell[singleCell$gene_name %in% ss,]
singleCell <-singleCell[!singleCell$gene_name %in% c("CRYBG3.1", "SRSF10.1"),]
ss <-ss[!ss%in%c("CRYBG3.1", "SRSF10.1")]
singleCell[singleCell$gene_name %in% ss,"gene_name"] <- gsub("\\..*","",ss)
rm(ss,ddsc,tm1)
colnames(unpaired) <- c("log2FC","p_val", "gene_name")
colnames(paired) <- c("log2FC","p_val", "gene_name")
colnames(singleCell)[c(2,3)] <- c("log2FC","p_val")
unn=unique(c(paired$gene_name, unpaired$gene_name, singleCell$gene_name)) #20878
comb <-matrix(NA,length(unn), 4)
colnames(comb) <-c("pv_sc","pv_unp", "pv_p","combined")
rownames(comb) <-unn
comb <-as.data.frame(comb)
comb[singleCell$gene_name,"pv_sc"] <- singleCell[,"p_val"]
comb[paired$gene_name,"pv_p"] <- paired[,"p_val"]
comb[unpaired$gene_name,"pv_unp"] <- unpaired[,"p_val"]
## remove genes with less than 2 p-values available
comb$num <-unlist(apply(comb,1,function(x) ifelse(sum(is.na(x)) > 2, x["num"] <-1, x["num"] <-0)))
comb <- comb[comb$num == 0,]
nrow(comb)
## [1] 14124
##prepare log2FC matrix, NA logs are set to 0
comb2 <-matrix(NA,length(unn), 3)
colnames(comb2) <-c("lg_sc","lg_unp", "lg_p")
rownames(comb2) <-unn
comb2 <-as.data.frame(comb2)
comb2[singleCell$gene_name,"lg_sc"] <- singleCell[,"log2FC"]
comb2[paired$gene_name,"lg_p"] <- paired[,"log2FC"]
comb2[unpaired$gene_name,"lg_unp"] <- unpaired[,"log2FC"]
comb2 <-comb2[rownames(comb),]
comb2[is.na(comb2)] <-0
comb2[1:10,]
## lg_sc lg_unp lg_p
## TMEM236 0.00000000 -0.06955428 -0.26787400
## CKS1B -0.43380952 -0.54173215 -0.41644124
## LSP1 0.02824830 0.90377762 0.47333724
## GOLGA7B 0.00000000 1.55892082 0.21265254
## ZNF33B -0.03770847 -1.35121872 -0.15997875
## POLR2J4 0.00000000 0.22663128 -0.06887927
## NPIPA7 0.00000000 0.33455492 -0.30219705
## CRYBG3 -0.05069620 0.40692647 0.11334711
## IDS -0.02367795 -0.24530223 0.09906970
## KIAA0391 0.00000000 0.01752132 -0.19913285
comb[1:10,]
## pv_sc pv_unp pv_p combined num
## TMEM236 NA 8.532968e-01 0.370640576 NA 0
## CKS1B 0.0008950385 6.264216e-07 0.005122658 NA 0
## LSP1 0.7988460535 2.025377e-07 0.024419136 NA 0
## GOLGA7B NA 4.204860e-04 0.251486402 NA 0
## ZNF33B 0.3086389326 9.853069e-11 0.235058259 NA 0
## POLR2J4 NA 3.500615e-01 0.592451405 NA 0
## NPIPA7 NA 3.583940e-01 0.002639714 NA 0
## CRYBG3 0.2460586015 6.063077e-06 0.307002996 NA 0
## IDS 0.8075623317 1.236367e-02 0.284627620 NA 0
## KIAA0391 NA 8.529354e-01 0.405906600 NA 0
Set genes with a -0.1 <= logFC <=0.1 to 0 and delete, where do
exist, 0 log genes in all experiments
for(i in 1:nrow(comb2)) for(k in 1:ncol(comb2[i,])) if(comb2[i,k] >= -0.1 & comb2[i,k] <= 0.1) comb2[i,k]<-0
comb2 <-comb2[!rowSums(comb2)==0,] #960 genes are removed
comb <-comb[rownames(comb2),]
Check for discordant genes
comb2$conc <-NA
for(i in 1:nrow(comb2)){
if((comb2[i,1] >= 0 & comb2[i,2] >=0 & comb2[i,3] >=0 )|(comb2[i,1] <=0 & comb2[i,2]<=0 & comb2[i,3] <=0)){
comb2[i,"conc"] <- 1
}else{
comb2[i,"conc"] <- 0
}
}
disc <-rownames(comb2[comb2$conc==0,])
length(disc)
## [1] 3879
Apply sum of logs (Fisher’s) method for pValue combination: for
discordant genes set combined pValue to 1
tmp <-apply(comb[,1:3], 1, function(x) sumlog(as.numeric(na.omit(x))))
comb$chisq <-unlist(tmp)[grep("\\.chisq",names(unlist(tmp)))]
comb$df<-unlist(tmp)[grep("\\.df",names(unlist(tmp)))]
comb$combined <-unlist(tmp)[grep("\\.p",names(unlist(tmp)))]
comb[rownames(comb) %in% disc, "combined"] <-1
comb2 <- comb2[rownames(comb),]
comb2$comb_pv <- comb$combined
comb2$md_log <-apply(comb2,1,function(x) median(x[1:3]))
comb2$name <-rownames(comb2)
combined <-comb2[,c(5:7)]
colnames(combined) <-c("pv_combined","median_logFC","gene")
univ <- rownames(combined) # universe for future enrichment
combined$pv_adj <-p.adjust(combined$pv_combined, method = "fdr")
combined[1:10,]
## pv_combined median_logFC gene pv_adj
## TMEM236 6.803429e-01 0.0000000 TMEM236 1.000000e+00
## CKS1B 1.093468e-09 -0.4338095 CKS1B 2.138843e-08
## LSP1 8.200047e-07 0.4733372 LSP1 9.040655e-06
## GOLGA7B 1.073799e-03 0.2126525 GOLGA7B 5.117847e-03
## ZNF33B 2.544684e-09 -0.1599787 ZNF33B 4.652531e-08
## POLR2J4 5.336534e-01 0.0000000 POLR2J4 8.361121e-01
## NPIPA7 1.000000e+00 0.0000000 NPIPA7 1.000000e+00
## CRYBG3 5.593369e-05 0.1133471 CRYBG3 3.875321e-04
## IDS 6.835336e-02 0.0000000 IDS 1.709353e-01
## KIAA0391 7.134405e-01 0.0000000 KIAA0391 1.000000e+00
Select consensus up-genes (39) and down-genes (5)
up <-combined[combined$median_logFC >= 1.5 & combined$pv_adj <= 0.01, ]
up <-up[order(up$pv_combined, decreasing = FALSE),]
dw <-combined[combined$median_logFC <= -1.5 & combined$pv_adj <= 0.01,]
dw <-dw[order(dw$pv_combined, decreasing = FALSE),]
consensus <-combined[c(rownames(up),rownames(dw)),]
print(consensus$gene)
## [1] "TSPAN2" "COL1A1" "TGFBI" "CRLF1" "ASPN" "TNFSF11" "GAS1"
## [8] "LOXL3" "CXCL14" "ADAMTS6" "SMOC1" "DYSF" "GFRA2" "ADAMTS2"
## [15] "THBS3" "MAMDC2" "COL15A1" "KCNN4" "OGN" "SLC2A12" "DNER"
## [22] "ACKR2" "TNC" "HTRA1" "TNNT3" "SEMA3D" "SYT8" "AQP1"
## [29] "TPPP3" "GALNT16" "TREM1" "THY1" "TMEM59L" "COL3A1" "SOX11"
## [36] "TNNI2" "S100A4" "CRABP2" "R3HDML" "STEAP4" "APOD" "MYOC"
## [43] "CHAC2" "BEX2"
up_genes <-up[,c(4,2,3,1)]
dw_genes <-dw[,c(4,2,3,1)]
up_genes[1:10,]
## pv_adj median_logFC gene pv_combined
## TSPAN2 7.989686e-78 1.685722 TSPAN2 6.069345e-82
## COL1A1 9.056705e-76 1.829887 COL1A1 1.375981e-79
## TGFBI 2.749927e-69 1.922828 TGFBI 6.266926e-73
## CRLF1 1.168456e-53 1.792391 CRLF1 4.438073e-57
## ASPN 1.005455e-52 2.096358 ASPN 4.582749e-56
## TNFSF11 8.214647e-51 3.040106 TNFSF11 4.368166e-54
## GAS1 1.928506e-46 1.521092 GAS1 1.464985e-49
## LOXL3 1.109822e-42 1.584828 LOXL3 9.273813e-46
## CXCL14 1.599133e-41 1.634375 CXCL14 1.579211e-44
## ADAMTS6 2.095952e-37 1.846219 ADAMTS6 2.388277e-40
dw_genes
## pv_adj median_logFC gene pv_combined
## STEAP4 1.830962e-35 -1.602582 STEAP4 2.642683e-38
## APOD 1.337696e-31 -2.194969 APOD 2.946914e-34
## MYOC 3.542249e-26 -1.631551 MYOC 1.426156e-28
## CHAC2 9.295129e-23 -2.192607 CHAC2 5.366377e-25
## BEX2 1.213766e-18 -1.591772 BEX2 1.124882e-20
rm(up,dw,comb2)
Save list
write_xlsx(list(up_genes,dw_genes),path=paste0(MEDIAsave,"up_dw_genes_consensus.xlsx"))
save(combined, file=paste0(MEDIAsave,"combined_full_list.RData"))
ORA up-genes consensus
oraGO <-enrichGO(up_genes$gene,hs,
universe=univ, ont="BP", keyType = "SYMBOL",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
pAdjustMethod = "fdr")
dotplot(oraGO, title="GO enhanced terms for Combined signature", showCategory = 30)
oraGO@result[1,2] <-"ECM org." #"extracellular matrix organization"
oraGO@result[2,2] <-"ECM str.org." #"extracellular structure organization"
oraGO@result[3,2] <-"Ext.enc.str.org." #"external encapsulating structure organization"
oraGO@result[5,2] <-"collagen fibril org." #"collagen fibril organization"
cnetplot(oraGO, cex.params = list(gene_label = 1.2, category_label =1),
color.params = list(edge=TRUE),
size="qvalue",showCategory = 6, max.overlaps=30)+
theme(plot.margin = unit(c(0.5,0,0,0.5), "cm"))
Session Info
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 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/liblapack.so.3
##
## 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
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ggplot2_3.4.3 clusterProfiler_4.4.4 enrichplot_1.21.0
## [4] org.Hs.eg.db_3.15.0 AnnotationDbi_1.58.0 IRanges_2.30.1
## [7] S4Vectors_0.34.0 Biobase_2.56.0 BiocGenerics_0.42.0
## [10] writexl_1.4.2 metap_1.8 readxl_1.4.3
##
## loaded via a namespace (and not attached):
## [1] shadowtext_0.1.2 fastmatch_1.1-4 sn_2.1.1
## [4] plyr_1.8.8 igraph_1.5.0 lazyeval_0.2.2
## [7] splines_4.2.1 BiocParallel_1.30.4 GenomeInfoDb_1.32.4
## [10] qqconf_1.3.2 TH.data_1.1-2 digest_0.6.35
## [13] yulab.utils_0.0.6 htmltools_0.5.6 GOSemSim_2.22.0
## [16] viridis_0.6.3 GO.db_3.15.0 fansi_1.0.4
## [19] magrittr_2.0.3 memoise_2.0.1 Biostrings_2.64.1
## [22] graphlayouts_1.0.0 sandwich_3.0-2 colorspace_2.1-0
## [25] blob_1.2.4 ggrepel_0.9.3 rbibutils_2.2.13
## [28] xfun_0.39 dplyr_1.1.3 crayon_1.5.2
## [31] RCurl_1.98-1.12 jsonlite_1.8.7 scatterpie_0.2.1
## [34] survival_3.5-7 zoo_1.8-12 ape_5.7-1
## [37] glue_1.6.2 polyclip_1.10-4 gtable_0.3.4
## [40] zlibbioc_1.42.0 XVector_0.36.0 scales_1.2.1
## [43] DOSE_3.22.1 mvtnorm_1.2-3 DBI_1.1.3
## [46] Rcpp_1.0.11 plotrix_3.8-2 viridisLite_0.4.2
## [49] gridGraphics_0.5-1 tidytree_0.4.2 bit_4.0.5
## [52] httr_1.4.7 fgsea_1.22.0 RColorBrewer_1.1-3
## [55] TFisher_0.2.0 pkgconfig_2.0.3 farver_2.1.1
## [58] sass_0.4.7 utf8_1.2.3 labeling_0.4.3
## [61] ggplotify_0.1.0 tidyselect_1.2.0 rlang_1.1.1
## [64] reshape2_1.4.4 munsell_0.5.0 cellranger_1.1.0
## [67] tools_4.2.1 cachem_1.0.8 downloader_0.4
## [70] cli_3.6.1 generics_0.1.3 RSQLite_2.3.1
## [73] mathjaxr_1.6-0 evaluate_0.21 stringr_1.5.0
## [76] fastmap_1.1.1 yaml_2.3.7 ggtree_3.6.2
## [79] knitr_1.43 bit64_4.0.5 tidygraph_1.2.3
## [82] purrr_1.0.1 KEGGREST_1.36.3 ggraph_2.1.0
## [85] nlme_3.1-163 aplot_0.1.10 DO.db_2.9
## [88] compiler_4.2.1 rstudioapi_0.15.0 png_0.1-8
## [91] treeio_1.20.2 tibble_3.2.1 tweenr_2.0.2
## [94] bslib_0.5.1 stringi_1.7.12 highr_0.10
## [97] lattice_0.21-8 Matrix_1.6-1 multtest_2.52.0
## [100] vctrs_0.6.3 mutoss_0.1-13 pillar_1.9.0
## [103] lifecycle_1.0.3 Rdpack_2.4 jquerylib_0.1.4
## [106] data.table_1.14.8 bitops_1.0-7 patchwork_1.1.3
## [109] qvalue_2.28.0 R6_2.5.1 gridExtra_2.3
## [112] codetools_0.2-19 MASS_7.3-60 withr_2.5.0
## [115] mnormt_2.1.1 multcomp_1.4-25 GenomeInfoDbData_1.2.8
## [118] parallel_4.2.1 grid_4.2.1 ggfun_0.1.0
## [121] tidyr_1.3.0 rmarkdown_2.24 ggforce_0.4.1
## [124] numDeriv_2016.8-1.1