Split-Splines with contrasts in Moanin (leaf, postflowering, Year3)
This document has last been compiled on 2021-12-10 13:31:58.
Data Entry
## Reading in metadata from file: results/Year3/data/leaf_meta.txt
## Reading in data from file: results/Year3/data/leaf_normCounts.txt## Total number of leaf samples: 299
## Total number of genes in leaf samples: 22318
## No additional bad samples to remove (probably removed during normalization)
## Removing the 46 samples identified as not part of the main experiment, from the leaf samples
## After filtering leaf samples: 253 samples, 22318 genesIn this document, we find genes that are globally DE, meaning that we search for any difference between the expression of the gene. We do this via splines, where for each gene we fit a spline function to the expression levels within each condition and compare whether those splines are statistically different between the conditions.
Here we use “split-splines” meaning we fit a different spline function to before and after the watering change (first time point after Day 56 for Preflowering and first time point after Day 63 for post-flowering).
(Note that currently we are not expanding the postflowering condition points into control like Nelle did…)
Create Moanin Object
## Dimensions of data matrix restricting to designated timepoints: 22318 x 82
## The unique time points are:
## 65, 66, 70, 77, 84, 91, 105
## The conditions and number of samples are:
##
## Control.BT642 Control.RT430 Postflowering.BT642 Postflowering.RT430
## 21 21 20 20## Moanin object on 82 samples containing the following information:
## Group variable given by 'Group' with the following levels:
## Control.BT642 Control.RT430 Postflowering.BT642 Postflowering.RT430
## 21 21 20 20
## Time variable given by 'Time.Point'
## Basis matrix with 16 basis_matrix functions
## Basis matrix was constructed with the following spline_formula
## ~Group:splines::bs(Time.Point, df = 3) + Group + 0
##
## Information about the data (a SummarizedExperiment object):
## class: SummarizedExperiment
## dim: 22318 82
## metadata(0):
## assays(1): ''
## rownames(22318): Sobic.001G000200.v3.1 Sobic.001G000400.v3.1 ...
## Sobic.K044413.v3.1 Sobic.K044418.v3.1
## rowData names(0):
## colnames(82): 0823187L02 0823187L19 ... 1002185L04 1002185L21
## colData names(66): Barcode Genotype ... Group WeeklyGroupWe separate out the results into genotype-differences, RT430, and BT642 results.
We save results into comma-deliminated file in results/Year3/DE_splitsplines in the following files:
## Saving matrix of results on 22318 genes to comma-deliminated file in results/Year3/DE_splitsplines.## [1] "leaf_GenoDiff_Postflowering_splitsplinesDE_all.csv"## [1] "leaf_RT430_Postflowering_splitsplinesDE_all.csv"## [1] "leaf_BT642_Postflowering_splitsplinesDE_all.csv"Summary of significant genes
## Number of significant genes:## BT642: 4508## RT430: 5227## GenoDiff: 331Visualization of BT642 Results
## Contrast matrix:
## Contrasts
## Levels Postflowering.BT642-Control.BT642
## Control.BT642 -1
## Control.RT430 0
## Postflowering.BT642 1
## Postflowering.RT430 0Volcano plots
We will plot volcano plots of the genes based on the calculated log-fold-change (lfc) and the adjusted p-value (qval).
volcano_plot_moanin(moanin_results_all, target, lfc_thres=1, pval_thres=0.05)
Splines for top 20 genes
Here we plot the spline fits for the top 20 genes.
top_de_genes <- plot_top_de_genes(moaninObj = moaninObj, moanin_results = moanin_results_all, target = target, n_top = 20,lfc_thres=1, pval_thres=0.05,drought=condition)## 189 selected after filtering based on p-value and log fold change
## Splines for the top 20 genes in Year3, BT642 Postflowering:top_de_genes## pval qval lfc
## Sobic.002G141800.v3.1 1.930941e-16 2.394045e-13 3.015771
## Sobic.009G116700.v3.1 5.307049e-13 1.076704e-10 2.905757
## Sobic.001G226600.v3.1 1.594625e-10 1.227146e-08 2.483257
## Sobic.004G115900.v3.1 7.157529e-14 2.249783e-11 2.270461
## Sobic.003G292400.v3.1 2.293563e-16 2.693971e-13 2.192735
## Sobic.010G133700.v3.1 1.593727e-13 4.391012e-11 2.163747
## Sobic.006G157700.v3.1 2.134697e-13 5.539538e-11 2.133975
## Sobic.001G089000.v3.1 7.465506e-11 6.239989e-09 2.131739
## Sobic.007G109800.v3.1 2.086134e-10 1.536510e-08 2.079154
## Sobic.004G295000.v3.1 7.397015e-12 9.487309e-10 2.063248
## Sobic.009G208500.v3.1 7.147762e-15 3.709688e-12 2.055623
## Sobic.001G483500.v3.1 3.059149e-12 4.491515e-10 2.051980
## Sobic.006G051100.v3.1 7.382290e-07 1.730573e-05 2.049536
## Sobic.004G089600.v3.1 9.268043e-14 2.795067e-11 2.031285
## Sobic.008G087500.v3.1 1.174970e-10 9.364927e-09 2.030903
## Sobic.007G003500.v3.1 1.799016e-11 1.930224e-09 2.024842
## Sobic.005G215200.v3.1 7.939196e-13 1.554202e-10 1.998953
## Sobic.002G316900.v3.1 4.496247e-15 2.711966e-12 1.956946
## Sobic.001G155900.v3.1 6.400338e-15 3.400865e-12 1.865645
## Sobic.001G418600.v3.1 5.474748e-06 9.575230e-05 1.849901
Visualization of RT430 Results
## Contrast matrix:
## Contrasts
## Levels Postflowering.RT430-Control.RT430
## Control.BT642 0
## Control.RT430 -1
## Postflowering.BT642 0
## Postflowering.RT430 1Volcano plots
We will plot volcano plots of the genes based on the calculated log-fold-change (lfc) and the adjusted p-value (qval).
volcano_plot_moanin(moanin_results_all, target, lfc_thres=1, pval_thres=0.05)
Splines for top 20 genes
Here we plot the spline fits for the top 20 genes.
top_de_genes <- plot_top_de_genes(moaninObj = moaninObj, moanin_results = moanin_results_all, target = target, n_top = 20,lfc_thres=1, pval_thres=0.05,drought=condition)## 250 selected after filtering based on p-value and log fold change
## Splines for the top 20 genes in Year3, RT430 Postflowering:top_de_genes## pval qval lfc
## Sobic.002G141800.v3.1 6.456995e-22 1.601119e-18 3.975225
## Sobic.009G116700.v3.1 2.660292e-16 1.060174e-13 3.425735
## Sobic.001G226600.v3.1 1.265442e-16 5.537424e-14 3.386727
## Sobic.007G109800.v3.1 1.709135e-17 1.155841e-14 3.162411
## Sobic.003G082000.v3.1 9.534062e-10 3.947526e-08 2.975549
## Sobic.003G149100.v3.1 2.468306e-15 6.800641e-13 2.711758
## Sobic.003G292400.v3.1 6.430588e-19 7.553233e-16 2.663890
## Sobic.009G215700.v3.1 3.501029e-14 6.404301e-12 2.656029
## Sobic.004G115900.v3.1 1.838256e-17 1.206599e-14 2.631958
## Sobic.004G295000.v3.1 2.998769e-15 8.063075e-13 2.631094
## Sobic.007G143200.v3.1 1.978313e-13 2.830129e-11 2.608977
## Sobic.010G133000.v3.1 1.813705e-24 8.095290e-21 2.572492
## Sobic.001G418600.v3.1 2.816673e-11 1.843393e-09 2.559749
## Sobic.003G397400.v3.1 1.422241e-11 1.044084e-09 2.531968
## Sobic.003G081900.v3.1 5.485579e-09 1.771659e-07 2.416898
## Sobic.010G092900.v3.1 2.154101e-16 8.902420e-14 2.388749
## Sobic.002G316900.v3.1 3.562135e-18 3.179847e-15 2.320404
## Sobic.001G438000.v3.1 4.087401e-06 5.947235e-05 2.293372
## Sobic.003G082100.v3.1 4.768830e-09 1.565088e-07 2.282603
## Sobic.003G173700.v3.1 1.056901e-12 1.116216e-10 2.178001
Visualization of Differences between genotypes
## Contrast matrix:
## Contrasts
## Levels Postflowering.BT642-Control.BT642-Postflowering.RT430+Control.RT430
## Control.BT642 -1
## Control.RT430 1
## Postflowering.BT642 1
## Postflowering.RT430 -1Volcano plots
We will plot volcano plots of the genes based on the calculated log-fold-change (lfc) and the adjusted p-value (qval).
volcano_plot_moanin(moanin_results_all, target, lfc_thres=1, pval_thres=0.05)
Splines for top 20 genes
Here we plot the spline fits for the top 20 genes.
top_de_genes <- plot_top_de_genes(moaninObj = moaninObj, moanin_results = moanin_results_all, target = target, n_top = 20,lfc_thres=1, pval_thres=0.05,drought=condition)## 18 selected after filtering based on p-value and log fold change
## No genes found matching criteria
## Splines for the top 18 genes in Year3, BT642_vs_RT430 Postflowering:top_de_genes## pval qval lfc
## Sobic.010G159500.v3.1 1.497769e-17 1.671285e-13 2.593039
## Sobic.003G097800.v3.1 4.389070e-16 1.399298e-12 2.457680
## Sobic.006G025800.v3.1 1.531894e-10 1.424470e-07 1.573401
## Sobic.002G109100.v3.1 1.718179e-06 3.912715e-04 1.508651
## Sobic.010G057800.v3.1 3.314273e-12 4.622790e-09 1.445160
## Sobic.005G082600.v3.1 5.192356e-05 6.737082e-03 1.346243
## Sobic.001G280700.v3.1 1.268997e-06 2.981075e-04 1.234420
## Sobic.006G243200.v3.1 7.611279e-08 2.878999e-05 1.188933
## Sobic.010G253500.v3.1 4.773851e-04 3.712126e-02 1.183594
## Sobic.003G097900.v3.1 4.189031e-08 1.763898e-05 1.173254
## Sobic.001G377500.v3.1 8.866132e-09 5.041153e-06 1.170536
## Sobic.003G096500.v3.1 1.563580e-16 5.815735e-13 1.164332
## Sobic.004G214800.v3.1 2.313077e-04 2.178099e-02 1.158371
## Sobic.001G299300.v3.1 3.682959e-09 2.417429e-06 1.148860
## Sobic.005G065000.v3.1 1.111151e-09 8.856267e-07 1.089148
## Sobic.003G246900.v3.1 2.623360e-04 2.429275e-02 1.011633
## Sobic.009G021300.v3.1 4.658749e-04 3.648046e-02 1.007578
## Sobic.003G288900.v3.1 2.776638e-04 2.518953e-02 1.006811
# Session Info
## [1] "2021-12-10 13:32:58 PST"## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04 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 splines stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] moanin_1.1.2 topGO_2.44.0
## [3] SparseM_1.81 GO.db_3.13.0
## [5] AnnotationDbi_1.56.1 graph_1.72.0
## [7] SummarizedExperiment_1.24.0 Biobase_2.54.0
## [9] GenomicRanges_1.46.0 GenomeInfoDb_1.30.0
## [11] IRanges_2.28.0 S4Vectors_0.32.2
## [13] BiocGenerics_0.40.0 MatrixGenerics_1.6.0
## [15] matrixStats_0.61.0 MASS_7.3-54
## [17] reshape2_1.4.4 forcats_0.5.1
## [19] stringr_1.4.0 purrr_0.3.4
## [21] readr_2.0.2 tidyr_1.1.4
## [23] tibble_3.1.6 tidyverse_1.3.1
## [25] dplyr_1.0.7 ggplot2_3.3.5
## [27] knitr_1.36 rmarkdown_2.11
## [29] SCF_4.1.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7 fs_1.5.0 bit64_4.0.5
## [4] lubridate_1.8.0 httr_1.4.2 tools_4.1.2
## [7] backports_1.3.0 bslib_0.3.1 utf8_1.2.2
## [10] R6_2.5.1 DBI_1.1.1 colorspace_2.0-2
## [13] withr_2.4.2 gridExtra_2.3 tidyselect_1.1.1
## [16] bit_4.0.4 compiler_4.1.2 cli_3.1.0
## [19] rvest_1.0.2 xml2_1.3.2 DelayedArray_0.20.0
## [22] sass_0.4.0 scales_1.1.1 digest_0.6.28
## [25] XVector_0.34.0 pkgconfig_2.0.3 htmltools_0.5.2
## [28] highr_0.9 limma_3.50.0 dbplyr_2.1.1
## [31] fastmap_1.1.0 rlang_0.4.12 readxl_1.3.1
## [34] rstudioapi_0.13 RSQLite_2.2.8 jquerylib_0.1.4
## [37] generics_0.1.1 jsonlite_1.7.2 gtools_3.9.2
## [40] RCurl_1.98-1.5 magrittr_2.0.1 GenomeInfoDbData_1.2.7
## [43] Matrix_1.3-4 Rcpp_1.0.7 munsell_0.5.0
## [46] fansi_0.5.0 viridis_0.6.2 lifecycle_1.0.1
## [49] edgeR_3.36.0 stringi_1.7.5 yaml_2.2.1
## [52] ClusterR_1.2.5 zlibbioc_1.40.0 plyr_1.8.6
## [55] blob_1.2.2 grid_4.1.2 crayon_1.4.2
## [58] lattice_0.20-45 Biostrings_2.60.2 haven_2.4.3
## [61] KEGGREST_1.34.0 hms_1.1.1 locfit_1.5-9.4
## [64] pillar_1.6.4 reprex_2.0.1 glue_1.5.0
## [67] evaluate_0.14 modelr_0.1.8 png_0.1-7
## [70] vctrs_0.3.8 tzdb_0.2.0 cellranger_1.1.0
## [73] gtable_0.3.0 assertthat_0.2.1 cachem_1.0.6
## [76] xfun_0.28 broom_0.7.10 viridisLite_0.4.0
## [79] memoise_2.0.0 gmp_0.6-2.1 ellipsis_0.3.2