SCLERORES - Identification of candidate genes
Thomas Bergmann
22 Februar, 2021
Abbreviations:
- B. villosa = B1896
- B. oleracea = B1909
1 Data import
Load libraries:
## Loading required package: S4Vectors## Loading required package: stats4## Loading required package: BiocGenerics## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which, which.max, which.min##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:base':
##
## expand.grid## Loading required package: IRanges##
## Attaching package: 'IRanges'## The following object is masked from 'package:grDevices':
##
## windows## Loading required package: GenomicRanges## Loading required package: GenomeInfoDb## Loading required package: SummarizedExperiment## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.## Loading required package: DelayedArray## Loading required package: matrixStats##
## Attaching package: 'matrixStats'## The following objects are masked from 'package:Biobase':
##
## anyMissing, rowMedians## Loading required package: BiocParallel##
## Attaching package: 'DelayedArray'## The following objects are masked from 'package:matrixStats':
##
## colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges## The following objects are masked from 'package:base':
##
## aperm, apply, rowsum## gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.## ## gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.##
## Attaching package: 'gdata'## The following object is masked from 'package:SummarizedExperiment':
##
## trim## The following object is masked from 'package:Biobase':
##
## combine## The following object is masked from 'package:GenomicRanges':
##
## trim## The following object is masked from 'package:IRanges':
##
## trim## The following objects are masked from 'package:S4Vectors':
##
## first, first<-## The following object is masked from 'package:BiocGenerics':
##
## combine## The following object is masked from 'package:stats4':
##
## nobs## The following object is masked from 'package:stats':
##
## nobs## The following object is masked from 'package:utils':
##
## object.size## The following object is masked from 'package:base':
##
## startsWith##
## Attaching package: 'gridExtra'## The following object is masked from 'package:gdata':
##
## combine## The following object is masked from 'package:Biobase':
##
## combine## The following object is masked from 'package:BiocGenerics':
##
## combine##
## Attaching package: 'ggpubr'## The following object is masked from 'package:cowplot':
##
## get_legend## Loading required package: grid## Loading required package: futile.logger##
## Attaching package: 'VennDiagram'## The following object is masked from 'package:ggpubr':
##
## rotate## Registered S3 method overwritten by 'gplots':
## method from
## reorder.factor gdata##
## Attaching package: 'gplots'## The following object is masked from 'package:gdata':
##
## reorder.factor## The following object is masked from 'package:IRanges':
##
## space## The following object is masked from 'package:S4Vectors':
##
## space## The following object is masked from 'package:stats':
##
## lowess##
## Attaching package: 'dplyr'## The following object is masked from 'package:gridExtra':
##
## combine## The following objects are masked from 'package:gdata':
##
## combine, first, last## The following object is masked from 'package:matrixStats':
##
## count## The following object is masked from 'package:Biobase':
##
## combine## The following objects are masked from 'package:GenomicRanges':
##
## intersect, setdiff, union## The following object is masked from 'package:GenomeInfoDb':
##
## intersect## The following objects are masked from 'package:IRanges':
##
## collapse, desc, intersect, setdiff, slice, union## The following objects are masked from 'package:S4Vectors':
##
## first, intersect, rename, setdiff, setequal, union## The following objects are masked from 'package:BiocGenerics':
##
## combine, intersect, setdiff, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## ========================================
## ComplexHeatmap version 2.2.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
## ========================================##
## Attaching package: 'officer'## The following object is masked from 'package:readxl':
##
## read_xlsx##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rows## Registering fonts with R## Scanning ttf files in C:\WINDOWS\Fonts ...## Extracting .afm files from .ttf files...## C:\Windows\Fonts\arial.ttf : ArialMT already registered in fonts database. Skipping.
## C:\Windows\Fonts\arialbd.ttf : Arial-BoldMT already registered in fonts database. Skipping.
## C:\Windows\Fonts\arialbi.ttf : Arial-BoldItalicMT already registered in fonts database. Skipping.
## C:\Windows\Fonts\ariali.ttf : Arial-ItalicMT already registered in fonts database. Skipping.
## Found FontName for 0 fonts.
## Scanning afm files in C:/Users/Thomas/Documents/R/win-library/3.6/extrafontdb/metricsImport and filter data:
## Data import
B1896_Genes <- read.table("./DESeq2/B1896_Genes.txt", sep="\t", header=TRUE, quote="")
B1909_Genes <- read.table("./DESeq2/B1909_Genes.txt", sep="\t", header=TRUE, quote="")
Interaction_Genes <- read.table("./DESeq2/Interaction_Genes.txt", sep="\t", header=TRUE, quote="")
GO_Annotations <- read.table("./Annotations/gene2GO.txt", sep="\t", quote="")
GO_Annotations <- GO_Annotations %>%
dplyr::select(V1, V5) %>%
dplyr::rename(Major_Iso = "V1", GO = "V5")
## Load geneID2GO-mappinf file
load(file="./Annotations/geneID2GO.RData")
## Filter for de novo transcripts
Trinity_B1896 <- B1896_Genes[grep("TRINITY", B1896_Genes$Gene_ID),]
##
Trinity_B1909 <- B1909_Genes[grep("TRINITY", B1909_Genes$Gene_ID),]
## Interaction
Trinity_Interaction <- Interaction_Genes[grep("TRINITY", Interaction_Genes$Gene_ID),]
Trinity_Interaction <- Trinity_Interaction %>% dplyr::filter(Significant == "Yes")Prepare venn diagram:
wf <- windowsFonts()
## Prepare the venn diagramm
B1896_DEGs_UP <- subset(Trinity_B1896, Regulation == "Up" & Significant == "Yes")
B1896_DEGs_DOWN <- subset(Trinity_B1896, Regulation == "Down" & Significant == "Yes")
B1909_DEGs_UP <- subset(Trinity_B1909, Regulation == "Up" & Significant == "Yes")
B1909_DEGs_DOWN <- subset(Trinity_B1909, Regulation == "Down" & Significant == "Yes")
## Venn diagram
UP <- VennDiagram::venn.diagram(list(B1896=B1896_DEGs_UP$Gene_ID, B1909=B1909_DEGs_UP$Gene_ID),
alpha = c(0.7,0.7), fill=c("darkslategray4", "deepskyblue2"), filename=NULL,
cex=0.8, cat.cex=0.8, cat.fontface="bold",
fontfamily=wf$Arial,
fontface="bold",
main.fontface="bold",
main.cex=1.1,
main.pos=c(0.5,1),
cat.fontfamily=wf$Arial,
main.fontfamily=wf$Arial,
cat.pos=c(-160,160),
cat.dist=c(0.04,0.04),
margin=0.05,
main="Upregulated")
DOWN <- VennDiagram::venn.diagram(list(B1896=B1896_DEGs_DOWN$Gene_ID, B1909=B1909_DEGs_DOWN$Gene_ID),
alpha = c(0.7,0.7), fill=c("darkslategray4", "deepskyblue2"), filename=NULL,
cex=0.8, cat.cex=0.8, cat.fontface="bold",
fontfamily=wf$Arial,
fontface="bold",
main.fontface="bold",
main.cex=1.1,
main.pos=c(0.5,1),
cat.fontfamily=wf$Arial,
main.fontfamily=wf$Arial,
cat.pos=c(-160,170),
cat.dist=c(0.04,0.04),
margin=0.03,
main="Downregulated")
Venn <- cowplot::plot_grid(grobTree(UP), grobTree(DOWN), ncol=2, nrow=1, align ="v", scale=0.9)Plot:
Preparation of tables:
## get IDs
venn.list <- list(B1896=B1896_DEGs_UP$Gene_ID, B1909=B1909_DEGs_UP$Gene_ID)
IDs_List <- gplots::venn(venn.list, show.plot=FALSE)
intersections_IDs <- attr(IDs_List, "intersections")
## Prepare background
All_IDs <- intersections_IDs[c(1)]
Resistant_IDs <- intersections_IDs[c(2)]
Susceptible_IDs <- intersections_IDs[c(3)]
## Convert to data frame
Resistant_B1896_IDs <- as.data.frame(Resistant_IDs) %>%
dplyr::rename(Gene_ID = "B1896")
Susceptible_B1909_IDs <- as.data.frame(Susceptible_IDs) %>%
dplyr::rename(Gene_ID = "B1909")
Common_IDs <- as.data.frame(All_IDs) %>%
dplyr::rename(Gene_ID = "B1896.B1909")
## Make common table
Background_IDs <- rbind(Resistant_B1896_IDs, Susceptible_B1909_IDs, Common_IDs)
## Add length info
Background_IDs <- dplyr::left_join(Background_IDs, Trinity_B1896, by = "Gene_ID") %>%
dplyr::select(Gene_ID, Transcript_Length)
Length <- as.numeric(Background_IDs$Transcript_Length)
## Convert
Background_IDs <- as.character(Background_IDs$Gene_ID)
Resistant_B1896_IDs <- as.character(Resistant_B1896_IDs$Gene_ID)
## Mark specific DEGs in background
B1896List <- as.integer(Background_IDs %in% Resistant_B1896_IDs)
names(B1896List) <- Background_IDs2 GO enrichment
Load GOseq-package:
## Loading required package: BiasedUrn## Loading required package: geneLenDataBase##
## Attaching package: 'geneLenDataBase'## The following object is masked from 'package:S4Vectors':
##
## unfactor## Calculate PWF:
## DEgenes bias.data pwf
## TRINITY_DN230_c0_g1_i15 1 1332 0.17791633
## TRINITY_DN12657_c0_g1_i3 1 828 0.26045442
## TRINITY_DN1285_c0_g2_i6 1 1623 0.12844429
## TRINITY_DN237_c0_g1_i19 1 1438 0.15924759
## TRINITY_DN6020_c1_g1_i1 1 1323 0.17950851
## TRINITY_DN2849_c0_g1_i8 1 1900 0.09284003Perform analysis and adjust P-values according to Benjamini & Hochberg:
## Perform analysis
B1896_Specific_GO=goseq::goseq(pwfB1896, gene2cat=geneID2GO, use_genes_without_cat=TRUE)## Using manually entered categories.## Calculating the p-values...## 'select()' returned 1:1 mapping between keys and columns## Function to adjust via FDR
FDRcalculater <- function(GOtable){
FDR_Over <- p.adjust(GOtable$over_represented_pvalue, method="BH")
as.data.frame(FDR_Over)
FDR_Under <- p.adjust(GOtable$under_represented_pvalue, method="BH")
as.data.frame(FDR_Under)
GOtableparsed <- cbind(GOtable, FDR_Over, FDR_Under)
return(GOtableparsed)
}
## Apply
B1896_Specific_GO <- FDRcalculater(B1896_Specific_GO)
## Prepare tables
B1896_Specific_GO <- B1896_Specific_GO %>%
dplyr::select(category, numDEInCat, numInCat, term, ontology, FDR_Over, FDR_Under) %>%
dplyr::rename(Category="category", DEGs = "numDEInCat", Genes = "numInCat", Term ="term", Ontology ="ontology")
## Remove 'NAs' row from the tables
B1896_Specific_GO <- B1896_Specific_GO %>% na.omit(cols="Category")
## Over stats
B1896_GOOverStats <- as.data.frame(table(B1896_Specific_GO$FDR_Over < 0.05))
colnames(B1896_GOOverStats) <- c("B1896", "Overrepresented")
## Under stats
B1896_GOUnderStats <- as.data.frame(table(B1896_Specific_GO$FDR_Under < 0.05))
colnames(B1896_GOUnderStats) <- c("B1896", "Underrepresented")
## Merge
B1896_GOStats <- cbind(B1896_GOOverStats[2,],B1896_GOUnderStats[2,])
rownames(B1896_GOStats) <- c("B1896")
B1896_GOStats <- subset(B1896_GOStats, select= -c(1,3))
## Stats
knitr::kable(B1896_GOStats) %>% kableExtra::kable_styling(bootstrap_options = c("striped"))| Overrepresented | Underrepresented | |
|---|---|---|
| B1896 | 170 | 1 |
2.1 Results
Subdivide the results by their GO category:
##
Specific_BP <- B1896_Specific_GO %>% filter(B1896_Specific_GO$Ontology == "BP")
Specific_CC <- B1896_Specific_GO %>% filter(B1896_Specific_GO$Ontology == "CC")
Specific_MF <- B1896_Specific_GO %>% filter(B1896_Specific_GO$Ontology == "MF")Export tables:
write.table(Specific_BP, "./GOseq/Specific_Upregulated_TrinityBP.txt", sep="\t", quote=FALSE, row.names=FALSE)
write.table(Specific_MF, "./GOseq/Specific_Upregulated_TrinityMF.txt", sep="\t", quote=FALSE, row.names=FALSE)
write.table(Specific_CC, "./GOseq/Specific_Upregulated_TrinityCC.txt", sep="\t", quote=FALSE, row.names=FALSE)Retrieve defense respone genes:
DEGs <- as.data.frame(Resistant_B1896_IDs)
colnames(DEGs) <- c("Major_Iso")
DEGs <- dplyr::left_join(DEGs, Trinity_B1896, by = "Major_Iso") %>%
dplyr::select(Major_Iso, Log2_FC) %>%
dplyr::rename(B1896_FC = "Log2_FC") %>%
dplyr::left_join(., Trinity_B1909, by = "Major_Iso") %>%
dplyr::rename(B1909_FC = "Log2_FC") %>%
dplyr::select(Major_Iso, B1896_FC, B1909_FC) %>%
dplyr::left_join(., GO_Annotations, by = "Major_Iso")
Defense_Response <- DEGs[grep("GO:0006952", DEGs$GO),]
## Prepare the data for the heatmap
rownames(Defense_Response) <- Defense_Response$Major_Iso
Defense_Response <- Defense_Response %>% dplyr::select(B1896_FC, B1909_FC)
Defense_Response <- as.matrix(Defense_Response, rownames.force=TRUE)
x11()Heatmap:
coul <- circlize::colorRamp2(c(-4,-2,0,2,4,6,8),c("darkslategray","deepskyblue2","white", "firebrick1","firebrick2","firebrick3","firebrick4"))
Heatmap_Defense_Response <- ComplexHeatmap::Heatmap(Defense_Response, col=coul,
show_row_names=TRUE, row_dend_width=unit(2,"cm"),
column_dend_height=unit(1,"cm"),
heatmap_legend_param=list(
title=expression(bold("LOG"[2]*"-Fold change")),
title_position="leftcenter-rot",
legend_height=unit(6,"cm"),
labels_gp=gpar(fontsize=10,fontface="bold"),
at=c(-4,-2,0,2,4,6,8),
border="black"),
border=TRUE,
column_title="",
row_title=NULL,
cluster_columns=FALSE,
column_km = 1,
row_gap = unit(0.2, "cm"),
row_dend_reorder=TRUE,
row_names_gp=gpar(fontsize=5.5, fontface="bold"),
row_names_rot=0,
column_names_rot=45,
column_names_side="top",
column_labels=c("BRA1896", "BRA1909"),
column_names_centered = TRUE,
column_names_gp=gpar(fontsize=8, fontface="bold"),
left_annotation = rowAnnotation(foo=anno_block(gp=gpar(fill="cyan2"),
labels=c("I","II","III","IV","V"),
labels_gp = gpar(col = "black", fontsize = 10, fontface="bold"))),
row_km=5,
row_km_repeats=100,
width=unit(1.75, "cm"),
height=unit(20,"cm"),
show_parent_dend_line=FALSE,
use_raster=FALSE,
raster_quality=20)
Candidates <- ComplexHeatmap::draw(Heatmap_Defense_Response, heatmap_legend_side="left", padding=unit(c(1,1,1,1), "cm"))
Prepare information:
## Check cluster
ListClusterI <- unlist(row_order(Heatmap_Defense_Response)[1])
names(ListClusterI) <- NULL
ListClusterII <- unlist(row_order(Heatmap_Defense_Response)[2])
names(ListClusterII) <- NULL
ListClusterIII <- unlist(row_order(Heatmap_Defense_Response)[3])
names(ListClusterIII) <- NULL
ListClusterIV <- unlist(row_order(Heatmap_Defense_Response)[4])
names(ListClusterIV) <- NULL
ListClusterV <- unlist(row_order(Heatmap_Defense_Response)[5])
names(ListClusterV) <- NULL
ClusterI <- Defense_Response[sort(ListClusterI),]
ClusterII <- Defense_Response[sort(ListClusterII),]
ClusterIII <- Defense_Response[sort(ListClusterIII),]
ClusterIV <- Defense_Response[sort(ListClusterIV),]
ClusterV <- Defense_Response[sort(ListClusterV),]Cluster I:
DT::datatable(as.data.frame(ClusterI, filter="top", options=list(scrollX=T))) %>%
DT::formatRound(c("B1896_FC","B1909_FC"), 3)Cluster II:
DT::datatable(as.data.frame(ClusterII, filter="top", options=list(scrollX=T))) %>%
DT::formatRound(c("B1896_FC","B1909_FC"), 3)Cluster III:
DT::datatable(as.data.frame(ClusterIII, filter="top", options=list(scrollX=T))) %>%
DT::formatRound(c("B1896_FC","B1909_FC"), 3)Cluster IV:
DT::datatable(as.data.frame(ClusterIV, filter="top", options=list(scrollX=T))) %>%
DT::formatRound(c("B1896_FC","B1909_FC"), 3)Cluster V:
DT::datatable(as.data.frame(ClusterV, filter="top", options=list(scrollX=T))) %>%
DT::formatRound(c("B1896_FC","B1909_FC"), 3)Finished!
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19041)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252
## [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
## [5] LC_TIME=German_Germany.1252
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] goseq_1.38.0 geneLenDataBase_1.22.0
## [3] BiasedUrn_1.07 extrafont_0.17
## [5] kableExtra_1.3.1 RColorBrewer_1.1-2
## [7] ggrepel_0.8.2 officer_0.3.15
## [9] ComplexHeatmap_2.2.0 dplyr_1.0.2
## [11] gplots_3.1.1 VennDiagram_1.6.20
## [13] futile.logger_1.4.3 ggpubr_0.4.0
## [15] gridExtra_2.3 readxl_1.3.1
## [17] gdata_2.18.0 cowplot_1.1.0
## [19] ggplot2_3.3.2 DESeq2_1.26.0
## [21] SummarizedExperiment_1.16.1 DelayedArray_0.12.3
## [23] BiocParallel_1.20.1 matrixStats_0.57.0
## [25] Biobase_2.46.0 GenomicRanges_1.38.0
## [27] GenomeInfoDb_1.22.1 IRanges_2.20.2
## [29] S4Vectors_0.24.4 BiocGenerics_0.32.0
##
## loaded via a namespace (and not attached):
## [1] uuid_0.1-4 backports_1.2.0 circlize_0.4.11
## [4] Hmisc_4.4-2 BiocFileCache_1.10.2 splines_3.6.3
## [7] crosstalk_1.1.0.1 digest_0.6.27 htmltools_0.5.0
## [10] GO.db_3.10.0 magrittr_2.0.1 checkmate_2.0.0
## [13] memoise_1.1.0 cluster_2.1.0 openxlsx_4.2.3
## [16] Biostrings_2.54.0 annotate_1.64.0 extrafontdb_1.0
## [19] askpass_1.1 prettyunits_1.1.1 jpeg_0.1-8.1
## [22] colorspace_2.0-0 rappdirs_0.3.1 blob_1.2.1
## [25] rvest_0.3.6 haven_2.3.1 xfun_0.19
## [28] jsonlite_1.7.2 crayon_1.3.4 RCurl_1.98-1.2
## [31] genefilter_1.68.0 survival_3.1-8 glue_1.4.2
## [34] gtable_0.3.0 zlibbioc_1.32.0 XVector_0.26.0
## [37] webshot_0.5.2 GetoptLong_1.0.4 car_3.0-10
## [40] Rttf2pt1_1.3.8 shape_1.4.5 abind_1.4-5
## [43] scales_1.1.1 futile.options_1.0.1 DBI_1.1.0
## [46] rstatix_0.6.0 Rcpp_1.0.5 progress_1.2.2
## [49] viridisLite_0.3.0 xtable_1.8-4 htmlTable_2.1.0
## [52] clue_0.3-58 foreign_0.8-75 bit_4.0.4
## [55] Formula_1.2-4 DT_0.16 htmlwidgets_1.5.3
## [58] httr_1.4.2 ellipsis_0.3.1 pkgconfig_2.0.3
## [61] XML_3.99-0.3 dbplyr_2.0.0 nnet_7.3-12
## [64] locfit_1.5-9.4 tidyselect_1.1.0 rlang_0.4.9
## [67] AnnotationDbi_1.48.0 munsell_0.5.0 cellranger_1.1.0
## [70] tools_3.6.3 generics_0.1.0 RSQLite_2.2.1
## [73] broom_0.7.2 evaluate_0.14 stringr_1.4.0
## [76] yaml_2.2.1 knitr_1.30 bit64_4.0.5
## [79] zip_2.1.1 caTools_1.18.0 purrr_0.3.4
## [82] nlme_3.1-144 formatR_1.7 xml2_1.3.2
## [85] biomaRt_2.42.1 compiler_3.6.3 rstudioapi_0.13
## [88] curl_4.3 png_0.1-7 ggsignif_0.6.0
## [91] tibble_3.0.4 geneplotter_1.64.0 stringi_1.5.3
## [94] highr_0.8 GenomicFeatures_1.38.2 forcats_0.5.0
## [97] lattice_0.20-38 Matrix_1.2-18 vctrs_0.3.5
## [100] pillar_1.4.7 lifecycle_0.2.0 GlobalOptions_0.1.2
## [103] data.table_1.13.4 bitops_1.0-6 rtracklayer_1.46.0
## [106] R6_2.5.0 latticeExtra_0.6-29 KernSmooth_2.23-16
## [109] rio_0.5.16 lambda.r_1.2.4 assertthat_0.2.1
## [112] gtools_3.8.2 openssl_1.4.3 rjson_0.2.20
## [115] withr_2.3.0 GenomicAlignments_1.22.1 Rsamtools_2.2.3
## [118] GenomeInfoDbData_1.2.2 mgcv_1.8-31 hms_0.5.3
## [121] rpart_4.1-15 tidyr_1.1.2 rmarkdown_2.5
## [124] carData_3.0-4 base64enc_0.1-3