NOTE: all plots are saved as objects to facilitate placement of individual plots into combined figures. The aesthetics of the plots have been modified to create these figures. If plotted directly in the Markdown document, the figures will likely look poor due to these changed aesthetics.
NOTE: commented out code is code that was not used in the analysis, but was kept for historical purposes. It does not need to be uncommented in order to replicate the analysis.
This pipeline was adapted in part from the training resources at https://github.com/hbctraining/DGE_workshop/tree/master/lessons.
Before starting, ensure that your working directory contains the following files and architecture:
## gencode_featurecounts.Rmatrix.txt
## metadata.csv
## gencode_Output_Files
Load the required libraries.
# BiocManager::install("DESeq2")
# BiocManager::install("GeneStructureTools")
# BiocManager::install("limma")
# BiocManager::install("apeglm")
# BiocManager::install("DEGreport")
# BiocManager::install("AnnotationHub")
# BiocManager::install("ensembldb")
# devtools::install_github("stephenturner/annotables")
library(DESeq2)
Loading required package: S4Vectors
Loading required package: stats4
Loading required package: BiocGenerics
Attaching package: ‘BiocGenerics’
The following objects are masked from ‘package:stats’:
IQR, mad, sd, var, xtabs
The following objects are masked from ‘package:base’:
anyDuplicated, aperm, 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.max, which.min
Attaching package: ‘S4Vectors’
The following object is masked from ‘package:utils’:
findMatches
The following objects are masked from ‘package:base’:
expand.grid, I, unname
Loading required package: IRanges
Loading required package: GenomicRanges
Loading required package: GenomeInfoDb
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats
Attaching package: ‘MatrixGenerics’
The following objects are masked from ‘package:matrixStats’:
colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins,
colCumprods, colCumsums, colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs, colMads,
colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds, colQuantiles, colRanges, colRanks,
colSdDiffs, colSds, colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads, colWeightedMeans,
colWeightedMedians, colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs,
rowIQRs, rowLogSumExps, rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2,
rowTabulates, rowVarDiffs, rowVars, rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
rowWeightedSds, rowWeightedVars
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")'.
Attaching package: ‘Biobase’
The following object is masked from ‘package:MatrixGenerics’:
rowMedians
The following objects are masked from ‘package:matrixStats’:
anyMissing, rowMedians
library(tidyverse)
── Attaching core tidyverse packages ───────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.2 ✔ readr 2.1.4
✔ forcats 1.0.0 ✔ stringr 1.5.0
✔ ggplot2 3.4.2 ✔ tibble 3.2.1
✔ lubridate 1.9.2 ✔ tidyr 1.3.0
✔ purrr 1.0.1 ── Conflicts ─────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ lubridate::%within%() masks IRanges::%within%()
✖ dplyr::collapse() masks IRanges::collapse()
✖ dplyr::combine() masks Biobase::combine(), BiocGenerics::combine()
✖ dplyr::count() masks matrixStats::count()
✖ dplyr::desc() masks IRanges::desc()
✖ tidyr::expand() masks S4Vectors::expand()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::first() masks S4Vectors::first()
✖ dplyr::lag() masks stats::lag()
✖ ggplot2::Position() masks BiocGenerics::Position(), base::Position()
✖ purrr::reduce() masks GenomicRanges::reduce(), IRanges::reduce()
✖ dplyr::rename() masks S4Vectors::rename()
✖ lubridate::second() masks S4Vectors::second()
✖ lubridate::second<-() masks S4Vectors::second<-()
✖ dplyr::slice() masks IRanges::slice()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(GeneStructureTools)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'data.table':
method from
print.data.table
library(limma)
Attaching package: 'limma'
The following object is masked from 'package:DESeq2':
plotMA
The following object is masked from 'package:BiocGenerics':
plotMA
library(ggplot2)
library(pheatmap)
library(apeglm)
library(ggrepel)
library(DEGreport)
library(AnnotationHub)
Loading required package: BiocFileCache
Loading required package: dbplyr
Attaching package: 'dbplyr'
The following objects are masked from 'package:dplyr':
ident, sql
Attaching package: 'AnnotationHub'
The following object is masked from 'package:Biobase':
cache
library(ensembldb)
Loading required package: GenomicFeatures
Loading required package: AnnotationDbi
Attaching package: 'AnnotationDbi'
The following object is masked from 'package:dplyr':
select
Loading required package: AnnotationFilter
Attaching package: 'ensembldb'
The following object is masked from 'package:dplyr':
filter
The following object is masked from 'package:stats':
filter
library(annotables)
library(pals)
library(patchwork)
library(ggpubr)
library(scales)
Attaching package: 'scales'
The following object is masked from 'package:purrr':
discard
The following object is masked from 'package:readr':
col_factor
library(ggiraphExtra)
library(flextable)
Attaching package: 'flextable'
The following objects are masked from 'package:ggpubr':
border, font, rotate
The following object is masked from 'package:purrr':
compose
The following object is masked from 'package:SummarizedExperiment':
width
The following object is masked from 'package:GenomicRanges':
width
The following object is masked from 'package:IRanges':
width
The following object is masked from 'package:S4Vectors':
width
The following object is masked from 'package:BiocGenerics':
width
library(pals)
library(patchwork)
library(ggiraph)
library(officer)
library(magrittr)
Attaching package: 'magrittr'
The following object is masked from 'package:AnnotationFilter':
not
The following object is masked from 'package:purrr':
set_names
The following object is masked from 'package:tidyr':
extract
The following object is masked from 'package:GenomicRanges':
subtract
citation("DESeq2")
To cite package 'DESeq2' in publications use:
Love, M.I., Huber, W., Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data
with DESeq2 Genome Biology 15(12):550 (2014)
A BibTeX entry for LaTeX users is
@Article{,
title = {Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2},
author = {Michael I. Love and Wolfgang Huber and Simon Anders},
year = {2014},
journal = {Genome Biology},
doi = {10.1186/s13059-014-0550-8},
volume = {15},
issue = {12},
pages = {550},
}
citation("tidyverse")
To cite package 'tidyverse' in publications use:
Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester
J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V,
Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the tidyverse." _Journal of Open
Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
A BibTeX entry for LaTeX users is
@Article{,
title = {Welcome to the {tidyverse}},
author = {Hadley Wickham and Mara Averick and Jennifer Bryan and Winston Chang and Lucy D'Agostino McGowan and Romain François and Garrett Grolemund and Alex Hayes and Lionel Henry and Jim Hester and Max Kuhn and Thomas Lin Pedersen and Evan Miller and Stephan Milton Bache and Kirill Müller and Jeroen Ooms and David Robinson and Dana Paige Seidel and Vitalie Spinu and Kohske Takahashi and Davis Vaughan and Claus Wilke and Kara Woo and Hiroaki Yutani},
year = {2019},
journal = {Journal of Open Source Software},
volume = {4},
number = {43},
pages = {1686},
doi = {10.21105/joss.01686},
}
citation("GeneStructureTools")
To cite package 'GeneStructureTools' in publications use:
Signal B (2023). _GeneStructureTools: Tools for spliced gene structure manipulation and analysis_.
doi:10.18129/B9.bioc.GeneStructureTools <https://doi.org/10.18129/B9.bioc.GeneStructureTools>, R
package version 1.20.0, <https://bioconductor.org/packages/GeneStructureTools>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {GeneStructureTools: Tools for spliced gene structure manipulation and analysis},
author = {Beth Signal},
year = {2023},
note = {R package version 1.20.0},
url = {https://bioconductor.org/packages/GeneStructureTools},
doi = {10.18129/B9.bioc.GeneStructureTools},
}
ATTENTION: This citation information has been auto-generated from the package DESCRIPTION file and may
need manual editing, see 'help("citation")'.
citation("limma")
To cite package 'limma' in publications use:
Ritchie, M.E., Phipson, B., Wu, D., Hu, Y., Law, C.W., Shi, W., and Smyth, G.K. (2015). limma powers
differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research
43(7), e47.
A BibTeX entry for LaTeX users is
@Article{,
author = {Matthew E Ritchie and Belinda Phipson and Di Wu and Yifang Hu and Charity W Law and Wei Shi and Gordon K Smyth},
title = {{limma} powers differential expression analyses for {RNA}-sequencing and microarray studies},
journal = {Nucleic Acids Research},
year = {2015},
volume = {43},
number = {7},
pages = {e47},
doi = {10.1093/nar/gkv007},
}
citation("ggplot2")
To cite ggplot2 in publications, please use
H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.
A BibTeX entry for LaTeX users is
@Book{,
author = {Hadley Wickham},
title = {ggplot2: Elegant Graphics for Data Analysis},
publisher = {Springer-Verlag New York},
year = {2016},
isbn = {978-3-319-24277-4},
url = {https://ggplot2.tidyverse.org},
}
citation("pheatmap")
To cite package 'pheatmap' in publications use:
Kolde R (2019). _pheatmap: Pretty Heatmaps_. R package version 1.0.12,
<https://CRAN.R-project.org/package=pheatmap>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {pheatmap: Pretty Heatmaps},
author = {Raivo Kolde},
year = {2019},
note = {R package version 1.0.12},
url = {https://CRAN.R-project.org/package=pheatmap},
}
ATTENTION: This citation information has been auto-generated from the package DESCRIPTION file and may
need manual editing, see 'help("citation")'.
citation("apeglm")
To cite package 'apeglm' in publications use:
Zhu, A., Ibrahim, J.G., Love, M.I. Heavy-tailed prior distributions for sequence count data: removing
the noise and preserving large differences Bioinformatics (2018)
A BibTeX entry for LaTeX users is
@Article{,
title = {Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences},
author = {Anqi Zhu and Joseph G. Ibrahim and Michael I. Love},
year = {2018},
journal = {Bioinformatics},
doi = {10.1093/bioinformatics/bty895},
}
citation("ggrepel")
To cite package 'ggrepel' in publications use:
Slowikowski K (2023). _ggrepel: Automatically Position Non-Overlapping Text Labels with 'ggplot2'_. R
package version 0.9.3, <https://CRAN.R-project.org/package=ggrepel>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {ggrepel: Automatically Position Non-Overlapping Text Labels with
'ggplot2'},
author = {Kamil Slowikowski},
year = {2023},
note = {R package version 0.9.3},
url = {https://CRAN.R-project.org/package=ggrepel},
}
citation("DEGreport")
To cite package 'DEGreport' in publications use:
Pantano L (2023). _DEGreport: Report of DEG analysis_. doi:10.18129/B9.bioc.DEGreport
<https://doi.org/10.18129/B9.bioc.DEGreport>, R package version 1.36.0,
<https://bioconductor.org/packages/DEGreport>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {DEGreport: Report of DEG analysis},
author = {Lorena Pantano},
year = {2023},
note = {R package version 1.36.0},
url = {https://bioconductor.org/packages/DEGreport},
doi = {10.18129/B9.bioc.DEGreport},
}
citation("AnnotationHub")
To cite package 'AnnotationHub' in publications use:
Morgan M, Shepherd L (2023). _AnnotationHub: Client to access AnnotationHub resources_.
doi:10.18129/B9.bioc.AnnotationHub <https://doi.org/10.18129/B9.bioc.AnnotationHub>, R package version
3.8.0, <https://bioconductor.org/packages/AnnotationHub>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {AnnotationHub: Client to access AnnotationHub resources},
author = {Martin Morgan and Lori Shepherd},
year = {2023},
note = {R package version 3.8.0},
url = {https://bioconductor.org/packages/AnnotationHub},
doi = {10.18129/B9.bioc.AnnotationHub},
}
citation("ensembldb")
To cite the ensembldb package please use:
Rainer J, Gatto L, Weichenberger CX (2019) ensembldb: an R package to create and use Ensembl-based
annotation resources. Bioinformatics. doi:10.1093/bioinformatics/btz031
A BibTeX entry for LaTeX users is
@Article{,
title = {ensembldb: an R package to create and use Ensembl-based annotation resources.},
author = {Johannes Rainer and Laurent Gatto and Christian X. Weichenberger},
year = {2019},
journal = {Bioinformatics},
doi = {10.1093/bioinformatics/btz031},
url = {https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btz031/5301311},
}
citation("annotables")
To cite package 'annotables' in publications use:
Turner S (2023). _annotables: Ensembl Annotation Tables_. R package version 0.2.0,
<https://github.com/stephenturner/annotables>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {annotables: Ensembl Annotation Tables},
author = {Stephen Turner},
year = {2023},
note = {R package version 0.2.0},
url = {https://github.com/stephenturner/annotables},
}
citation("pals")
To cite package 'pals' in publications use:
Wright K (2021). _pals: Color Palettes, Colormaps, and Tools to Evaluate Them_. R package version 1.7,
<https://CRAN.R-project.org/package=pals>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {pals: Color Palettes, Colormaps, and Tools to Evaluate Them},
author = {Kevin Wright},
year = {2021},
note = {R package version 1.7},
url = {https://CRAN.R-project.org/package=pals},
}
citation("patchwork")
To cite package 'patchwork' in publications use:
Pedersen T (2022). _patchwork: The Composer of Plots_. R package version 1.1.2,
<https://CRAN.R-project.org/package=patchwork>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {patchwork: The Composer of Plots},
author = {Thomas Lin Pedersen},
year = {2022},
note = {R package version 1.1.2},
url = {https://CRAN.R-project.org/package=patchwork},
}
citation("ggpubr")
To cite package 'ggpubr' in publications use:
Kassambara A (2023). _ggpubr: 'ggplot2' Based Publication Ready Plots_. R package version 0.6.0,
<https://CRAN.R-project.org/package=ggpubr>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {ggpubr: 'ggplot2' Based Publication Ready Plots},
author = {Alboukadel Kassambara},
year = {2023},
note = {R package version 0.6.0},
url = {https://CRAN.R-project.org/package=ggpubr},
}
citation("scales")
To cite package 'scales' in publications use:
Wickham H, Seidel D (2022). _scales: Scale Functions for Visualization_. R package version 1.2.1,
<https://CRAN.R-project.org/package=scales>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {scales: Scale Functions for Visualization},
author = {Hadley Wickham and Dana Seidel},
year = {2022},
note = {R package version 1.2.1},
url = {https://CRAN.R-project.org/package=scales},
}
citation("ggiraphExtra")
To cite package 'ggiraphExtra' in publications use:
Moon K (2020). _ggiraphExtra: Make Interactive 'ggplot2'. Extension to 'ggplot2' and 'ggiraph'_. R
package version 0.3.0, <https://CRAN.R-project.org/package=ggiraphExtra>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {ggiraphExtra: Make Interactive 'ggplot2'. Extension to 'ggplot2' and 'ggiraph'},
author = {Keon-Woong Moon},
year = {2020},
note = {R package version 0.3.0},
url = {https://CRAN.R-project.org/package=ggiraphExtra},
}
citation("flextable")
To cite package 'flextable' in publications use:
Gohel D, Skintzos P (2023). _flextable: Functions for Tabular Reporting_. R package version 0.9.2,
<https://CRAN.R-project.org/package=flextable>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {flextable: Functions for Tabular Reporting},
author = {David Gohel and Panagiotis Skintzos},
year = {2023},
note = {R package version 0.9.2},
url = {https://CRAN.R-project.org/package=flextable},
}
citation("pals")
To cite package 'pals' in publications use:
Wright K (2021). _pals: Color Palettes, Colormaps, and Tools to Evaluate Them_. R package version 1.7,
<https://CRAN.R-project.org/package=pals>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {pals: Color Palettes, Colormaps, and Tools to Evaluate Them},
author = {Kevin Wright},
year = {2021},
note = {R package version 1.7},
url = {https://CRAN.R-project.org/package=pals},
}
citation("patchwork")
To cite package 'patchwork' in publications use:
Pedersen T (2022). _patchwork: The Composer of Plots_. R package version 1.1.2,
<https://CRAN.R-project.org/package=patchwork>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {patchwork: The Composer of Plots},
author = {Thomas Lin Pedersen},
year = {2022},
note = {R package version 1.1.2},
url = {https://CRAN.R-project.org/package=patchwork},
}
citation("ggiraph")
To cite package 'ggiraph' in publications use:
Gohel D, Skintzos P (2023). _ggiraph: Make 'ggplot2' Graphics Interactive_. R package version 0.8.7,
<https://CRAN.R-project.org/package=ggiraph>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {ggiraph: Make 'ggplot2' Graphics Interactive},
author = {David Gohel and Panagiotis Skintzos},
year = {2023},
note = {R package version 0.8.7},
url = {https://CRAN.R-project.org/package=ggiraph},
}
citation("officer")
To cite package 'officer' in publications use:
Gohel D (2023). _officer: Manipulation of Microsoft Word and PowerPoint Documents_. R package version
0.6.2, <https://CRAN.R-project.org/package=officer>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {officer: Manipulation of Microsoft Word and PowerPoint Documents},
author = {David Gohel},
year = {2023},
note = {R package version 0.6.2},
url = {https://CRAN.R-project.org/package=officer},
}
citation("magrittr")
To cite package 'magrittr' in publications use:
Bache S, Wickham H (2022). _magrittr: A Forward-Pipe Operator for R_. R package version 2.0.3,
<https://CRAN.R-project.org/package=magrittr>.
A BibTeX entry for LaTeX users is
@Manual{,
title = {magrittr: A Forward-Pipe Operator for R},
author = {Stefan Milton Bache and Hadley Wickham},
year = {2022},
note = {R package version 2.0.3},
url = {https://CRAN.R-project.org/package=magrittr},
}
sessionInfo()
R version 4.3.0 (2023-04-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 8.9 (Ootpa)
Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so; LAPACK version 3.9.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C 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 LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: America/Chicago
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] magrittr_2.0.3 officer_0.6.2 ggiraph_0.8.7 flextable_0.9.2
[5] ggiraphExtra_0.3.0 scales_1.2.1 ggpubr_0.6.0 patchwork_1.1.2
[9] pals_1.7 annotables_0.2.0 ensembldb_2.24.0 AnnotationFilter_1.24.0
[13] GenomicFeatures_1.52.0 AnnotationDbi_1.62.1 AnnotationHub_3.8.0 BiocFileCache_2.8.0
[17] dbplyr_2.3.2 DEGreport_1.36.0 ggrepel_0.9.3 apeglm_1.22.1
[21] pheatmap_1.0.12 limma_3.56.1 GeneStructureTools_1.20.0 lubridate_1.9.2
[25] forcats_1.0.0 stringr_1.5.0 dplyr_1.1.2 purrr_1.0.1
[29] readr_2.1.4 tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.2
[33] tidyverse_2.0.0 DESeq2_1.40.1 SummarizedExperiment_1.30.1 Biobase_2.60.0
[37] MatrixGenerics_1.12.0 matrixStats_1.0.0 GenomicRanges_1.52.0 GenomeInfoDb_1.36.0
[41] IRanges_2.34.0 S4Vectors_0.38.1 BiocGenerics_0.46.0
loaded via a namespace (and not attached):
[1] ProtGenerics_1.32.0 bitops_1.0-7 insight_0.19.3
[4] httr_1.4.6 RColorBrewer_1.1-3 doParallel_1.0.17
[7] numDeriv_2016.8-1.1 tools_4.3.0 backports_1.4.1
[10] sjlabelled_1.2.0 utf8_1.2.3 R6_2.5.1
[13] lazyeval_0.2.2 mgcv_1.8-42 Gviz_1.44.0
[16] GetoptLong_1.0.5 withr_2.5.0 prettyunits_1.1.1
[19] gridExtra_2.3 textshaping_0.3.6 cli_3.6.1
[22] logging_0.10-108 mvtnorm_1.1-3 askpass_1.1
[25] Rsamtools_2.16.0 systemfonts_1.0.4 foreign_0.8-84
[28] gfonts_0.2.0 stringdist_0.9.10 dichromat_2.0-0.1
[31] BSgenome_1.68.0 bbmle_1.0.25 maps_3.4.1
[34] httpcode_0.3.0 rstudioapi_0.14 RSQLite_2.3.1
[37] generics_0.1.3 shape_1.4.6 BiocIO_1.10.0
[40] zip_2.3.0 car_3.1-2 Matrix_1.5-4.1
[43] interp_1.1-4 fansi_1.0.4 abind_1.4-5
[46] lifecycle_1.0.3 yaml_2.3.7 edgeR_3.42.4
[49] carData_3.0-5 grid_4.3.0 blob_1.2.4
[52] promises_1.2.0.1 crayon_1.5.2 ppcor_1.1
[55] bdsmatrix_1.3-6 lattice_0.21-8 cowplot_1.1.1
[58] KEGGREST_1.40.0 mapproj_1.2.11 pillar_1.9.0
[61] knitr_1.43 ComplexHeatmap_2.16.0 rjson_0.2.21
[64] codetools_0.2-19 glue_1.6.2 fontLiberation_0.1.0
[67] data.table_1.14.8 vctrs_0.6.3 png_0.1-8
[70] gtable_0.3.3 emdbook_1.3.12 cachem_1.0.8
[73] xfun_0.39 S4Arrays_1.0.4 mime_0.12
[76] ConsensusClusterPlus_1.64.0 coda_0.19-4 iterators_1.0.14
[79] interactiveDisplayBase_1.38.0 ellipsis_0.3.2 nlme_3.1-162
[82] fontquiver_0.2.1 bit64_4.0.5 progress_1.2.2
[85] filelock_1.0.2 rpart_4.1.19 colorspace_2.1-0
[88] DBI_1.1.3 Hmisc_5.1-0 nnet_7.3-19
[91] mnormt_2.1.1 tidyselect_1.2.0 bit_4.0.5
[94] compiler_4.3.0 curl_5.0.0 htmlTable_2.4.1
[97] BSgenome.Mmusculus.UCSC.mm10_1.4.3 xml2_1.3.4 fontBitstreamVera_0.1.1
[100] ggdendro_0.1.23 DelayedArray_0.26.3 rtracklayer_1.60.0
[103] checkmate_2.2.0 psych_2.3.3 rappdirs_0.3.3
[106] digest_0.6.32 rmarkdown_2.22 XVector_0.40.0
[109] htmltools_0.5.5 pkgconfig_2.0.3 jpeg_0.1-10
[112] base64enc_0.1-3 fastmap_1.1.1 rlang_1.1.1
[115] GlobalOptions_0.1.2 htmlwidgets_1.6.2 shiny_1.7.4
[118] jsonlite_1.8.5 BiocParallel_1.34.2 VariantAnnotation_1.46.0
[121] RCurl_1.98-1.12 Formula_1.2-5 GenomeInfoDbData_1.2.10
[124] munsell_0.5.0 Rcpp_1.0.10 gdtools_0.3.3
[127] stringi_1.7.12 zlibbioc_1.46.0 MASS_7.3-60
[130] plyr_1.8.8 parallel_4.3.0 sjmisc_2.8.9
[133] deldir_1.0-9 Biostrings_2.68.1 splines_4.3.0
[136] hms_1.1.3 circlize_0.4.15 locfit_1.5-9.7
[139] uuid_1.1-0 ggsignif_0.6.4 reshape2_1.4.4
[142] biomaRt_2.56.1 crul_1.4.0 BiocVersion_3.17.1
[145] XML_3.99-0.14 evaluate_0.21 mycor_0.1.1
[148] latticeExtra_0.6-30 biovizBase_1.48.0 BiocManager_1.30.20
[151] tzdb_0.4.0 foreach_1.5.2 httpuv_1.6.11
[154] openssl_2.0.6 reshape_0.8.9 clue_0.3-64
[157] broom_1.0.4 xtable_1.8-4 restfulr_0.0.15
[160] rstatix_0.7.2 later_1.3.1 ragg_1.2.5
[163] memoise_2.0.1 GenomicAlignments_1.36.0 cluster_2.1.4
[166] timechange_0.2.0
Reads in the counts and metadata files into R.
counts=read.table(file="PM_gencode_featurecounts.Rmatrix.txt", header=TRUE)
metadata_PB=read.csv(file="metadata_PB_untrimmed.csv", header=TRUE, row.names=1)
metadata_PL=read.csv(file="metadata_PL_untrimmed.csv", header=TRUE, row.names=1)
metadata_PM=read.csv(file="metadata_PM_untrimmed.csv", header=TRUE, row.names=1)
metadata_PS=read.csv(file="metadata_PS_untrimmed.csv", header=TRUE, row.names=1)
metadata_all=read.csv(file="metadata_untrimmed.csv", header=TRUE, row.names=1)
If it is necessary to remove the version numbers from the Ensembl IDs, the following code can be run. Depending on the sample set, this may result in duplicate rows. If so, the data set needs to be de-duplicated as well.
counts$Geneid=removeVersion(counts$Geneid)
counts=plyr::ddply(counts, "Geneid", plyr::numcolwise(sum)) %>%
remove_rownames %>%
column_to_rownames(var="Geneid")
Subsets the counts object based on tissue types.
counts_PB=counts[,grepl("PB_RN_BA", colnames(counts))]
counts_PL=counts[,grepl("PL_RN_BA", colnames(counts))]
counts_PM=counts[,grepl("PM_RN_BA", colnames(counts))]
counts_PS=counts[,grepl("PS_RN_BA", colnames(counts))]
Sets the column names of the counts matrix to the row names of the metadata matrix and then confirms they are the same.
# PB
colnames(counts_PB)=rownames(metadata_PB)
all(colnames(counts_PB) %in% rownames(metadata_PB))
[1] TRUE
all(colnames(counts_PB)==rownames(metadata_PB))
[1] TRUE
# PL
colnames(counts_PL)=rownames(metadata_PL)
all(colnames(counts_PL) %in% rownames(metadata_PL))
[1] TRUE
all(colnames(counts_PL)==rownames(metadata_PL))
[1] TRUE
# PM
colnames(counts_PM)=rownames(metadata_PM)
all(colnames(counts_PM) %in% rownames(metadata_PM))
[1] TRUE
all(colnames(counts_PM)==rownames(metadata_PM))
[1] TRUE
# PS
colnames(counts_PS)=rownames(metadata_PS)
all(colnames(counts_PS) %in% rownames(metadata_PS))
[1] TRUE
all(colnames(counts_PS)==rownames(metadata_PS))
[1] TRUE
Generates the DESeqDataSet (dds) objects. The design parameter must correspond with the condition column name in the metadata file. This is the variable(s) by which the differential analysis will be performed. This command will generate a warning that some variables are characters; this is normal and is expected. Ensures that the reference level is set to the control variable.
# PB
dds_PB=DESeqDataSetFromMatrix(countData=counts_PB, colData=metadata_PB, design=~ THC)
Warning: some variables in design formula are characters, converting to factors
dds_PB$group=relevel(dds_PB$THC, ref="negative")
# PL
dds_PL=DESeqDataSetFromMatrix(countData=counts_PL, colData=metadata_PL, design=~ THC)
Warning: some variables in design formula are characters, converting to factors
dds_PL$group=relevel(dds_PL$THC, ref="negative")
# PM
dds_PM=DESeqDataSetFromMatrix(countData=counts_PM, colData=metadata_PM, design=~ THC)
Warning: some variables in design formula are characters, converting to factors
dds_PM$group=relevel(dds_PM$THC, ref="negative")
# PS
dds_PS=DESeqDataSetFromMatrix(countData=counts_PS, colData=metadata_PS, design=~ THC)
Warning: some variables in design formula are characters, converting to factors
dds_PS$group=relevel(dds_PS$THC, ref="negative")
Performs normalization on the dds objects and then writes the normalized counts to a txt file.
# PB
dds_PB=estimateSizeFactors(dds_PB)
sizeFactors(dds_PB)
PM_01_PB_RN_BA_220603 PM_03_PB_RN_BA_220607 PM_04_PB_RN_BA_220607 PM_05_PB_RN_BA_220606 PM_06_PB_RN_BA_220606
1.9927487 0.6533585 1.5179348 0.0765266 1.2907222
PM_07_PB_RN_BA_220608 PM_08_PB_RN_BA_220608 PM_09_PB_RN_BA_220608 PM_10_PB_RN_BA_220608 PM_11_PB_RN_BA_220608
1.3418045 0.6808323 1.4021222 0.7180940 1.7293236
PM_12_PB_RN_BA_220608 PM_13_PB_RN_BA_220608 PM_14_PB_RN_BA_220606 PM_15_PB_RN_BA_220606 PM_16_PB_RN_BA_220607
0.3475635 1.4749431 1.4372030 1.2532442 1.4619712
PM_17_PB_RN_BA_220718 PM_18_PB_RN_BA_220608 PM_19_PB_RN_BA_220607 PM_20_PB_RN_BA_220607 PM_22_PB_RN_BA_220603
0.6044018 1.9538787 1.9826549 1.0651692 1.1004654
PM_23_PB_RN_BA_220606 PM_24_PB_RN_BA_220606 PM_25_PB_RN_BA_220607 PM_26_PB_RN_BA_220608 PM_27_PB_RN_BA_220603
1.0663400 1.3711046 1.0738322 1.3351838 1.0590360
PM_28_PB_RN_BA_220607 PM_29_PB_RN_BA_220608 PM_30_PB_RN_BA_220603 PM_31_PB_RN_BA_220603 PM_32_PB_RN_BA_220603
1.4048117 0.2811831 0.2582569 0.7841821 1.0484875
PM_35_PB_RN_BA_220606 PM_36_PB_RN_BA_220603 PM_37_PB_RN_BA_220606 PM_38_PB_RN_BA_220609 PM_39_PB_RN_BA_220609
0.1375362 1.5660193 1.6073976 1.5682022 1.8185940
PM_40_PB_RN_BA_220603 PM_41_PB_RN_BA_220609 PM_42_PB_RN_BA_220609 PM_43_PB_RN_BA_220609 PM_44_PB_RN_BA_220609
1.4885983 0.3197314 1.1066387 0.8736601 1.4449211
PM_45_PB_RN_BA_220609 PM_46_PB_RN_BA_220609 PM_47_PB_RN_BA_220609 PM_48_PB_RN_BA_220718 PM_49_PB_RN_BA_220718
1.4280626 1.7844170 1.0444030 1.7825095 0.6698575
PM_50_PB_RN_BA_220718 PM_52_PB_RN_BA_220916 PM_53_PB_RN_BA_220916 PM_54_PB_RN_BA_220916 PM_55_PB_RN_BA_220916
1.5285517 1.3123926 0.5113159 0.8603980 0.7516655
PM_56_PB_RN_BA_220916 PM_57_PB_RN_BA_220916 PM_58_PB_RN_BA_220916
1.5282431 1.7239848 1.6094768
normalized_counts_PB=counts(dds_PB, normalized=TRUE)
write.table(normalized_counts_PB, file="gencode_Output_Files/PB_normalizedcounts.txt",
sep="\t", quote=FALSE, col.names=NA)
# PL
dds_PL=estimateSizeFactors(dds_PL)
sizeFactors(dds_PL)
PM_01_PL_RN_BA_220622 PM_03_PL_RN_BA_220622 PM_04_PL_RN_BA_220622 PM_06_PL_RN_BA_220622 PM_07_PL_RN_BA_220622
1.40297862 0.46367867 2.85696861 0.05537568 2.80900092
PM_08_PL_RN_BA_220622 PM_09_PL_RN_BA_220622 PM_10_PL_RN_BA_220622 PM_11_PL_RN_BA_220622 PM_12_PL_RN_BA_220622
1.58957104 2.00337729 1.17858581 2.09246939 0.36782415
PM_13_PL_RN_BA_220622 PM_14_PL_RN_BA_220720 PM_15_PL_RN_BA_220622 PM_16_PL_RN_BA_220622 PM_17_PL_RN_BA_220622
2.66309782 3.69152590 1.46677811 2.49919006 0.35155524
PM_18_PL_RN_BA_220622 PM_19_PL_RN_BA_220622 PM_20_PL_RN_BA_220622 PM_21_PL_RN_BA_220622 PM_22_PL_RN_BA_220720
0.05415938 2.05344480 0.03270570 2.45919055 0.02814269
PM_24_PL_RN_BA_220622 PM_25_PL_RN_BA_220622 PM_26_PL_RN_BA_220623 PM_27_PL_RN_BA_220623 PM_28_PL_RN_BA_220623
0.09620581 0.86529255 2.24341313 2.02138616 1.50732385
PM_29_PL_RN_BA_220623 PM_30_PL_RN_BA_220720 PM_31_PL_RN_BA_220623 PM_32_PL_RN_BA_220720 PM_33_PL_RN_BA_220623
0.41538316 1.85433641 0.68014153 1.32914721 0.56766953
PM_36_PL_RN_BA_220623 PM_37_PL_RN_BA_220623 PM_38_PL_RN_BA_220623 PM_39_PL_RN_BA_220720 PM_40_PL_RN_BA_220624
1.33532765 2.62394196 3.02302763 0.35657666 3.11030096
PM_41_PL_RN_BA_220624 PM_42_PL_RN_BA_220624 PM_43_PL_RN_BA_220720 PM_44_PL_RN_BA_220624 PM_45_PL_RN_BA_220624
1.05301075 1.51788168 0.30244477 0.85706746 2.18470276
PM_46_PL_RN_BA_220624 PM_47_PL_RN_BA_220624 PM_48_PL_RN_BA_220720 PM_49_PL_RN_BA_220720 PM_50_PL_RN_BA_220720
1.13637430 1.40400835 1.60445764 2.78561615 2.15553491
PM_51_PL_RN_BA_220920 PM_52_PL_RN_BA_220920 PM_54_PL_RN_BA_220920 PM_55_PL_RN_BA_220920 PM_56_PL_RN_BA_220923
2.82774982 0.89931011 1.50238908 0.28197242 1.50707690
PM_57_PL_RN_BA_220920 PM_58_PL_RN_BA_220923
2.06168949 3.80973657
normalized_counts_PL=counts(dds_PL, normalized=TRUE)
write.table(normalized_counts_PL, file="gencode_Output_Files/PL_normalizedcounts.txt",
sep="\t", quote=FALSE, col.names=NA)
# PM
dds_PM=estimateSizeFactors(dds_PM)
sizeFactors(dds_PM)
PM_01_PM_RN_BA_220613 PM_03_PM_RN_BA_220613 PM_04_PM_RN_BA_220613 PM_06_PM_RN_BA_220613 PM_07_PM_RN_BA_220613
1.5088832 1.4255698 1.6204817 1.1970366 0.9800432
PM_08_PM_RN_BA_220613 PM_09_PM_RN_BA_220613 PM_10_PM_RN_BA_220613 PM_11_PM_RN_BA_220613 PM_12_PM_RN_BA_220613
1.2978390 1.4868746 1.1757672 1.1081042 1.1678733
PM_13_PM_RN_BA_220613 PM_14_PM_RN_BA_220614 PM_15_PM_RN_BA_220614 PM_16_PM_RN_BA_220614 PM_17_PM_RN_BA_220719
1.1631088 1.8022065 1.0390151 1.2341993 0.6961912
PM_18_PM_RN_BA_220614 PM_19_PM_RN_BA_220614 PM_20_PM_RN_BA_220614 PM_21_PM_RN_BA_220614 PM_22_PM_RN_BA_220614
0.9084388 1.4110597 0.2821348 1.3634059 0.7242727
PM_23_PM_RN_BA_220719 PM_24_PM_RN_BA_220614 PM_25_PM_RN_BA_220614 PM_26_PM_RN_BA_220615 PM_27_PM_RN_BA_220615
0.4507395 1.3024728 1.3538257 1.4185118 1.2038617
PM_28_PM_RN_BA_220615 PM_29_PM_RN_BA_220615 PM_30_PM_RN_BA_220615 PM_31_PM_RN_BA_220615 PM_33_PM_RN_BA_220615
1.3040438 0.6821155 1.1197386 1.5807151 0.1739217
PM_36_PM_RN_BA_220615 PM_37_PM_RN_BA_220615 PM_38_PM_RN_BA_220615 PM_39_PM_RN_BA_220615 PM_40_PM_RN_BA_220616
1.6355448 0.4399030 1.1170166 1.2976556 0.1718784
PM_41_PM_RN_BA_220616 PM_42_PM_RN_BA_220616 PM_43_PM_RN_BA_220616 PM_44_PM_RN_BA_220616 PM_45_PM_RN_BA_220616
0.8429614 0.7767747 0.9374023 1.2488331 1.3410268
PM_46_PM_RN_BA_220616 PM_47_PM_RN_BA_220616 PM_48_PM_RN_BA_220719 PM_49_PM_RN_BA_220719 PM_50_PM_RN_BA_220719
1.1273183 1.1636958 1.4509711 1.2990044 0.5381660
PM_51_PM_RN_BA_220920 PM_52_PM_RN_BA_220920 PM_53_PM_RN_BA_220920 PM_54_PM_RN_BA_220920 PM_55_PM_RN_BA_220920
1.4074151 1.0996795 2.1120325 1.5010878 0.1360339
PM_56_PM_RN_BA_220920 PM_57_PM_RN_BA_220920 PM_58_PM_RN_BA_220920
1.6941223 1.5719031 1.5294640
normalized_counts_PM=counts(dds_PM, normalized=TRUE)
write.table(normalized_counts_PM, file="gencode_Output_Files/PM_normalizedcounts.txt",
sep="\t", quote=FALSE, col.names=NA)
# PS
dds_PS=estimateSizeFactors(dds_PS)
sizeFactors(dds_PS)
PM_01_PS_RN_BA_220627 PM_02_PS_RN_BA_220627 PM_03_PS_RN_BA_220627 PM_04_PS_RN_BA_220627 PM_06_PS_RN_BA_220627
1.058600541 0.658393088 1.222699830 1.834500045 4.842249194
PM_07_PS_RN_BA_220627 PM_08_PS_RN_BA_220627 PM_11_PS_RN_BA_220627 PM_13_PS_RN_BA_220627 PM_14_PS_RN_BA_220628
8.126935887 0.398138600 7.237634747 4.384803112 6.038798058
PM_15_PS_RN_BA_220714 PM_16_PS_RN_BA_220628 PM_17_PS_RN_BA_220628 PM_18_PS_RN_BA_220628 PM_22_PS_RN_BA_220706
8.271116876 6.656457794 8.197012251 0.642441448 0.067976354
PM_23_PS_RN_BA_220628 PM_24_PS_RN_BA_220628 PM_25_PS_RN_BA_220628 PM_27_PS_RN_BA_220628 PM_29_PS_RN_BA_220628
12.952010726 0.208645715 0.090946482 0.643588097 0.004792714
PM_30_PS_RN_BA_220706 PM_31_PS_RN_BA_220630 PM_33_PS_RN_BA_220630 PM_35_PS_RN_BA_220630 PM_36_PS_RN_BA_220630
1.432265234 0.110024402 0.043707832 7.299222374 3.312055627
PM_37_PS_RN_BA_220630 PM_38_PS_RN_BA_220630 PM_39_PS_RN_BA_220630 PM_40_PS_RN_BA_220630 PM_42_PS_RN_BA_220630
1.806424448 1.457886987 1.219826843 0.076546734 0.052348316
PM_43_PS_RN_BA_220630 PM_44_PS_RN_BA_220630 PM_47_PS_RN_BA_220706 PM_48_PS_RN_BA_220714 PM_49_PS_RN_BA_220714
0.272828739 6.442905955 1.594276066 0.361234164 0.550746554
PM_50_PS_RN_BA_220714 PM_52_PS_RN_BA_220921 PM_53_PS_RN_BA_220921 PM_54_PS_RN_BA_220921 PM_55_PS_RN_BA_220921
5.975779874 2.006948025 0.479613646 3.565880650 0.126282891
PM_56_PS_RN_BA_220921 PM_57_PS_RN_BA_220921 PM_58_PS_RN_BA_220921
4.111315224 0.544084792 5.130149915
normalized_counts_PS=counts(dds_PS, normalized=TRUE)
write.table(normalized_counts_PS, file="gencode_Output_Files/PS_normalizedcounts.txt",
sep="\t", quote=FALSE, col.names=NA)
Transforms counts for data visualization.
rld_PB=rlog(dds_PB, blind=TRUE)
rlog() may take a long time with 50 or more samples,
vst() is a much faster transformation
rld_PL=rlog(dds_PL, blind=TRUE)
rlog() may take a long time with 50 or more samples,
vst() is a much faster transformation
rld_PM=rlog(dds_PM, blind=TRUE)
rlog() may take a long time with 50 or more samples,
vst() is a much faster transformation
rld_PS=rlog(dds_PS, blind=TRUE)
rlog() may take a few minutes with 30 or more samples,
vst() is a much faster transformation
Plots the PCA.
# PB
pca1=plotPCA(rld_PB, intgroup=c("THC"))+
geom_point(alpha=0.75, color="black", pch=21, size=3, show.legend=TRUE)+
# ggrepel::geom_text_repel(aes(label=name), color="black", force=3, min.segment.length=0.5, box.padding=0.5,
# max.overlaps=10, show.legend=FALSE)+
scale_color_manual(limits=c("positive", "negative"),
labels=c("Positive", "Negative"),
values=as.vector(polychrome(2)))+
theme_minimal()+
theme(axis.text.y=element_text(color="black", face="bold", size=rel(1.25)),
axis.text.x=element_text(color="black", face="bold", size=rel(1.25)),
axis.title.x=element_text(color="black", face="bold", size=rel(1.1)),
axis.title.y=element_text(color="black", face="bold", size=rel(1.1)),
title=element_text(color="black", face="bold", size=rel(0.8)),
panel.border=element_rect(color="black", linewidth=1, fill=NA),
strip.background=element_rect(color="black", linewidth=1),
strip.text=element_text(color="black", face="bold", size=rel(0.8)),
axis.line=element_blank(),
axis.ticks=element_line(linewidth=1),
# legend.title=element_text(color="black", face="bold", size=rel(1)),
legend.title=element_blank(),
legend.text=element_text(color="black", face="bold", size=rel(0.8)),
legend.position="bottom",
# plot.title=element_text(color="black", face="bold", size=rel(1)))+
plot.title=element_blank())+
labs(title="Postmortem Brain", color="THC")
# PL
pca2=plotPCA(rld_PL, intgroup=c("THC"))+
geom_point(alpha=0.75, color="black", pch=21, size=3, show.legend=TRUE)+
# ggrepel::geom_text_repel(aes(label=name), color="black", force=3, min.segment.length=0.5, box.padding=0.5,
# max.overlaps=10, show.legend=FALSE)+
scale_color_manual(limits=c("positive", "negative"),
labels=c("Positive", "Negative"),
values=as.vector(polychrome(2)))+
theme_minimal()+
theme(axis.text.y=element_text(color="black", face="bold", size=rel(1.25)),
axis.text.x=element_text(color="black", face="bold", size=rel(1.25)),
axis.title.x=element_text(color="black", face="bold", size=rel(1.1)),
axis.title.y=element_text(color="black", face="bold", size=rel(1.1)),
title=element_text(color="black", face="bold", size=rel(0.8)),
panel.border=element_rect(color="black", linewidth=1, fill=NA),
strip.background=element_rect(color="black", linewidth=1),
strip.text=element_text(color="black", face="bold", size=rel(0.8)),
axis.line=element_blank(),
axis.ticks=element_line(linewidth=1),
# legend.title=element_text(color="black", face="bold", size=rel(1)),
legend.title=element_blank(),
legend.text=element_text(color="black", face="bold", size=rel(0.8)),
legend.position="bottom",
# plot.title=element_text(color="black", face="bold", size=rel(1)))+
plot.title=element_blank())+
labs(title="Postmortem Lung", color="THC")
# PM
pca3=plotPCA(rld_PM, intgroup=c("THC"))+
geom_point(alpha=0.75, color="black", pch=21, size=3, show.legend=TRUE)+
# ggrepel::geom_text_repel(aes(label=name), color="black", force=3, min.segment.length=0.5, box.padding=0.5,
# max.overlaps=10, show.legend=FALSE)+
scale_color_manual(limits=c("positive", "negative"),
labels=c("Positive", "Negative"),
values=as.vector(polychrome(2)))+
theme_minimal()+
theme(axis.text.y=element_text(color="black", face="bold", size=rel(1.25)),
axis.text.x=element_text(color="black", face="bold", size=rel(1.25)),
axis.title.x=element_text(color="black", face="bold", size=rel(1.1)),
axis.title.y=element_text(color="black", face="bold", size=rel(1.1)),
title=element_text(color="black", face="bold", size=rel(0.8)),
panel.border=element_rect(color="black", linewidth=1, fill=NA),
strip.background=element_rect(color="black", linewidth=1),
strip.text=element_text(color="black", face="bold", size=rel(0.8)),
axis.line=element_blank(),
axis.ticks=element_line(linewidth=1),
# legend.title=element_text(color="black", face="bold", size=rel(1)),
legend.title=element_blank(),
legend.text=element_text(color="black", face="bold", size=rel(0.8)),
legend.position="bottom",
# plot.title=element_text(color="black", face="bold", size=rel(1)))+
plot.title=element_blank())+
labs(title="Postmortem Muscle", color="THC")
# PS
pca4=plotPCA(rld_PS, intgroup=c("THC"))+
geom_point(alpha=0.75, color="black", pch=21, size=3, show.legend=TRUE)+
# ggrepel::geom_text_repel(aes(label=name), color="black", force=3, min.segment.length=0.5, box.padding=0.5,
# max.overlaps=10, show.legend=FALSE)+
scale_color_manual(limits=c("positive", "negative"),
labels=c("Positive", "Negative"),
values=as.vector(polychrome(2)))+
theme_minimal()+
theme(axis.text.y=element_text(color="black", face="bold", size=rel(1.25)),
axis.text.x=element_text(color="black", face="bold", size=rel(1.25)),
axis.title.x=element_text(color="black", face="bold", size=rel(1.1)),
axis.title.y=element_text(color="black", face="bold", size=rel(1.1)),
title=element_text(color="black", face="bold", size=rel(0.8)),
panel.border=element_rect(color="black", linewidth=1, fill=NA),
strip.background=element_rect(color="black", linewidth=1),
strip.text=element_text(color="black", face="bold", size=rel(0.8)),
axis.line=element_blank(),
axis.ticks=element_line(linewidth=1),
# legend.title=element_text(color="black", face="bold", size=rel(1)),
legend.title=element_blank(),
legend.text=element_text(color="black", face="bold", size=rel(0.8)),
legend.position="bottom",
# plot.title=element_text(color="black", face="bold", size=rel(1)))+
plot.title=element_blank())+
labs(title="Postmortem Blood", color="THC")
Creating figures of the PCAs.
pca_patchwork=(pca1 + pca2 + pca3 + pca4 + plot_layout(ncol=2, guides="collect") & theme(legend.position="bottom"))
pca_patchwork + plot_annotation(tag_levels="A")
ggsave("PM_pca_plots.png", width=7, height=8.36, unit="in", dpi=320)
Extracts the rlog matrix from rld, computes pairwise correlation values for the samples and then checks its outputs, and then plots the heatmaps.
# PB
rld_mat_PB=assay(rld_PB)
rld_cor_PB=cor(rld_mat_PB)
head(rld_cor_PB)
PM_01_PB_RN_BA_220603 PM_03_PB_RN_BA_220607 PM_04_PB_RN_BA_220607 PM_05_PB_RN_BA_220606
PM_01_PB_RN_BA_220603 1.0000000 0.9901091 0.9965209 0.9755158
PM_03_PB_RN_BA_220607 0.9901091 1.0000000 0.9917846 0.9802242
PM_04_PB_RN_BA_220607 0.9965209 0.9917846 1.0000000 0.9756877
PM_05_PB_RN_BA_220606 0.9755158 0.9802242 0.9756877 1.0000000
PM_06_PB_RN_BA_220606 0.9855471 0.9855993 0.9858795 0.9851248
PM_07_PB_RN_BA_220608 0.9863301 0.9839060 0.9872548 0.9746333
PM_06_PB_RN_BA_220606 PM_07_PB_RN_BA_220608 PM_08_PB_RN_BA_220608 PM_09_PB_RN_BA_220608
PM_01_PB_RN_BA_220603 0.9855471 0.9863301 0.9785581 0.9842733
PM_03_PB_RN_BA_220607 0.9855993 0.9839060 0.9806932 0.9809368
PM_04_PB_RN_BA_220607 0.9858795 0.9872548 0.9784769 0.9841310
PM_05_PB_RN_BA_220606 0.9851248 0.9746333 0.9861917 0.9817424
PM_06_PB_RN_BA_220606 1.0000000 0.9830799 0.9933317 0.9933161
PM_07_PB_RN_BA_220608 0.9830799 1.0000000 0.9786033 0.9830570
PM_10_PB_RN_BA_220608 PM_11_PB_RN_BA_220608 PM_12_PB_RN_BA_220608 PM_13_PB_RN_BA_220608
PM_01_PB_RN_BA_220603 0.9717560 0.9929799 0.9740851 0.9755814
PM_03_PB_RN_BA_220607 0.9716203 0.9877791 0.9813001 0.9723104
PM_04_PB_RN_BA_220607 0.9728802 0.9929132 0.9746039 0.9750949
PM_05_PB_RN_BA_220606 0.9708778 0.9807026 0.9909002 0.9784128
PM_06_PB_RN_BA_220606 0.9801115 0.9932426 0.9876757 0.9890929
PM_07_PB_RN_BA_220608 0.9701029 0.9873095 0.9738313 0.9762802
PM_14_PB_RN_BA_220606 PM_15_PB_RN_BA_220606 PM_16_PB_RN_BA_220607 PM_17_PB_RN_BA_220718
PM_01_PB_RN_BA_220603 0.9775747 0.9747433 0.9765312 0.9843602
PM_03_PB_RN_BA_220607 0.9747705 0.9706778 0.9725660 0.9864632
PM_04_PB_RN_BA_220607 0.9769762 0.9741001 0.9759970 0.9843375
PM_05_PB_RN_BA_220606 0.9805473 0.9778020 0.9787977 0.9840879
PM_06_PB_RN_BA_220606 0.9937848 0.9915884 0.9921455 0.9923032
PM_07_PB_RN_BA_220608 0.9775195 0.9754976 0.9770170 0.9808634
PM_18_PB_RN_BA_220608 PM_19_PB_RN_BA_220607 PM_20_PB_RN_BA_220607 PM_22_PB_RN_BA_220603
PM_01_PB_RN_BA_220603 0.9957606 0.9871995 0.9849784 0.9821417
PM_03_PB_RN_BA_220607 0.9886438 0.9843303 0.9845239 0.9793116
PM_04_PB_RN_BA_220607 0.9943797 0.9860169 0.9846571 0.9818114
PM_05_PB_RN_BA_220606 0.9777466 0.9806215 0.9832981 0.9793994
PM_06_PB_RN_BA_220606 0.9891921 0.9905778 0.9931690 0.9939133
PM_07_PB_RN_BA_220608 0.9868246 0.9834152 0.9819857 0.9802699
PM_23_PB_RN_BA_220606 PM_24_PB_RN_BA_220606 PM_25_PB_RN_BA_220607 PM_26_PB_RN_BA_220608
PM_01_PB_RN_BA_220603 0.9810947 0.9958426 0.9873861 0.9785448
PM_03_PB_RN_BA_220607 0.9775491 0.9908163 0.9834349 0.9754719
PM_04_PB_RN_BA_220607 0.9806034 0.9959215 0.9874594 0.9778489
PM_05_PB_RN_BA_220606 0.9804694 0.9795154 0.9825426 0.9806273
PM_06_PB_RN_BA_220606 0.9925625 0.9900548 0.9942281 0.9939817
PM_07_PB_RN_BA_220608 0.9803718 0.9876731 0.9840954 0.9785637
PM_27_PB_RN_BA_220603 PM_28_PB_RN_BA_220607 PM_29_PB_RN_BA_220608 PM_30_PB_RN_BA_220603
PM_01_PB_RN_BA_220603 0.9909775 0.9855913 0.9879817 0.9782437
PM_03_PB_RN_BA_220607 0.9889029 0.9843359 0.9893007 0.9799113
PM_04_PB_RN_BA_220607 0.9905886 0.9852213 0.9883848 0.9778960
PM_05_PB_RN_BA_220606 0.9824394 0.9820877 0.9866963 0.9858239
PM_06_PB_RN_BA_220606 0.9907489 0.9914526 0.9918916 0.9914125
PM_07_PB_RN_BA_220608 0.9855766 0.9831244 0.9849305 0.9780425
PM_31_PB_RN_BA_220603 PM_32_PB_RN_BA_220603 PM_35_PB_RN_BA_220606 PM_36_PB_RN_BA_220603
PM_01_PB_RN_BA_220603 0.9754693 0.9793567 0.9773495 0.9950925
PM_03_PB_RN_BA_220607 0.9724432 0.9778665 0.9804925 0.9900565
PM_04_PB_RN_BA_220607 0.9752859 0.9788968 0.9777990 0.9947534
PM_05_PB_RN_BA_220606 0.9790986 0.9818489 0.9839451 0.9805891
PM_06_PB_RN_BA_220606 0.9919468 0.9923454 0.9904186 0.9913215
PM_07_PB_RN_BA_220608 0.9768913 0.9796217 0.9777062 0.9875329
PM_37_PB_RN_BA_220606 PM_38_PB_RN_BA_220609 PM_39_PB_RN_BA_220609 PM_40_PB_RN_BA_220603
PM_01_PB_RN_BA_220603 0.9738679 0.9791838 0.9968656 0.9797049
PM_03_PB_RN_BA_220607 0.9711035 0.9779420 0.9909328 0.9776203
PM_04_PB_RN_BA_220607 0.9735817 0.9786025 0.9961988 0.9794711
PM_05_PB_RN_BA_220606 0.9783617 0.9818949 0.9774080 0.9823178
PM_06_PB_RN_BA_220606 0.9920999 0.9907470 0.9880266 0.9926192
PM_07_PB_RN_BA_220608 0.9744059 0.9794491 0.9865220 0.9804671
PM_41_PB_RN_BA_220609 PM_42_PB_RN_BA_220609 PM_43_PB_RN_BA_220609 PM_44_PB_RN_BA_220609
PM_01_PB_RN_BA_220603 0.9750532 0.9954317 0.9857752 0.9807207
PM_03_PB_RN_BA_220607 0.9823896 0.9889208 0.9867488 0.9787779
PM_04_PB_RN_BA_220607 0.9758833 0.9945386 0.9853340 0.9803736
PM_05_PB_RN_BA_220606 0.9883736 0.9771400 0.9850938 0.9824471
PM_06_PB_RN_BA_220606 0.9866711 0.9889500 0.9917196 0.9945050
PM_07_PB_RN_BA_220608 0.9780788 0.9857168 0.9814694 0.9807023
PM_45_PB_RN_BA_220609 PM_46_PB_RN_BA_220609 PM_47_PB_RN_BA_220609 PM_48_PB_RN_BA_220718
PM_01_PB_RN_BA_220603 0.9930457 0.9930573 0.9919979 0.9944350
PM_03_PB_RN_BA_220607 0.9873628 0.9918957 0.9871815 0.9927051
PM_04_PB_RN_BA_220607 0.9926572 0.9931280 0.9915813 0.9957659
PM_05_PB_RN_BA_220606 0.9810813 0.9812422 0.9790655 0.9779724
PM_06_PB_RN_BA_220606 0.9923422 0.9901902 0.9918409 0.9858069
PM_07_PB_RN_BA_220608 0.9869426 0.9861186 0.9854367 0.9870139
PM_49_PB_RN_BA_220718 PM_50_PB_RN_BA_220718 PM_52_PB_RN_BA_220916 PM_53_PB_RN_BA_220916
PM_01_PB_RN_BA_220603 0.9780281 0.9869637 0.9814291 0.9820308
PM_03_PB_RN_BA_220607 0.9782839 0.9854810 0.9754347 0.9811132
PM_04_PB_RN_BA_220607 0.9780459 0.9875752 0.9799466 0.9822948
PM_05_PB_RN_BA_220606 0.9848959 0.9799129 0.9729781 0.9815129
PM_06_PB_RN_BA_220606 0.9944583 0.9900513 0.9844777 0.9922553
PM_07_PB_RN_BA_220608 0.9780897 0.9880291 0.9817195 0.9863578
PM_54_PB_RN_BA_220916 PM_55_PB_RN_BA_220916 PM_56_PB_RN_BA_220916 PM_57_PB_RN_BA_220916
PM_01_PB_RN_BA_220603 0.9924716 0.9885575 0.9964062 0.9786740
PM_03_PB_RN_BA_220607 0.9866341 0.9877329 0.9901431 0.9756041
PM_04_PB_RN_BA_220607 0.9915561 0.9879424 0.9959651 0.9779961
PM_05_PB_RN_BA_220606 0.9776259 0.9825420 0.9769205 0.9808700
PM_06_PB_RN_BA_220606 0.9882380 0.9907674 0.9882553 0.9934185
PM_07_PB_RN_BA_220608 0.9848664 0.9824158 0.9866356 0.9790188
PM_58_PB_RN_BA_220916
PM_01_PB_RN_BA_220603 0.9958764
PM_03_PB_RN_BA_220607 0.9907376
PM_04_PB_RN_BA_220607 0.9961337
PM_05_PB_RN_BA_220606 0.9789579
PM_06_PB_RN_BA_220606 0.9892049
PM_07_PB_RN_BA_220608 0.9876831
pheatmap(rld_cor_PB,
# color=colorRampPalette(c("", "", ""))(200), ## The color palette.
display_numbers=matrix(ifelse(rld_cor_PB>0.99995, "*", ""), nrow(rld_cor_PB)), fontsize_row=5)
# PL
rld_mat_PL=assay(rld_PL)
rld_cor_PL=cor(rld_mat_PL)
head(rld_cor_PL)
PM_01_PL_RN_BA_220622 PM_03_PL_RN_BA_220622 PM_04_PL_RN_BA_220622 PM_06_PL_RN_BA_220622
PM_01_PL_RN_BA_220622 1.0000000 0.9910972 0.9908423 0.9788182
PM_03_PL_RN_BA_220622 0.9910972 1.0000000 0.9891178 0.9800172
PM_04_PL_RN_BA_220622 0.9908423 0.9891178 1.0000000 0.9667516
PM_06_PL_RN_BA_220622 0.9788182 0.9800172 0.9667516 1.0000000
PM_07_PL_RN_BA_220622 0.9885582 0.9876787 0.9956143 0.9652953
PM_08_PL_RN_BA_220622 0.9939659 0.9902991 0.9935189 0.9759458
PM_07_PL_RN_BA_220622 PM_08_PL_RN_BA_220622 PM_09_PL_RN_BA_220622 PM_10_PL_RN_BA_220622
PM_01_PL_RN_BA_220622 0.9885582 0.9939659 0.9895170 0.9912643
PM_03_PL_RN_BA_220622 0.9876787 0.9902991 0.9878119 0.9879389
PM_04_PL_RN_BA_220622 0.9956143 0.9935189 0.9956073 0.9867542
PM_06_PL_RN_BA_220622 0.9652953 0.9759458 0.9678336 0.9767337
PM_07_PL_RN_BA_220622 1.0000000 0.9932818 0.9952379 0.9866332
PM_08_PL_RN_BA_220622 0.9932818 1.0000000 0.9930632 0.9925234
PM_11_PL_RN_BA_220622 PM_12_PL_RN_BA_220622 PM_13_PL_RN_BA_220622 PM_14_PL_RN_BA_220720
PM_01_PL_RN_BA_220622 0.9889402 0.9907880 0.9869112 0.9764878
PM_03_PL_RN_BA_220622 0.9857152 0.9904311 0.9869310 0.9699733
PM_04_PL_RN_BA_220622 0.9951765 0.9894851 0.9961892 0.9724355
PM_06_PL_RN_BA_220622 0.9648406 0.9724260 0.9637285 0.9591717
PM_07_PL_RN_BA_220622 0.9957924 0.9884486 0.9952194 0.9700284
PM_08_PL_RN_BA_220622 0.9924971 0.9897543 0.9896300 0.9761248
PM_15_PL_RN_BA_220622 PM_16_PL_RN_BA_220622 PM_17_PL_RN_BA_220622 PM_18_PL_RN_BA_220622
PM_01_PL_RN_BA_220622 0.9832401 0.9862336 0.9864013 0.9910573
PM_03_PL_RN_BA_220622 0.9827567 0.9859469 0.9861141 0.9892271
PM_04_PL_RN_BA_220622 0.9911250 0.9962128 0.9853669 0.9857015
PM_06_PL_RN_BA_220622 0.9610647 0.9638893 0.9721627 0.9803519
PM_07_PL_RN_BA_220622 0.9914297 0.9952613 0.9865335 0.9850350
PM_08_PL_RN_BA_220622 0.9860369 0.9908062 0.9862337 0.9904943
PM_19_PL_RN_BA_220622 PM_20_PL_RN_BA_220622 PM_21_PL_RN_BA_220622 PM_22_PL_RN_BA_220720
PM_01_PL_RN_BA_220622 0.9893315 0.9895382 0.9844077 0.9802875
PM_03_PL_RN_BA_220622 0.9888412 0.9880382 0.9853940 0.9816627
PM_04_PL_RN_BA_220622 0.9971654 0.9867954 0.9937606 0.9788522
PM_06_PL_RN_BA_220622 0.9669540 0.9738614 0.9622983 0.9717783
PM_07_PL_RN_BA_220622 0.9967776 0.9854025 0.9945551 0.9774128
PM_08_PL_RN_BA_220622 0.9931623 0.9905607 0.9885249 0.9803582
PM_24_PL_RN_BA_220622 PM_25_PL_RN_BA_220622 PM_26_PL_RN_BA_220623 PM_27_PL_RN_BA_220623
PM_01_PL_RN_BA_220622 0.9910214 0.9922485 0.9928030 0.9790891
PM_03_PL_RN_BA_220622 0.9857184 0.9889739 0.9885774 0.9774138
PM_04_PL_RN_BA_220622 0.9863849 0.9889040 0.9953340 0.9835981
PM_06_PL_RN_BA_220622 0.9737984 0.9778924 0.9693387 0.9605199
PM_07_PL_RN_BA_220622 0.9872452 0.9884993 0.9939113 0.9819150
PM_08_PL_RN_BA_220622 0.9924648 0.9917105 0.9937963 0.9800863
PM_28_PL_RN_BA_220623 PM_29_PL_RN_BA_220623 PM_30_PL_RN_BA_220720 PM_31_PL_RN_BA_220623
PM_01_PL_RN_BA_220622 0.9933404 0.9894508 0.9909944 0.9904141
PM_03_PL_RN_BA_220622 0.9866288 0.9900216 0.9878813 0.9911445
PM_04_PL_RN_BA_220622 0.9885498 0.9833870 0.9929728 0.9913255
PM_06_PL_RN_BA_220622 0.9740784 0.9845542 0.9696201 0.9748131
PM_07_PL_RN_BA_220622 0.9899298 0.9812471 0.9942135 0.9896455
PM_08_PL_RN_BA_220622 0.9939345 0.9890249 0.9937233 0.9905585
PM_32_PL_RN_BA_220720 PM_33_PL_RN_BA_220623 PM_36_PL_RN_BA_220623 PM_37_PL_RN_BA_220623
PM_01_PL_RN_BA_220622 0.9899136 0.9842380 0.9895537 0.9949305
PM_03_PL_RN_BA_220622 0.9859173 0.9853551 0.9855686 0.9931103
PM_04_PL_RN_BA_220622 0.9896973 0.9765829 0.9909755 0.9904631
PM_06_PL_RN_BA_220622 0.9703771 0.9822009 0.9683041 0.9828744
PM_07_PL_RN_BA_220622 0.9879403 0.9780689 0.9934616 0.9893501
PM_08_PL_RN_BA_220622 0.9918477 0.9862346 0.9922855 0.9940937
PM_38_PL_RN_BA_220623 PM_39_PL_RN_BA_220720 PM_40_PL_RN_BA_220624 PM_41_PL_RN_BA_220624
PM_01_PL_RN_BA_220622 0.9891455 0.9898819 0.9902533 0.9908972
PM_03_PL_RN_BA_220622 0.9881962 0.9883877 0.9890439 0.9902173
PM_04_PL_RN_BA_220622 0.9966812 0.9860993 0.9975169 0.9933809
PM_06_PL_RN_BA_220622 0.9655941 0.9751445 0.9657767 0.9686879
PM_07_PL_RN_BA_220622 0.9944172 0.9864230 0.9948443 0.9937233
PM_08_PL_RN_BA_220622 0.9908344 0.9885173 0.9919258 0.9928422
PM_42_PL_RN_BA_220624 PM_43_PL_RN_BA_220720 PM_44_PL_RN_BA_220624 PM_45_PL_RN_BA_220624
PM_01_PL_RN_BA_220622 0.9940023 0.9716941 0.9855997 0.9894707
PM_03_PL_RN_BA_220622 0.9924331 0.9754358 0.9851053 0.9896770
PM_04_PL_RN_BA_220622 0.9923110 0.9591240 0.9730816 0.9973428
PM_06_PL_RN_BA_220622 0.9789120 0.9830714 0.9879308 0.9659792
PM_07_PL_RN_BA_220622 0.9884774 0.9577194 0.9712228 0.9946679
PM_08_PL_RN_BA_220622 0.9921656 0.9679535 0.9824618 0.9916550
PM_46_PL_RN_BA_220624 PM_47_PL_RN_BA_220624 PM_48_PL_RN_BA_220720 PM_49_PL_RN_BA_220720
PM_01_PL_RN_BA_220622 0.9899534 0.9907368 0.9853260 0.9915726
PM_03_PL_RN_BA_220622 0.9903483 0.9916739 0.9857574 0.9881633
PM_04_PL_RN_BA_220622 0.9825910 0.9948399 0.9934015 0.9960451
PM_06_PL_RN_BA_220622 0.9846138 0.9690746 0.9611666 0.9668113
PM_07_PL_RN_BA_220622 0.9798687 0.9925659 0.9935555 0.9934250
PM_08_PL_RN_BA_220622 0.9875212 0.9916737 0.9886173 0.9931984
PM_50_PL_RN_BA_220720 PM_51_PL_RN_BA_220920 PM_52_PL_RN_BA_220920 PM_54_PL_RN_BA_220920
PM_01_PL_RN_BA_220622 0.9937004 0.9885754 0.9864976 0.9927054
PM_03_PL_RN_BA_220622 0.9907624 0.9883850 0.9852506 0.9921812
PM_04_PL_RN_BA_220622 0.9940794 0.9962668 0.9914533 0.9873088
PM_06_PL_RN_BA_220622 0.9735155 0.9667552 0.9666389 0.9831078
PM_07_PL_RN_BA_220622 0.9925445 0.9939848 0.9939015 0.9852994
PM_08_PL_RN_BA_220622 0.9931208 0.9906672 0.9918310 0.9909328
PM_55_PL_RN_BA_220920 PM_56_PL_RN_BA_220923 PM_57_PL_RN_BA_220920 PM_58_PL_RN_BA_220923
PM_01_PL_RN_BA_220622 0.9768707 0.9939164 0.9931044 0.9874121
PM_03_PL_RN_BA_220622 0.9800790 0.9914302 0.9884081 0.9867016
PM_04_PL_RN_BA_220622 0.9668560 0.9857882 0.9929644 0.9954676
PM_06_PL_RN_BA_220622 0.9841286 0.9849932 0.9716605 0.9656613
PM_07_PL_RN_BA_220622 0.9655719 0.9837626 0.9933536 0.9966771
PM_08_PL_RN_BA_220622 0.9749519 0.9909946 0.9929169 0.9914692
pheatmap(rld_cor_PL,
# color=colorRampPalette(c("", "", ""))(200), ## The color palette.
display_numbers=matrix(ifelse(rld_cor_PL>0.99995, "*", ""), nrow(rld_cor_PL)), fontsize_row=5)
# PM
rld_mat_PM=assay(rld_PM)
rld_cor_PM=cor(rld_mat_PM)
head(rld_cor_PM)
PM_01_PM_RN_BA_220613 PM_03_PM_RN_BA_220613 PM_04_PM_RN_BA_220613 PM_06_PM_RN_BA_220613
PM_01_PM_RN_BA_220613 1.0000000 0.9973757 0.9985904 0.9987619
PM_03_PM_RN_BA_220613 0.9973757 1.0000000 0.9965159 0.9967169
PM_04_PM_RN_BA_220613 0.9985904 0.9965159 1.0000000 0.9986604
PM_06_PM_RN_BA_220613 0.9987619 0.9967169 0.9986604 1.0000000
PM_07_PM_RN_BA_220613 0.9982264 0.9965426 0.9982525 0.9985847
PM_08_PM_RN_BA_220613 0.9983503 0.9969611 0.9984054 0.9986196
PM_07_PM_RN_BA_220613 PM_08_PM_RN_BA_220613 PM_09_PM_RN_BA_220613 PM_10_PM_RN_BA_220613
PM_01_PM_RN_BA_220613 0.9982264 0.9983503 0.9988498 0.9978903
PM_03_PM_RN_BA_220613 0.9965426 0.9969611 0.9971380 0.9962465
PM_04_PM_RN_BA_220613 0.9982525 0.9984054 0.9985171 0.9978920
PM_06_PM_RN_BA_220613 0.9985847 0.9986196 0.9983856 0.9977571
PM_07_PM_RN_BA_220613 1.0000000 0.9983868 0.9981028 0.9975851
PM_08_PM_RN_BA_220613 0.9983868 1.0000000 0.9981590 0.9979193
PM_11_PM_RN_BA_220613 PM_12_PM_RN_BA_220613 PM_13_PM_RN_BA_220613 PM_14_PM_RN_BA_220614
PM_01_PM_RN_BA_220613 0.9987021 0.9981612 0.9987513 0.9990086
PM_03_PM_RN_BA_220613 0.9971445 0.9967268 0.9969281 0.9973247
PM_04_PM_RN_BA_220613 0.9987208 0.9984499 0.9985884 0.9988004
PM_06_PM_RN_BA_220613 0.9988347 0.9986200 0.9985122 0.9988176
PM_07_PM_RN_BA_220613 0.9986176 0.9983615 0.9983141 0.9985533
PM_08_PM_RN_BA_220613 0.9986957 0.9987288 0.9981208 0.9985569
PM_15_PM_RN_BA_220614 PM_16_PM_RN_BA_220614 PM_17_PM_RN_BA_220719 PM_18_PM_RN_BA_220614
PM_01_PM_RN_BA_220613 0.9987772 0.9982192 0.9982712 0.9985237
PM_03_PM_RN_BA_220613 0.9968441 0.9964356 0.9970057 0.9970611
PM_04_PM_RN_BA_220613 0.9987165 0.9983014 0.9979029 0.9983203
PM_06_PM_RN_BA_220613 0.9984843 0.9982547 0.9980772 0.9984894
PM_07_PM_RN_BA_220613 0.9981030 0.9980492 0.9975320 0.9984014
PM_08_PM_RN_BA_220613 0.9980320 0.9983454 0.9979423 0.9982200
PM_19_PM_RN_BA_220614 PM_20_PM_RN_BA_220614 PM_21_PM_RN_BA_220614 PM_22_PM_RN_BA_220614
PM_01_PM_RN_BA_220613 0.9983051 0.9923845 0.9984041 0.9974719
PM_03_PM_RN_BA_220613 0.9964955 0.9912391 0.9964212 0.9962789
PM_04_PM_RN_BA_220613 0.9982412 0.9916291 0.9982942 0.9971858
PM_06_PM_RN_BA_220613 0.9977185 0.9922680 0.9978915 0.9978763
PM_07_PM_RN_BA_220613 0.9974616 0.9918770 0.9974821 0.9979052
PM_08_PM_RN_BA_220613 0.9974662 0.9936361 0.9977280 0.9985899
PM_23_PM_RN_BA_220719 PM_24_PM_RN_BA_220614 PM_25_PM_RN_BA_220614 PM_26_PM_RN_BA_220615
PM_01_PM_RN_BA_220613 0.9920253 0.9988268 0.9982624 0.9982536
PM_03_PM_RN_BA_220613 0.9903870 0.9971421 0.9969494 0.9964475
PM_04_PM_RN_BA_220613 0.9921391 0.9989218 0.9980160 0.9980085
PM_06_PM_RN_BA_220613 0.9931568 0.9987724 0.9983870 0.9978310
PM_07_PM_RN_BA_220613 0.9930612 0.9984133 0.9979679 0.9974655
PM_08_PM_RN_BA_220613 0.9945537 0.9985643 0.9988006 0.9976686
PM_27_PM_RN_BA_220615 PM_28_PM_RN_BA_220615 PM_29_PM_RN_BA_220615 PM_30_PM_RN_BA_220615
PM_01_PM_RN_BA_220613 0.9980094 0.9975200 0.9974110 0.9976490
PM_03_PM_RN_BA_220613 0.9960510 0.9961710 0.9966255 0.9964691
PM_04_PM_RN_BA_220613 0.9980859 0.9974715 0.9974742 0.9977576
PM_06_PM_RN_BA_220613 0.9981876 0.9972144 0.9978113 0.9981469
PM_07_PM_RN_BA_220613 0.9979672 0.9966964 0.9977005 0.9981801
PM_08_PM_RN_BA_220613 0.9984926 0.9970235 0.9983968 0.9989239
PM_31_PM_RN_BA_220615 PM_33_PM_RN_BA_220615 PM_36_PM_RN_BA_220615 PM_37_PM_RN_BA_220615
PM_01_PM_RN_BA_220613 0.9977794 0.9973311 0.9986959 0.9985054
PM_03_PM_RN_BA_220613 0.9957756 0.9972537 0.9974029 0.9971031
PM_04_PM_RN_BA_220613 0.9978578 0.9970381 0.9982711 0.9983527
PM_06_PM_RN_BA_220613 0.9981783 0.9972737 0.9982182 0.9984263
PM_07_PM_RN_BA_220613 0.9978319 0.9970234 0.9982121 0.9982932
PM_08_PM_RN_BA_220613 0.9985372 0.9970778 0.9980394 0.9981640
PM_38_PM_RN_BA_220615 PM_39_PM_RN_BA_220615 PM_40_PM_RN_BA_220616 PM_41_PM_RN_BA_220616
PM_01_PM_RN_BA_220613 0.9984232 0.9982751 0.9969698 0.9984334
PM_03_PM_RN_BA_220613 0.9963453 0.9968883 0.9964748 0.9966821
PM_04_PM_RN_BA_220613 0.9984750 0.9979667 0.9969069 0.9985163
PM_06_PM_RN_BA_220613 0.9985518 0.9974829 0.9974060 0.9986491
PM_07_PM_RN_BA_220613 0.9983308 0.9972670 0.9973491 0.9983732
PM_08_PM_RN_BA_220613 0.9985519 0.9971138 0.9979958 0.9983455
PM_42_PM_RN_BA_220616 PM_43_PM_RN_BA_220616 PM_44_PM_RN_BA_220616 PM_45_PM_RN_BA_220616
PM_01_PM_RN_BA_220613 0.9976977 0.9938385 0.9972014 0.9986992
PM_03_PM_RN_BA_220613 0.9960416 0.9923238 0.9958658 0.9970957
PM_04_PM_RN_BA_220613 0.9979002 0.9938706 0.9974646 0.9986205
PM_06_PM_RN_BA_220613 0.9984540 0.9944048 0.9982074 0.9987075
PM_07_PM_RN_BA_220613 0.9982859 0.9946017 0.9977906 0.9983220
PM_08_PM_RN_BA_220613 0.9985283 0.9948631 0.9982159 0.9985586
PM_46_PM_RN_BA_220616 PM_47_PM_RN_BA_220616 PM_48_PM_RN_BA_220719 PM_49_PM_RN_BA_220719
PM_01_PM_RN_BA_220613 0.9987134 0.9962924 0.9984630 0.9976713
PM_03_PM_RN_BA_220613 0.9974159 0.9968620 0.9973488 0.9965995
PM_04_PM_RN_BA_220613 0.9984216 0.9953985 0.9981605 0.9976977
PM_06_PM_RN_BA_220613 0.9981395 0.9959300 0.9981725 0.9980737
PM_07_PM_RN_BA_220613 0.9979053 0.9959119 0.9979720 0.9981773
PM_08_PM_RN_BA_220613 0.9981801 0.9969624 0.9985422 0.9988883
PM_50_PM_RN_BA_220719 PM_51_PM_RN_BA_220920 PM_52_PM_RN_BA_220920 PM_53_PM_RN_BA_220920
PM_01_PM_RN_BA_220613 0.9982543 0.9962985 0.9979297 0.9985335
PM_03_PM_RN_BA_220613 0.9968178 0.9946202 0.9962293 0.9975803
PM_04_PM_RN_BA_220613 0.9979204 0.9964866 0.9981757 0.9983194
PM_06_PM_RN_BA_220613 0.9983608 0.9972383 0.9985343 0.9986150
PM_07_PM_RN_BA_220613 0.9980038 0.9970392 0.9985421 0.9984662
PM_08_PM_RN_BA_220613 0.9985258 0.9980670 0.9986091 0.9985065
PM_54_PM_RN_BA_220920 PM_55_PM_RN_BA_220920 PM_56_PM_RN_BA_220920 PM_57_PM_RN_BA_220920
PM_01_PM_RN_BA_220613 0.9971791 0.9912566 0.9985562 0.9967177
PM_03_PM_RN_BA_220613 0.9951933 0.9907744 0.9962658 0.9960258
PM_04_PM_RN_BA_220613 0.9973352 0.9914641 0.9985717 0.9958937
PM_06_PM_RN_BA_220613 0.9973462 0.9923627 0.9984169 0.9951978
PM_07_PM_RN_BA_220613 0.9970108 0.9926741 0.9979450 0.9949202
PM_08_PM_RN_BA_220613 0.9978954 0.9930917 0.9981147 0.9950165
PM_58_PM_RN_BA_220920
PM_01_PM_RN_BA_220613 0.9977199
PM_03_PM_RN_BA_220613 0.9966960
PM_04_PM_RN_BA_220613 0.9972143
PM_06_PM_RN_BA_220613 0.9967604
PM_07_PM_RN_BA_220613 0.9969378
PM_08_PM_RN_BA_220613 0.9966285
pheatmap(rld_cor_PM,
# color=colorRampPalette(c("", "", ""))(200), ## The color palette.
display_numbers=matrix(ifelse(rld_cor_PM>0.99995, "*", ""), nrow(rld_cor_PM)), fontsize_row=5)
# PS
rld_mat_PS=assay(rld_PS)
rld_cor_PS=cor(rld_mat_PS)
head(rld_cor_PS)
PM_01_PS_RN_BA_220627 PM_02_PS_RN_BA_220627 PM_03_PS_RN_BA_220627 PM_04_PS_RN_BA_220627
PM_01_PS_RN_BA_220627 1.0000000 0.9606450 0.8869079 0.9811236
PM_02_PS_RN_BA_220627 0.9606450 1.0000000 0.8948609 0.9613389
PM_03_PS_RN_BA_220627 0.8869079 0.8948609 1.0000000 0.8837305
PM_04_PS_RN_BA_220627 0.9811236 0.9613389 0.8837305 1.0000000
PM_06_PS_RN_BA_220627 0.9831240 0.9534672 0.8675734 0.9784720
PM_07_PS_RN_BA_220627 0.9837870 0.9539435 0.8616793 0.9879060
PM_06_PS_RN_BA_220627 PM_07_PS_RN_BA_220627 PM_08_PS_RN_BA_220627 PM_11_PS_RN_BA_220627
PM_01_PS_RN_BA_220627 0.9831240 0.9837870 0.7755913 0.9831499
PM_02_PS_RN_BA_220627 0.9534672 0.9539435 0.7936407 0.9512955
PM_03_PS_RN_BA_220627 0.8675734 0.8616793 0.9555562 0.8592704
PM_04_PS_RN_BA_220627 0.9784720 0.9879060 0.7719992 0.9876247
PM_06_PS_RN_BA_220627 1.0000000 0.9857576 0.7482938 0.9859800
PM_07_PS_RN_BA_220627 0.9857576 1.0000000 0.7426523 0.9953816
PM_13_PS_RN_BA_220627 PM_14_PS_RN_BA_220628 PM_15_PS_RN_BA_220714 PM_16_PS_RN_BA_220628
PM_01_PS_RN_BA_220627 0.9759046 0.9849443 0.9826522 0.9865720
PM_02_PS_RN_BA_220627 0.9654009 0.9559280 0.9517803 0.9515115
PM_03_PS_RN_BA_220627 0.9183266 0.8654967 0.8595112 0.8713983
PM_04_PS_RN_BA_220627 0.9758798 0.9864787 0.9877557 0.9833672
PM_06_PS_RN_BA_220627 0.9632580 0.9821050 0.9850207 0.9865273
PM_07_PS_RN_BA_220627 0.9705848 0.9909780 0.9945672 0.9898921
PM_17_PS_RN_BA_220628 PM_18_PS_RN_BA_220628 PM_22_PS_RN_BA_220706 PM_23_PS_RN_BA_220628
PM_01_PS_RN_BA_220627 0.9830409 0.9222055 0.7924112 0.9852402
PM_02_PS_RN_BA_220627 0.9532591 0.9268354 0.8256418 0.9555273
PM_03_PS_RN_BA_220627 0.8583692 0.9658484 0.9265663 0.8810603
PM_04_PS_RN_BA_220627 0.9778826 0.9232192 0.7952704 0.9807152
PM_06_PS_RN_BA_220627 0.9869623 0.9038526 0.7637763 0.9850447
PM_07_PS_RN_BA_220627 0.9873818 0.9045864 0.7658308 0.9867451
PM_24_PS_RN_BA_220628 PM_25_PS_RN_BA_220628 PM_27_PS_RN_BA_220628 PM_29_PS_RN_BA_220628
PM_01_PS_RN_BA_220627 0.7705787 0.9428463 0.8199152 0.8684170
PM_02_PS_RN_BA_220627 0.7869594 0.9483247 0.8388339 0.8788922
PM_03_PS_RN_BA_220627 0.9489078 0.9262650 0.9715437 0.9646375
PM_04_PS_RN_BA_220627 0.7657202 0.9370881 0.8184634 0.8670901
PM_06_PS_RN_BA_220627 0.7423229 0.9241191 0.7945800 0.8472343
PM_07_PS_RN_BA_220627 0.7361024 0.9248893 0.7907046 0.8471786
PM_30_PS_RN_BA_220706 PM_31_PS_RN_BA_220630 PM_33_PS_RN_BA_220630 PM_35_PS_RN_BA_220630
PM_01_PS_RN_BA_220627 0.9799911 0.9266582 0.8227289 0.9811261
PM_02_PS_RN_BA_220627 0.9571301 0.9265018 0.8463876 0.9491906
PM_03_PS_RN_BA_220627 0.8880815 0.9404344 0.9393521 0.8687451
PM_04_PS_RN_BA_220627 0.9669150 0.9287156 0.8204771 0.9775745
PM_06_PS_RN_BA_220627 0.9665778 0.9145174 0.7961030 0.9900084
PM_07_PS_RN_BA_220627 0.9696289 0.9142789 0.7957224 0.9865509
PM_36_PS_RN_BA_220630 PM_37_PS_RN_BA_220630 PM_38_PS_RN_BA_220630 PM_39_PS_RN_BA_220630
PM_01_PS_RN_BA_220627 0.9822193 0.9613822 0.9772613 0.9560031
PM_02_PS_RN_BA_220627 0.9535603 0.9580780 0.9583482 0.9669563
PM_03_PS_RN_BA_220627 0.8771271 0.8948083 0.8950081 0.8897202
PM_04_PS_RN_BA_220627 0.9833107 0.9544933 0.9834811 0.9504566
PM_06_PS_RN_BA_220627 0.9798069 0.9492800 0.9787246 0.9487470
PM_07_PS_RN_BA_220627 0.9858837 0.9523343 0.9809602 0.9450771
PM_40_PS_RN_BA_220630 PM_42_PS_RN_BA_220630 PM_43_PS_RN_BA_220630 PM_44_PS_RN_BA_220630
PM_01_PS_RN_BA_220627 0.7608749 0.8575425 0.7754604 0.9790492
PM_02_PS_RN_BA_220627 0.8005038 0.8866383 0.7990796 0.9587052
PM_03_PS_RN_BA_220627 0.8929633 0.9311835 0.9523295 0.8659574
PM_04_PS_RN_BA_220627 0.7698169 0.8616466 0.7750679 0.9804302
PM_06_PS_RN_BA_220627 0.7335992 0.8342816 0.7472615 0.9895171
PM_07_PS_RN_BA_220627 0.7395331 0.8390196 0.7443485 0.9859969
PM_47_PS_RN_BA_220706 PM_48_PS_RN_BA_220714 PM_49_PS_RN_BA_220714 PM_50_PS_RN_BA_220714
PM_01_PS_RN_BA_220627 0.8937449 0.9550159 0.7635917 0.9821767
PM_02_PS_RN_BA_220627 0.9010620 0.9608715 0.7783919 0.9596088
PM_03_PS_RN_BA_220627 0.9886646 0.9292460 0.9524269 0.8791126
PM_04_PS_RN_BA_220627 0.8920354 0.9551523 0.7581308 0.9899734
PM_06_PS_RN_BA_220627 0.8753964 0.9453726 0.7362177 0.9830774
PM_07_PS_RN_BA_220627 0.8716289 0.9441045 0.7287314 0.9915495
PM_52_PS_RN_BA_220921 PM_53_PS_RN_BA_220921 PM_54_PS_RN_BA_220921 PM_55_PS_RN_BA_220921
PM_01_PS_RN_BA_220627 0.9643769 0.7397205 0.9767794 0.7744572
PM_02_PS_RN_BA_220627 0.9581497 0.7510873 0.9769598 0.7982340
PM_03_PS_RN_BA_220627 0.9320928 0.9361341 0.8959031 0.9461651
PM_04_PS_RN_BA_220627 0.9641459 0.7316286 0.9765499 0.7739316
PM_06_PS_RN_BA_220627 0.9520287 0.7127788 0.9753158 0.7466492
PM_07_PS_RN_BA_220627 0.9525584 0.7025515 0.9728629 0.7436672
PM_56_PS_RN_BA_220921 PM_57_PS_RN_BA_220921 PM_58_PS_RN_BA_220921
PM_01_PS_RN_BA_220627 0.9843753 0.9088365 0.9800834
PM_02_PS_RN_BA_220627 0.9556750 0.9198265 0.9604107
PM_03_PS_RN_BA_220627 0.8687782 0.9714176 0.8868092
PM_04_PS_RN_BA_220627 0.9837168 0.9056669 0.9882764
PM_06_PS_RN_BA_220627 0.9845354 0.8874210 0.9786215
PM_07_PS_RN_BA_220627 0.9907088 0.8851451 0.9888248
pheatmap(rld_cor_PS,
# color=colorRampPalette(c("", "", ""))(200), ## The color palette.
display_numbers=matrix(ifelse(rld_cor_PS>0.99995, "*", ""), nrow(rld_cor_PS)), fontsize_row=5)
Performs the DESeq2 analysis.
dds_PB=DESeq(dds_PB)
using pre-existing size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 872 genes
-- DESeq argument 'minReplicatesForReplace' = 7
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
dds_PL=DESeq(dds_PL)
using pre-existing size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 553 genes
-- DESeq argument 'minReplicatesForReplace' = 7
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
dds_PM=DESeq(dds_PM)
using pre-existing size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 1023 genes
-- DESeq argument 'minReplicatesForReplace' = 7
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
dds_PS=DESeq(dds_PS)
using pre-existing size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 5513 genes
-- DESeq argument 'minReplicatesForReplace' = 7
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
Step-by-step breakdown to the DESeq2 analysis: Step 1: estimate size factors (check size factors, total number of raw counts per sample, total number of normalized counts per sample)
# PB
sizeFactors(dds_PB)
PM_01_PB_RN_BA_220603 PM_03_PB_RN_BA_220607 PM_04_PB_RN_BA_220607 PM_05_PB_RN_BA_220606 PM_06_PB_RN_BA_220606
1.9927487 0.6533585 1.5179348 0.0765266 1.2907222
PM_07_PB_RN_BA_220608 PM_08_PB_RN_BA_220608 PM_09_PB_RN_BA_220608 PM_10_PB_RN_BA_220608 PM_11_PB_RN_BA_220608
1.3418045 0.6808323 1.4021222 0.7180940 1.7293236
PM_12_PB_RN_BA_220608 PM_13_PB_RN_BA_220608 PM_14_PB_RN_BA_220606 PM_15_PB_RN_BA_220606 PM_16_PB_RN_BA_220607
0.3475635 1.4749431 1.4372030 1.2532442 1.4619712
PM_17_PB_RN_BA_220718 PM_18_PB_RN_BA_220608 PM_19_PB_RN_BA_220607 PM_20_PB_RN_BA_220607 PM_22_PB_RN_BA_220603
0.6044018 1.9538787 1.9826549 1.0651692 1.1004654
PM_23_PB_RN_BA_220606 PM_24_PB_RN_BA_220606 PM_25_PB_RN_BA_220607 PM_26_PB_RN_BA_220608 PM_27_PB_RN_BA_220603
1.0663400 1.3711046 1.0738322 1.3351838 1.0590360
PM_28_PB_RN_BA_220607 PM_29_PB_RN_BA_220608 PM_30_PB_RN_BA_220603 PM_31_PB_RN_BA_220603 PM_32_PB_RN_BA_220603
1.4048117 0.2811831 0.2582569 0.7841821 1.0484875
PM_35_PB_RN_BA_220606 PM_36_PB_RN_BA_220603 PM_37_PB_RN_BA_220606 PM_38_PB_RN_BA_220609 PM_39_PB_RN_BA_220609
0.1375362 1.5660193 1.6073976 1.5682022 1.8185940
PM_40_PB_RN_BA_220603 PM_41_PB_RN_BA_220609 PM_42_PB_RN_BA_220609 PM_43_PB_RN_BA_220609 PM_44_PB_RN_BA_220609
1.4885983 0.3197314 1.1066387 0.8736601 1.4449211
PM_45_PB_RN_BA_220609 PM_46_PB_RN_BA_220609 PM_47_PB_RN_BA_220609 PM_48_PB_RN_BA_220718 PM_49_PB_RN_BA_220718
1.4280626 1.7844170 1.0444030 1.7825095 0.6698575
PM_50_PB_RN_BA_220718 PM_52_PB_RN_BA_220916 PM_53_PB_RN_BA_220916 PM_54_PB_RN_BA_220916 PM_55_PB_RN_BA_220916
1.5285517 1.3123926 0.5113159 0.8603980 0.7516655
PM_56_PB_RN_BA_220916 PM_57_PB_RN_BA_220916 PM_58_PB_RN_BA_220916
1.5282431 1.7239848 1.6094768
colSums(counts(dds_PB))
PM_01_PB_RN_BA_220603 PM_03_PB_RN_BA_220607 PM_04_PB_RN_BA_220607 PM_05_PB_RN_BA_220606 PM_06_PB_RN_BA_220606
32610758 26068281 24419181 62243285 28494916
PM_07_PB_RN_BA_220608 PM_08_PB_RN_BA_220608 PM_09_PB_RN_BA_220608 PM_10_PB_RN_BA_220608 PM_11_PB_RN_BA_220608
23656591 24315300 22753815 15612562 26675103
PM_12_PB_RN_BA_220608 PM_13_PB_RN_BA_220608 PM_14_PB_RN_BA_220606 PM_15_PB_RN_BA_220606 PM_16_PB_RN_BA_220607
56362531 23912393 26536451 22286059 24583380
PM_17_PB_RN_BA_220718 PM_18_PB_RN_BA_220608 PM_19_PB_RN_BA_220607 PM_20_PB_RN_BA_220607 PM_22_PB_RN_BA_220603
18031630 34029949 34548855 27721241 23266383
PM_23_PB_RN_BA_220606 PM_24_PB_RN_BA_220606 PM_25_PB_RN_BA_220607 PM_26_PB_RN_BA_220608 PM_27_PB_RN_BA_220603
19886101 21047296 18830666 24774561 18055672
PM_28_PB_RN_BA_220607 PM_29_PB_RN_BA_220608 PM_30_PB_RN_BA_220603 PM_31_PB_RN_BA_220603 PM_32_PB_RN_BA_220603
26475749 9276923 11816620 15226024 21441796
PM_35_PB_RN_BA_220606 PM_36_PB_RN_BA_220603 PM_37_PB_RN_BA_220606 PM_38_PB_RN_BA_220609 PM_39_PB_RN_BA_220609
10168756 29010123 30849112 30061293 28351551
PM_40_PB_RN_BA_220603 PM_41_PB_RN_BA_220609 PM_42_PB_RN_BA_220609 PM_43_PB_RN_BA_220609 PM_44_PB_RN_BA_220609
25347844 23677004 18959347 22222748 29381822
PM_45_PB_RN_BA_220609 PM_46_PB_RN_BA_220609 PM_47_PB_RN_BA_220609 PM_48_PB_RN_BA_220718 PM_49_PB_RN_BA_220718
21823454 33633562 19382950 30596372 19165835
PM_50_PB_RN_BA_220718 PM_52_PB_RN_BA_220916 PM_53_PB_RN_BA_220916 PM_54_PB_RN_BA_220916 PM_55_PB_RN_BA_220916
27842815 20237292 13092332 13189986 19135364
PM_56_PB_RN_BA_220916 PM_57_PB_RN_BA_220916 PM_58_PB_RN_BA_220916
24219328 30964760 26185826
colSums(counts(dds_PB, normalized=TRUE))
PM_01_PB_RN_BA_220603 PM_03_PB_RN_BA_220607 PM_04_PB_RN_BA_220607 PM_05_PB_RN_BA_220606 PM_06_PB_RN_BA_220606
16364712 39898892 16087108 813354863 22076722
PM_07_PB_RN_BA_220608 PM_08_PB_RN_BA_220608 PM_09_PB_RN_BA_220608 PM_10_PB_RN_BA_220608 PM_11_PB_RN_BA_220608
17630431 35714080 16228125 21741670 15425166
PM_12_PB_RN_BA_220608 PM_13_PB_RN_BA_220608 PM_14_PB_RN_BA_220606 PM_15_PB_RN_BA_220606 PM_16_PB_RN_BA_220607
162164702 16212418 18463954 17782694 16815229
PM_17_PB_RN_BA_220718 PM_18_PB_RN_BA_220608 PM_19_PB_RN_BA_220607 PM_20_PB_RN_BA_220607 PM_22_PB_RN_BA_220603
29833845 17416613 17425551 26025201 21142312
PM_23_PB_RN_BA_220606 PM_24_PB_RN_BA_220606 PM_25_PB_RN_BA_220607 PM_26_PB_RN_BA_220608 PM_27_PB_RN_BA_220603
18648931 15350613 17535948 18555169 17049158
PM_28_PB_RN_BA_220607 PM_29_PB_RN_BA_220608 PM_30_PB_RN_BA_220603 PM_31_PB_RN_BA_220603 PM_32_PB_RN_BA_220603
18846475 32992459 45755290 19416439 20450216
PM_35_PB_RN_BA_220606 PM_36_PB_RN_BA_220603 PM_37_PB_RN_BA_220606 PM_38_PB_RN_BA_220609 PM_39_PB_RN_BA_220609
73935128 18524754 19191960 19169271 15589819
PM_40_PB_RN_BA_220603 PM_41_PB_RN_BA_220609 PM_42_PB_RN_BA_220609 PM_43_PB_RN_BA_220609 PM_44_PB_RN_BA_220609
17027995 74052800 17132373 25436377 20334550
PM_45_PB_RN_BA_220609 PM_46_PB_RN_BA_220609 PM_47_PB_RN_BA_220609 PM_48_PB_RN_BA_220718 PM_49_PB_RN_BA_220718
15281861 18848488 18558879 17164774 28611808
PM_50_PB_RN_BA_220718 PM_52_PB_RN_BA_220916 PM_53_PB_RN_BA_220916 PM_54_PB_RN_BA_220916 PM_55_PB_RN_BA_220916
18215161 15420151 25605172 15330099 25457286
PM_56_PB_RN_BA_220916 PM_57_PB_RN_BA_220916 PM_58_PB_RN_BA_220916
15847824 17961156 16269776
# PL
sizeFactors(dds_PL)
PM_01_PL_RN_BA_220622 PM_03_PL_RN_BA_220622 PM_04_PL_RN_BA_220622 PM_06_PL_RN_BA_220622 PM_07_PL_RN_BA_220622
1.40297862 0.46367867 2.85696861 0.05537568 2.80900092
PM_08_PL_RN_BA_220622 PM_09_PL_RN_BA_220622 PM_10_PL_RN_BA_220622 PM_11_PL_RN_BA_220622 PM_12_PL_RN_BA_220622
1.58957104 2.00337729 1.17858581 2.09246939 0.36782415
PM_13_PL_RN_BA_220622 PM_14_PL_RN_BA_220720 PM_15_PL_RN_BA_220622 PM_16_PL_RN_BA_220622 PM_17_PL_RN_BA_220622
2.66309782 3.69152590 1.46677811 2.49919006 0.35155524
PM_18_PL_RN_BA_220622 PM_19_PL_RN_BA_220622 PM_20_PL_RN_BA_220622 PM_21_PL_RN_BA_220622 PM_22_PL_RN_BA_220720
0.05415938 2.05344480 0.03270570 2.45919055 0.02814269
PM_24_PL_RN_BA_220622 PM_25_PL_RN_BA_220622 PM_26_PL_RN_BA_220623 PM_27_PL_RN_BA_220623 PM_28_PL_RN_BA_220623
0.09620581 0.86529255 2.24341313 2.02138616 1.50732385
PM_29_PL_RN_BA_220623 PM_30_PL_RN_BA_220720 PM_31_PL_RN_BA_220623 PM_32_PL_RN_BA_220720 PM_33_PL_RN_BA_220623
0.41538316 1.85433641 0.68014153 1.32914721 0.56766953
PM_36_PL_RN_BA_220623 PM_37_PL_RN_BA_220623 PM_38_PL_RN_BA_220623 PM_39_PL_RN_BA_220720 PM_40_PL_RN_BA_220624
1.33532765 2.62394196 3.02302763 0.35657666 3.11030096
PM_41_PL_RN_BA_220624 PM_42_PL_RN_BA_220624 PM_43_PL_RN_BA_220720 PM_44_PL_RN_BA_220624 PM_45_PL_RN_BA_220624
1.05301075 1.51788168 0.30244477 0.85706746 2.18470276
PM_46_PL_RN_BA_220624 PM_47_PL_RN_BA_220624 PM_48_PL_RN_BA_220720 PM_49_PL_RN_BA_220720 PM_50_PL_RN_BA_220720
1.13637430 1.40400835 1.60445764 2.78561615 2.15553491
PM_51_PL_RN_BA_220920 PM_52_PL_RN_BA_220920 PM_54_PL_RN_BA_220920 PM_55_PL_RN_BA_220920 PM_56_PL_RN_BA_220923
2.82774982 0.89931011 1.50238908 0.28197242 1.50707690
PM_57_PL_RN_BA_220920 PM_58_PL_RN_BA_220923
2.06168949 3.80973657
colSums(counts(dds_PL))
PM_01_PL_RN_BA_220622 PM_03_PL_RN_BA_220622 PM_04_PL_RN_BA_220622 PM_06_PL_RN_BA_220622 PM_07_PL_RN_BA_220622
16972118 6310158 28707632 2338314 32026511
PM_08_PL_RN_BA_220622 PM_09_PL_RN_BA_220622 PM_10_PL_RN_BA_220622 PM_11_PL_RN_BA_220622 PM_12_PL_RN_BA_220622
18320239 23015351 15700609 23659162 5421460
PM_13_PL_RN_BA_220622 PM_14_PL_RN_BA_220720 PM_15_PL_RN_BA_220622 PM_16_PL_RN_BA_220622 PM_17_PL_RN_BA_220622
28759574 47565046 18565452 27103526 7078831
PM_18_PL_RN_BA_220622 PM_19_PL_RN_BA_220622 PM_20_PL_RN_BA_220622 PM_21_PL_RN_BA_220622 PM_22_PL_RN_BA_220720
1241466 22100675 1027920 28509685 2501194
PM_24_PL_RN_BA_220622 PM_25_PL_RN_BA_220622 PM_26_PL_RN_BA_220623 PM_27_PL_RN_BA_220623 PM_28_PL_RN_BA_220623
2117416 12270997 21533902 21784421 20048331
PM_29_PL_RN_BA_220623 PM_30_PL_RN_BA_220720 PM_31_PL_RN_BA_220623 PM_32_PL_RN_BA_220720 PM_33_PL_RN_BA_220623
4401669 24962592 8000917 13171711 7748093
PM_36_PL_RN_BA_220623 PM_37_PL_RN_BA_220623 PM_38_PL_RN_BA_220623 PM_39_PL_RN_BA_220720 PM_40_PL_RN_BA_220624
17344782 28610826 28423069 7481797 30004923
PM_41_PL_RN_BA_220624 PM_42_PL_RN_BA_220624 PM_43_PL_RN_BA_220720 PM_44_PL_RN_BA_220624 PM_45_PL_RN_BA_220624
13261902 17371986 4187739 9030104 22140191
PM_46_PL_RN_BA_220624 PM_47_PL_RN_BA_220624 PM_48_PL_RN_BA_220720 PM_49_PL_RN_BA_220720 PM_50_PL_RN_BA_220720
13415213 14740941 20500476 27367289 22692850
PM_51_PL_RN_BA_220920 PM_52_PL_RN_BA_220920 PM_54_PL_RN_BA_220920 PM_55_PL_RN_BA_220920 PM_56_PL_RN_BA_220923
28682485 14346107 16788965 3608214 19865015
PM_57_PL_RN_BA_220920 PM_58_PL_RN_BA_220923
23736194 43400131
colSums(counts(dds_PL, normalized=TRUE))
PM_01_PL_RN_BA_220622 PM_03_PL_RN_BA_220622 PM_04_PL_RN_BA_220622 PM_06_PL_RN_BA_220622 PM_07_PL_RN_BA_220622
12097204 13608903 10048284 42226371 11401389
PM_08_PL_RN_BA_220622 PM_09_PL_RN_BA_220622 PM_10_PL_RN_BA_220622 PM_11_PL_RN_BA_220622 PM_12_PL_RN_BA_220622
11525272 11488276 13321566 11306814 14739271
PM_13_PL_RN_BA_220622 PM_14_PL_RN_BA_220720 PM_15_PL_RN_BA_220622 PM_16_PL_RN_BA_220622 PM_17_PL_RN_BA_220622
10799293 12884928 12657301 10844924 20135757
PM_18_PL_RN_BA_220622 PM_19_PL_RN_BA_220622 PM_20_PL_RN_BA_220622 PM_21_PL_RN_BA_220622 PM_22_PL_RN_BA_220720
22922456 10762732 31429389 11593117 88875427
PM_24_PL_RN_BA_220622 PM_25_PL_RN_BA_220622 PM_26_PL_RN_BA_220623 PM_27_PL_RN_BA_220623 PM_28_PL_RN_BA_220623
22009232 14181327 9598723 10776971 13300613
PM_29_PL_RN_BA_220623 PM_30_PL_RN_BA_220720 PM_31_PL_RN_BA_220623 PM_32_PL_RN_BA_220720 PM_33_PL_RN_BA_220623
10596648 13461739 11763606 9909896 13648950
PM_36_PL_RN_BA_220623 PM_37_PL_RN_BA_220623 PM_38_PL_RN_BA_220623 PM_39_PL_RN_BA_220720 PM_40_PL_RN_BA_220624
12989158 10903757 9402186 20982296 9646952
PM_41_PL_RN_BA_220624 PM_42_PL_RN_BA_220624 PM_43_PL_RN_BA_220720 PM_44_PL_RN_BA_220624 PM_45_PL_RN_BA_220624
12594270 11444888 13846293 10536048 10134189
PM_46_PL_RN_BA_220624 PM_47_PL_RN_BA_220624 PM_48_PL_RN_BA_220720 PM_49_PL_RN_BA_220720 PM_50_PL_RN_BA_220720
11805277 10499183 12777200 9824501 10527712
PM_51_PL_RN_BA_220920 PM_52_PL_RN_BA_220920 PM_54_PL_RN_BA_220920 PM_55_PL_RN_BA_220920 PM_56_PL_RN_BA_220923
10143219 15952347 11174845 12796337 13181155
PM_57_PL_RN_BA_220920 PM_58_PL_RN_BA_220923
11512982 11391898
# PM
sizeFactors(dds_PM)
PM_01_PM_RN_BA_220613 PM_03_PM_RN_BA_220613 PM_04_PM_RN_BA_220613 PM_06_PM_RN_BA_220613 PM_07_PM_RN_BA_220613
1.5088832 1.4255698 1.6204817 1.1970366 0.9800432
PM_08_PM_RN_BA_220613 PM_09_PM_RN_BA_220613 PM_10_PM_RN_BA_220613 PM_11_PM_RN_BA_220613 PM_12_PM_RN_BA_220613
1.2978390 1.4868746 1.1757672 1.1081042 1.1678733
PM_13_PM_RN_BA_220613 PM_14_PM_RN_BA_220614 PM_15_PM_RN_BA_220614 PM_16_PM_RN_BA_220614 PM_17_PM_RN_BA_220719
1.1631088 1.8022065 1.0390151 1.2341993 0.6961912
PM_18_PM_RN_BA_220614 PM_19_PM_RN_BA_220614 PM_20_PM_RN_BA_220614 PM_21_PM_RN_BA_220614 PM_22_PM_RN_BA_220614
0.9084388 1.4110597 0.2821348 1.3634059 0.7242727
PM_23_PM_RN_BA_220719 PM_24_PM_RN_BA_220614 PM_25_PM_RN_BA_220614 PM_26_PM_RN_BA_220615 PM_27_PM_RN_BA_220615
0.4507395 1.3024728 1.3538257 1.4185118 1.2038617
PM_28_PM_RN_BA_220615 PM_29_PM_RN_BA_220615 PM_30_PM_RN_BA_220615 PM_31_PM_RN_BA_220615 PM_33_PM_RN_BA_220615
1.3040438 0.6821155 1.1197386 1.5807151 0.1739217
PM_36_PM_RN_BA_220615 PM_37_PM_RN_BA_220615 PM_38_PM_RN_BA_220615 PM_39_PM_RN_BA_220615 PM_40_PM_RN_BA_220616
1.6355448 0.4399030 1.1170166 1.2976556 0.1718784
PM_41_PM_RN_BA_220616 PM_42_PM_RN_BA_220616 PM_43_PM_RN_BA_220616 PM_44_PM_RN_BA_220616 PM_45_PM_RN_BA_220616
0.8429614 0.7767747 0.9374023 1.2488331 1.3410268
PM_46_PM_RN_BA_220616 PM_47_PM_RN_BA_220616 PM_48_PM_RN_BA_220719 PM_49_PM_RN_BA_220719 PM_50_PM_RN_BA_220719
1.1273183 1.1636958 1.4509711 1.2990044 0.5381660
PM_51_PM_RN_BA_220920 PM_52_PM_RN_BA_220920 PM_53_PM_RN_BA_220920 PM_54_PM_RN_BA_220920 PM_55_PM_RN_BA_220920
1.4074151 1.0996795 2.1120325 1.5010878 0.1360339
PM_56_PM_RN_BA_220920 PM_57_PM_RN_BA_220920 PM_58_PM_RN_BA_220920
1.6941223 1.5719031 1.5294640
colSums(counts(dds_PM))
PM_01_PM_RN_BA_220613 PM_03_PM_RN_BA_220613 PM_04_PM_RN_BA_220613 PM_06_PM_RN_BA_220613 PM_07_PM_RN_BA_220613
46376872 55014579 43656352 35010196 31094263
PM_08_PM_RN_BA_220613 PM_09_PM_RN_BA_220613 PM_10_PM_RN_BA_220613 PM_11_PM_RN_BA_220613 PM_12_PM_RN_BA_220613
35346930 47533431 39954181 33786727 30577498
PM_13_PM_RN_BA_220613 PM_14_PM_RN_BA_220614 PM_15_PM_RN_BA_220614 PM_16_PM_RN_BA_220614 PM_17_PM_RN_BA_220719
39998832 56743258 32843218 34735488 23196778
PM_18_PM_RN_BA_220614 PM_19_PM_RN_BA_220614 PM_20_PM_RN_BA_220614 PM_21_PM_RN_BA_220614 PM_22_PM_RN_BA_220614
30517047 47242029 7128026 37241719 19605904
PM_23_PM_RN_BA_220719 PM_24_PM_RN_BA_220614 PM_25_PM_RN_BA_220614 PM_26_PM_RN_BA_220615 PM_27_PM_RN_BA_220615
8539845 41821260 38103175 39848440 31824294
PM_28_PM_RN_BA_220615 PM_29_PM_RN_BA_220615 PM_30_PM_RN_BA_220615 PM_31_PM_RN_BA_220615 PM_33_PM_RN_BA_220615
39570810 18318893 27867163 37619368 7699167
PM_36_PM_RN_BA_220615 PM_37_PM_RN_BA_220615 PM_38_PM_RN_BA_220615 PM_39_PM_RN_BA_220615 PM_40_PM_RN_BA_220616
52832937 14847347 30443900 50012859 5653052
PM_41_PM_RN_BA_220616 PM_42_PM_RN_BA_220616 PM_43_PM_RN_BA_220616 PM_44_PM_RN_BA_220616 PM_45_PM_RN_BA_220616
25155474 21967858 29089204 33313500 37297150
PM_46_PM_RN_BA_220616 PM_47_PM_RN_BA_220616 PM_48_PM_RN_BA_220719 PM_49_PM_RN_BA_220719 PM_50_PM_RN_BA_220719
38935423 28303024 44317920 35130201 16121758
PM_51_PM_RN_BA_220920 PM_52_PM_RN_BA_220920 PM_53_PM_RN_BA_220920 PM_54_PM_RN_BA_220920 PM_55_PM_RN_BA_220920
32136838 34048376 64047561 35428780 4583139
PM_56_PM_RN_BA_220920 PM_57_PM_RN_BA_220920 PM_58_PM_RN_BA_220920
43632109 63412833 62072224
colSums(counts(dds_PM, normalized=TRUE))
PM_01_PM_RN_BA_220613 PM_03_PM_RN_BA_220613 PM_04_PM_RN_BA_220613 PM_06_PM_RN_BA_220613 PM_07_PM_RN_BA_220613
30735892 38591291 26940355 29247390 31727440
PM_08_PM_RN_BA_220613 PM_09_PM_RN_BA_220613 PM_10_PM_RN_BA_220613 PM_11_PM_RN_BA_220613 PM_12_PM_RN_BA_220613
27235219 31968689 33981372 30490569 26182204
PM_13_PM_RN_BA_220613 PM_14_PM_RN_BA_220614 PM_15_PM_RN_BA_220614 PM_16_PM_RN_BA_220614 PM_17_PM_RN_BA_220719
34389588 31485437 31609952 28144149 33319550
PM_18_PM_RN_BA_220614 PM_19_PM_RN_BA_220614 PM_20_PM_RN_BA_220614 PM_21_PM_RN_BA_220614 PM_22_PM_RN_BA_220614
33592850 33479822 25264610 27315211 27069782
PM_23_PM_RN_BA_220719 PM_24_PM_RN_BA_220614 PM_25_PM_RN_BA_220614 PM_26_PM_RN_BA_220615 PM_27_PM_RN_BA_220615
18946298 32109124 28144817 28091723 26435175
PM_28_PM_RN_BA_220615 PM_29_PM_RN_BA_220615 PM_30_PM_RN_BA_220615 PM_31_PM_RN_BA_220615 PM_33_PM_RN_BA_220615
30344694 26855997 24887203 23798956 44268004
PM_36_PM_RN_BA_220615 PM_37_PM_RN_BA_220615 PM_38_PM_RN_BA_220615 PM_39_PM_RN_BA_220615 PM_40_PM_RN_BA_220616
32302959 33751408 27254654 38540934 32889841
PM_41_PM_RN_BA_220616 PM_42_PM_RN_BA_220616 PM_43_PM_RN_BA_220616 PM_44_PM_RN_BA_220616 PM_45_PM_RN_BA_220616
29841788 28280860 31031718 26675702 27812383
PM_46_PM_RN_BA_220616 PM_47_PM_RN_BA_220616 PM_48_PM_RN_BA_220719 PM_49_PM_RN_BA_220719 PM_50_PM_RN_BA_220719
34538092 24321668 30543627 27043944 29956850
PM_51_PM_RN_BA_220920 PM_52_PM_RN_BA_220920 PM_53_PM_RN_BA_220920 PM_54_PM_RN_BA_220920 PM_55_PM_RN_BA_220920
22833944 30962092 30325084 23602071 33691145
PM_56_PM_RN_BA_220920 PM_57_PM_RN_BA_220920 PM_58_PM_RN_BA_220920
25754994 40341438 40584299
# PS
sizeFactors(dds_PS)
PM_01_PS_RN_BA_220627 PM_02_PS_RN_BA_220627 PM_03_PS_RN_BA_220627 PM_04_PS_RN_BA_220627 PM_06_PS_RN_BA_220627
1.058600541 0.658393088 1.222699830 1.834500045 4.842249194
PM_07_PS_RN_BA_220627 PM_08_PS_RN_BA_220627 PM_11_PS_RN_BA_220627 PM_13_PS_RN_BA_220627 PM_14_PS_RN_BA_220628
8.126935887 0.398138600 7.237634747 4.384803112 6.038798058
PM_15_PS_RN_BA_220714 PM_16_PS_RN_BA_220628 PM_17_PS_RN_BA_220628 PM_18_PS_RN_BA_220628 PM_22_PS_RN_BA_220706
8.271116876 6.656457794 8.197012251 0.642441448 0.067976354
PM_23_PS_RN_BA_220628 PM_24_PS_RN_BA_220628 PM_25_PS_RN_BA_220628 PM_27_PS_RN_BA_220628 PM_29_PS_RN_BA_220628
12.952010726 0.208645715 0.090946482 0.643588097 0.004792714
PM_30_PS_RN_BA_220706 PM_31_PS_RN_BA_220630 PM_33_PS_RN_BA_220630 PM_35_PS_RN_BA_220630 PM_36_PS_RN_BA_220630
1.432265234 0.110024402 0.043707832 7.299222374 3.312055627
PM_37_PS_RN_BA_220630 PM_38_PS_RN_BA_220630 PM_39_PS_RN_BA_220630 PM_40_PS_RN_BA_220630 PM_42_PS_RN_BA_220630
1.806424448 1.457886987 1.219826843 0.076546734 0.052348316
PM_43_PS_RN_BA_220630 PM_44_PS_RN_BA_220630 PM_47_PS_RN_BA_220706 PM_48_PS_RN_BA_220714 PM_49_PS_RN_BA_220714
0.272828739 6.442905955 1.594276066 0.361234164 0.550746554
PM_50_PS_RN_BA_220714 PM_52_PS_RN_BA_220921 PM_53_PS_RN_BA_220921 PM_54_PS_RN_BA_220921 PM_55_PS_RN_BA_220921
5.975779874 2.006948025 0.479613646 3.565880650 0.126282891
PM_56_PS_RN_BA_220921 PM_57_PS_RN_BA_220921 PM_58_PS_RN_BA_220921
4.111315224 0.544084792 5.130149915
colSums(counts(dds_PS))
PM_01_PS_RN_BA_220627 PM_02_PS_RN_BA_220627 PM_03_PS_RN_BA_220627 PM_04_PS_RN_BA_220627 PM_06_PS_RN_BA_220627
3526585 11237092 3717571 9209448 20212599
PM_07_PS_RN_BA_220627 PM_08_PS_RN_BA_220627 PM_11_PS_RN_BA_220627 PM_13_PS_RN_BA_220627 PM_14_PS_RN_BA_220628
22166010 2972093 20998687 10983762 17703875
PM_15_PS_RN_BA_220714 PM_16_PS_RN_BA_220628 PM_17_PS_RN_BA_220628 PM_18_PS_RN_BA_220628 PM_22_PS_RN_BA_220706
23441702 19927882 21618079 2387938 377162
PM_23_PS_RN_BA_220628 PM_24_PS_RN_BA_220628 PM_25_PS_RN_BA_220628 PM_27_PS_RN_BA_220628 PM_29_PS_RN_BA_220628
30967209 3413665 765909 2435442 3363
PM_30_PS_RN_BA_220706 PM_31_PS_RN_BA_220630 PM_33_PS_RN_BA_220630 PM_35_PS_RN_BA_220630 PM_36_PS_RN_BA_220630
5795113 393254 421510 21057303 11732248
PM_37_PS_RN_BA_220630 PM_38_PS_RN_BA_220630 PM_39_PS_RN_BA_220630 PM_40_PS_RN_BA_220630 PM_42_PS_RN_BA_220630
11255323 4391577 15039306 446009 212373
PM_43_PS_RN_BA_220630 PM_44_PS_RN_BA_220630 PM_47_PS_RN_BA_220706 PM_48_PS_RN_BA_220714 PM_49_PS_RN_BA_220714
1284454 21334116 4501689 3702169 2662658
PM_50_PS_RN_BA_220714 PM_52_PS_RN_BA_220921 PM_53_PS_RN_BA_220921 PM_54_PS_RN_BA_220921 PM_55_PS_RN_BA_220921
16254131 8875700 2631656 18026037 712422
PM_56_PS_RN_BA_220921 PM_57_PS_RN_BA_220921 PM_58_PS_RN_BA_220921
13346141 1650953 16560610
colSums(counts(dds_PS, normalized=TRUE))
PM_01_PS_RN_BA_220627 PM_02_PS_RN_BA_220627 PM_03_PS_RN_BA_220627 PM_04_PS_RN_BA_220627 PM_06_PS_RN_BA_220627
3331365.2 17067451.4 3040460.9 5020140.5 4174217.0
PM_07_PS_RN_BA_220627 PM_08_PS_RN_BA_220627 PM_11_PS_RN_BA_220627 PM_13_PS_RN_BA_220627 PM_14_PS_RN_BA_220628
2727474.5 7464970.7 2901319.0 2504961.3 2931688.5
PM_15_PS_RN_BA_220714 PM_16_PS_RN_BA_220628 PM_17_PS_RN_BA_220628 PM_18_PS_RN_BA_220628 PM_22_PS_RN_BA_220706
2834164.0 2993766.7 2637312.0 3716973.8 5548429.4
PM_23_PS_RN_BA_220628 PM_24_PS_RN_BA_220628 PM_25_PS_RN_BA_220628 PM_27_PS_RN_BA_220628 PM_29_PS_RN_BA_220628
2390919.0 16361059.7 8421535.2 3784162.6 701690.2
PM_30_PS_RN_BA_220706 PM_31_PS_RN_BA_220630 PM_33_PS_RN_BA_220630 PM_35_PS_RN_BA_220630 PM_36_PS_RN_BA_220630
4046117.2 3574243.5 9643809.5 2884869.4 3542285.9
PM_37_PS_RN_BA_220630 PM_38_PS_RN_BA_220630 PM_39_PS_RN_BA_220630 PM_40_PS_RN_BA_220630 PM_42_PS_RN_BA_220630
6230718.9 3012289.0 12329049.9 5826623.5 4056921.4
PM_43_PS_RN_BA_220630 PM_44_PS_RN_BA_220630 PM_47_PS_RN_BA_220706 PM_48_PS_RN_BA_220714 PM_49_PS_RN_BA_220714
4707913.1 3311256.8 2823657.1 10248668.0 4834634.0
PM_50_PS_RN_BA_220714 PM_52_PS_RN_BA_220921 PM_53_PS_RN_BA_220921 PM_54_PS_RN_BA_220921 PM_55_PS_RN_BA_220921
2720001.6 4422486.2 5487033.2 5055143.1 5641476.8
PM_56_PS_RN_BA_220921 PM_57_PS_RN_BA_220921 PM_58_PS_RN_BA_220921
3246197.5 3034367.1 3228094.7
The following is quoted directly from https://github.com/hbctraining/DGE_workshop/tree/master/lessons.
Step 2: estimate gene-wide dispersion. The DESeq2 dispersion estimates are inversely related to the mean and directly related to variance. Based on this relationship, the dispersion is higher for small mean counts and lower for large mean counts. Therefore, the dispersion estimates reflect the variance in gene expression for a given mean value. For low mean counts, the variance estimates have a much larger spread; therefore, the dispersion estimates will differ much more between genes with small means. With only a few (3-6) replicates per group, the estimates of variation for each gene are often unreliable (due to the large differences in dispersion for genes with similar means). To address this problem, DESeq2 shares information across genes to generate more accurate estimates of variation based on the mean expression level of the gene using a method called ‘shrinkage’. DESeq2 assumes that genes with similar expression levels have similar dispersion.
Step 3: fit curve to gene-wise dispersion estimates. This curve is displayed as a red line in the figure below, which plots the estimate for the expected dispersion value for genes of a given expression strength. Each black dot is a gene with an associated mean expression level and maximum likelihood estimation (MLE) of the dispersion (Step 1).
Step 4: shrink gene-wise dispersion estimates toward the values predicted by the curve. This shrinkage method is particularly important to reduce false positives in the differential expression analysis. Genes with low dispersion estimates are shrunken towards the curve, and the more accurate, higher shrunken values are output for fitting of the model and differential expression testing. Dispersion estimates that are slightly above the curve are also shrunk toward the curve for better dispersion estimation; however, genes with extremely high dispersion values are not. This is due to the likelihood that the gene does not follow the modeling assumptions and has higher variability than others for biological or technical reasons. Shrinking the values toward the curve could result in false positives, so these values are not shrunken.
Plot dispersion estimates.
plotDispEsts(dds_PB)
plotDispEsts(dds_PL)
plotDispEsts(dds_PM)
plotDispEsts(dds_PS)
Shrunken log2 foldchanges (LFC).
NOTE: Shrinking the log2 fold changes will not change the total number of genes that are identified as significantly differentially expressed. The shrinkage of fold change is to help with downstream assessment of results. For example, if you wanted to subset your significant genes based on fold change for further evaluation, you may want to use shrunken values. Additionally, for functional analysis tools such as GSEA which require fold change values as input you would want to provide shrunken values. In the most recent versions of DESeq2, the shrinkage of LFC estimates is not performed by default. This means that the log2 foldchanges would be the same as those calculated by: log2(normalized_counts_group1 / normalized_counts_group2)
Lists all the coefficient names for apeglm LFC shrinkage.
resultsNames(dds_PB)
[1] "Intercept" "THC_positive_vs_negative"
resultsNames(dds_PL)
[1] "Intercept" "THC_positive_vs_negative"
resultsNames(dds_PM)
[1] "Intercept" "THC_positive_vs_negative"
resultsNames(dds_PS)
[1] "Intercept" "THC_positive_vs_negative"
Builds results tables for the unshrunken LFC and apeglm LFC shrinkage.
contrast=c("THC", "positive", "negative")
# PB
res_table_unshrunken_PB=results(dds_PB, contrast=contrast, alpha=0.1)
res_table_apeglm_PB=lfcShrink(dds_PB, coef="THC_positive_vs_negative",
res=res_table_unshrunken_PB, type="apeglm")
using 'apeglm' for LFC shrinkage. If used in published research, please cite:
Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
sequence count data: removing the noise and preserving large differences.
Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
Warning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precision
# PL
res_table_unshrunken_PL=results(dds_PL, contrast=contrast, alpha=0.1)
res_table_apeglm_PL=lfcShrink(dds_PL, coef="THC_positive_vs_negative",
res=res_table_unshrunken_PL, type="apeglm")
using 'apeglm' for LFC shrinkage. If used in published research, please cite:
Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
sequence count data: removing the noise and preserving large differences.
Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
# PM
res_table_unshrunken_PM=results(dds_PM, contrast=contrast, alpha=0.1)
res_table_apeglm_PM=lfcShrink(dds_PM, coef="THC_positive_vs_negative",
res=res_table_unshrunken_PM, type="apeglm")
using 'apeglm' for LFC shrinkage. If used in published research, please cite:
Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
sequence count data: removing the noise and preserving large differences.
Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
Warning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precision
# PS
res_table_unshrunken_PS=results(dds_PS, contrast=contrast, alpha=0.1)
res_table_apeglm_PS=lfcShrink(dds_PS, coef="THC_positive_vs_negative",
res=res_table_unshrunken_PS, type="apeglm")
using 'apeglm' for LFC shrinkage. If used in published research, please cite:
Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
sequence count data: removing the noise and preserving large differences.
Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
Warning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precisionWarning: the line search routine failed, possibly due to insufficient numeric precision
Creates the MA plot for the apeglm results tables.
options(ggrepel.max.overlaps=Inf)
# PB
# ggmaplot(res_table_unshrunken_PB, main=expression("positive" %->% "negative"),
# fdr=0.05, fc=1.5, size=0.4,
# genenames=as.vector(res_table_unshrunken_PB$rownames),
# legend="top", top=10, select.top.method="padj",
# font.label=c("bold", 11), label.rectangle=TRUE,
# font.legend="bold",
# font.main="bold",
# ggtheme=ggplot2::theme_minimal())
maplot1=ggmaplot(res_table_apeglm_PB,
# main=expression("positive" %->% "negative"),
main="Brain, THC Positive vs Negative",
fdr=0.1, fc=1.5, size=0.4,
legend="top", select.top.method="padj",
top=0,
# font.label=c("bold", 8), label.rectangle=TRUE,
# font.legend="bold",
font.main=c("bold", 10),
ggtheme=ggplot2::theme_minimal())
# PL
# ggmaplot(res_table_unshrunken_PL, main=expression("positive" %->% "negative"),
# fdr=0.05, fc=1.5, size=0.4,
# genenames=as.vector(res_table_unshrunken_PL$rownames),
# legend="top", top=10, select.top.method="padj",
# font.label=c("bold", 11), label.rectangle=TRUE,
# font.legend="bold",
# font.main="bold",
# ggtheme=ggplot2::theme_minimal())
maplot2=ggmaplot(res_table_apeglm_PL,
# main=expression("positive" %->% "negative"),
main="Lung, THC Positive vs Negative",
fdr=0.1, fc=1.5, size=0.4,
# genenames=as.vector(res_table_apeglm_PL$rownames),
legend="top", select.top.method="padj",
top=0,
# font.label=c("bold", 8), label.rectangle=TRUE,
# font.legend="bold",
font.main=c("bold", 10),
ggtheme=ggplot2::theme_minimal())
# PM
# ggmaplot(res_table_unshrunken_PM, main=expression("positive" %->% "negative"),
# fdr=0.05, fc=1.5, size=0.4,
# genenames=as.vector(res_table_unshrunken_PM$rownames),
# legend="top", top=10, select.top.method="padj",
# font.label=c("bold", 11), label.rectangle=TRUE,
# font.legend="bold",
# font.main="bold",
# ggtheme=ggplot2::theme_minimal())
maplot3=ggmaplot(res_table_apeglm_PM,
# main=expression("positive" %->% "negative"),
main="Muscle, THC Positive vs Negative",
fdr=0.1, fc=1.5, size=0.4,
# genenames=as.vector(res_table_apeglm_PM$rownames),
legend="top", select.top.method="padj",
top=0,
# font.label=c("bold", 8), label.rectangle=TRUE,
# font.legend="bold",
font.main=c("bold", 10),
ggtheme=ggplot2::theme_minimal())
# PS
# ggmaplot(res_table_unshrunken_PS, main=expression("positive" %->% "negative"),
# fdr=0.05, fc=1.5, size=0.4,
# genenames=as.vector(res_table_unshrunken_PS$rownames),
# legend="top", top=10, select.top.method="padj",
# font.label=c("bold", 11), label.rectangle=TRUE,
# font.legend="bold",
# font.main="bold",
# ggtheme=ggplot2::theme_minimal())
maplot4=ggmaplot(res_table_apeglm_PS,
# main=expression("positive" %->% "negative"),
main="Blood, THC Positive vs Negative",
fdr=0.1, fc=1.5, size=0.4,
# genenames=as.vector(res_table_apeglm_PS$rownames),
legend="top", select.top.method="padj",
top=0,
# font.label=c("bold", 8), label.rectangle=TRUE,
# font.legend="bold",
font.main=c("bold", 10),
ggtheme=ggplot2::theme_minimal())
Creating a figure of all THC MA plots.
ma_THC_patchwork=maplot1 + maplot2 + maplot3 + maplot4 + plot_layout(ncol=2)
ma_THC_patchwork + plot_annotation(tag_levels="A")
ggsave("PM_ma_THC_plots.png", width=7, height=10, unit="in", dpi=320)
Renames the results table; downstream analysis will be using the apeglm shrunken LFC values.
res_table_PB=res_table_apeglm_PB
res_table_PL=res_table_apeglm_PL
res_table_PM=res_table_apeglm_PM
res_table_PS=res_table_apeglm_PS
Results exploration.
# PB
class(res_table_PB)
[1] "DESeqResults"
attr(,"package")
[1] "DESeq2"
mcols(res_table_PB, use.names=TRUE)
DataFrame with 5 rows and 2 columns
type description
<character> <character>
baseMean intermediate mean of normalized c..
log2FoldChange results log2 fold change (MA..
lfcSE results posterior SD: THC po..
pvalue results Wald test p-value: T..
padj results BH adjusted p-values
res_table_PB %>% data.frame()
# PL
class(res_table_PL)
[1] "DESeqResults"
attr(,"package")
[1] "DESeq2"
mcols(res_table_PL, use.names=TRUE)
DataFrame with 5 rows and 2 columns
type description
<character> <character>
baseMean intermediate mean of normalized c..
log2FoldChange results log2 fold change (MA..
lfcSE results posterior SD: THC po..
pvalue results Wald test p-value: T..
padj results BH adjusted p-values
res_table_PL %>% data.frame()
# PM
class(res_table_PM)
[1] "DESeqResults"
attr(,"package")
[1] "DESeq2"
mcols(res_table_PM, use.names=TRUE)
DataFrame with 5 rows and 2 columns
type description
<character> <character>
baseMean intermediate mean of normalized c..
log2FoldChange results log2 fold change (MA..
lfcSE results posterior SD: THC po..
pvalue results Wald test p-value: T..
padj results BH adjusted p-values
res_table_PM %>% data.frame()
# PS
class(res_table_PS)
[1] "DESeqResults"
attr(,"package")
[1] "DESeq2"
mcols(res_table_PS, use.names=TRUE)
DataFrame with 5 rows and 2 columns
type description
<character> <character>
baseMean intermediate mean of normalized c..
log2FoldChange results log2 fold change (MA..
lfcSE results posterior SD: THC po..
pvalue results Wald test p-value: T..
padj results BH adjusted p-values
res_table_PS %>% data.frame()
Summarizing results.
summary(res_table_PB)
out of 50232 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 0, 0%
LFC < 0 (down) : 1, 0.002%
outliers [1] : 0, 0%
low counts [2] : 10, 0.02%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
summary(res_table_PL)
out of 60183 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 22, 0.037%
LFC < 0 (down) : 2, 0.0033%
outliers [1] : 0, 0%
low counts [2] : 38420, 64%
(mean count < 8)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
summary(res_table_PM)
out of 55776 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 3, 0.0054%
LFC < 0 (down) : 1, 0.0018%
outliers [1] : 0, 0%
low counts [2] : 17, 0.03%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
summary(res_table_PS)
out of 60084 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 0, 0%
LFC < 0 (down) : 0, 0%
outliers [1] : 0, 0%
low counts [2] : 1, 0.0017%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
Extracting significant differentially expressed genes.
## Set thresholds.
padj.cutoff=0.1
lfc.cutoff=0.58
# PB
## Convert the results table to a tibble.
res_table_tb_PB=res_table_PB %>%
data.frame() %>%
rownames_to_column(var='gene') %>%
as_tibble()
## Subset the tibble use pre-defined thresholds and then show the table.
sig_PB=res_table_tb_PB %>%
dplyr::filter(padj<padj.cutoff & abs(log2FoldChange) > lfc.cutoff) %>%
dplyr::arrange(padj)
sig_PB %>%
count()
sig_PB %>%
dplyr::filter(log2FoldChange>0) %>%
count()
sig_PB %>%
dplyr::filter(log2FoldChange<0) %>%
count()
## Write the full and subsetted results to .txt and .csv files.
write.table(res_table_tb_PB, "gencode_Output_Files/PB_THC_genes.txt", quote=FALSE, row.names=FALSE)
write.csv(res_table_tb_PB, "gencode_Output_Files/PB_THC_genes.csv", quote=FALSE, row.names=FALSE)
write.table(sig_PB, "gencode_Output_Files/PB_THC_sig_DEGs.txt", quote=FALSE, row.names=FALSE) # IPA requires a .txt file
write.csv(sig_PB, "gencode_Output_Files/PB_THC_sig_DEGs.csv", quote=FALSE, row.names=FALSE) # .csv files work better for R
# PL
## Convert the results table to a tibble.
res_table_tb_PL=res_table_PL %>%
data.frame() %>%
rownames_to_column(var='gene') %>%
as_tibble()
## Subset the tibble use pre-defined thresholds and then show the table.
sig_PL=res_table_tb_PL %>%
dplyr::filter(padj<padj.cutoff & abs(log2FoldChange) > lfc.cutoff) %>%
dplyr::arrange(padj)
sig_PL %>%
count()
sig_PL %>%
dplyr::filter(log2FoldChange>0) %>%
count()
sig_PL %>%
dplyr::filter(log2FoldChange<0) %>%
count()
## Write the full and subsetted results to .txt and .csv files.
write.table(res_table_tb_PL, "gencode_Output_Files/PL_THC_genes.txt", quote=FALSE, row.names=FALSE)
write.csv(res_table_tb_PL, "gencode_Output_Files/PL_THC_genes.csv", quote=FALSE, row.names=FALSE)
write.table(sig_PL, "gencode_Output_Files/PL_THC_sig_DEGs.txt", quote=FALSE, row.names=FALSE) # IPA requires a .txt file
write.csv(sig_PL, "gencode_Output_Files/PL_THC_sig_DEGs.csv", quote=FALSE, row.names=FALSE) # .csv files work better for R
# PM
## Convert the results table to a tibble.
res_table_tb_PM=res_table_PM %>%
data.frame() %>%
rownames_to_column(var='gene') %>%
as_tibble()
## Subset the tibble use pre-defined thresholds and then show the table.
sig_PM=res_table_tb_PM %>%
dplyr::filter(padj<padj.cutoff & abs(log2FoldChange) > lfc.cutoff) %>%
dplyr::arrange(padj)
sig_PM %>%
count()
sig_PM %>%
dplyr::filter(log2FoldChange>0) %>%
count()
sig_PM %>%
dplyr::filter(log2FoldChange<0) %>%
count()
## Write the full and subsetted results to .txt and .csv files.
write.table(res_table_tb_PM, "gencode_Output_Files/PM_THC_genes.txt", quote=FALSE, row.names=FALSE)
write.csv(res_table_tb_PM, "gencode_Output_Files/PM_THC_genes.csv", quote=FALSE, row.names=FALSE)
write.table(sig_PM, "gencode_Output_Files/PM_THC_sig_DEGs.txt", quote=FALSE, row.names=FALSE) # IPA requires a .txt file
write.csv(sig_PM, "gencode_Output_Files/PM_THC_sig_DEGs.csv", quote=FALSE, row.names=FALSE) # .csv files work better for R
# PS
## Convert the results table to a tibble.
res_table_tb_PS=res_table_PS %>%
data.frame() %>%
rownames_to_column(var='gene') %>%
as_tibble()
## Subset the tibble use pre-defined thresholds and then show the table.
sig_PS=res_table_tb_PS %>%
dplyr::filter(padj<padj.cutoff & abs(log2FoldChange) > lfc.cutoff) %>%
dplyr::arrange(padj)
sig_PS %>%
count()
sig_PS %>%
dplyr::filter(log2FoldChange>0) %>%
count()
sig_PS %>%
dplyr::filter(log2FoldChange<0) %>%
count()
## Write the full and subsetted results to .txt and .csv files.
write.table(res_table_tb_PS, "gencode_Output_Files/PS_THC_genes.txt", quote=FALSE, row.names=FALSE)
write.csv(res_table_tb_PS, "gencode_Output_Files/PS_THC_genes.csv", quote=FALSE, row.names=FALSE)
write.table(sig_PS, "gencode_Output_Files/PS_THC_sig_DEGs.txt", quote=FALSE, row.names=FALSE) # IPA requires a .txt file
write.csv(sig_PS, "gencode_Output_Files/PS_THC_sig_DEGs.csv", quote=FALSE, row.names=FALSE) # .csv files work better for R
Reads in the normalized expression data table.
normalizedcounts_PB=read.table(file="gencode_Output_Files/PB_normalizedcounts.txt", header=TRUE, row.names=1)
normalizedcounts_PL=read.table(file="gencode_Output_Files/PL_normalizedcounts.txt", header=TRUE, row.names=1)
normalizedcounts_PM=read.table(file="gencode_Output_Files/PM_normalizedcounts.txt", header=TRUE, row.names=1)
normalizedcounts_PS=read.table(file="gencode_Output_Files/PS_normalizedcounts.txt", header=TRUE, row.names=1)
Create tibbles including row names from metadata and normalized expression data.
# PB
meta_PB=metadata_PB %>%
rownames_to_column(var="samplename") %>%
as_tibble()
normalized_counts_PB=normalizedcounts_PB %>%
rownames_to_column(var="gene") %>%
as_tibble()
# PL
meta_PL=metadata_PL %>%
rownames_to_column(var="samplename") %>%
as_tibble()
normalized_counts_PL=normalizedcounts_PL %>%
rownames_to_column(var="gene") %>%
as_tibble()
# PM
meta_PM=metadata_PM %>%
rownames_to_column(var="samplename") %>%
as_tibble()
normalized_counts_PM=normalizedcounts_PM %>%
rownames_to_column(var="gene") %>%
as_tibble()
# PS
meta_PS=metadata_PS %>%
rownames_to_column(var="samplename") %>%
as_tibble()
normalized_counts_PS=normalizedcounts_PS %>%
rownames_to_column(var="gene") %>%
as_tibble()
Using ggplot2 to plot multiple genes (top 20).
# # PB
# ## Order results by padj values.
# top20_sig_genes_PB=res_table_tb_PB %>%
# arrange(padj) %>% ## Arrange rows by padj values.
# pull(gene) %>% ## Extract character vector of ordered genes.
# head(n=20) ## Extract the top 20 genes.
# ## Normalized counts for top 20 significant genes.
# top20_sig_norm_PB=normalized_counts_PB %>%
# dplyr::filter(gene %in% top20_sig_genes_PB)
# ## Gathering the columns to have normalized to a single column.
# gathered_top20_sig_PB=top20_sig_norm_PB %>%
# gather(colnames(top20_sig_norm_PB)[2:50], key="samplename", value="normalized_counts")
# ## Combine metadata with melted normalized counts data.
# gathered_top20_sig_PB=inner_join(meta_PB, gathered_top20_sig_PB)
# ## Plot using ggplot2.
# ggplot()+
# geom_point(data=gathered_top20_sig_PB, aes(x=gene, y=normalized_counts, fill=THC),
# alpha=0.75, color="black", pch=21, size=2)+
# scale_fill_manual(limits=c("positive", "negative"),
# labels=c("Positive", "Negative"),
# values=as.vector(polychrome(2)))+
# scale_y_log10()+
# theme_minimal()+
# theme(axis.text.y=element_text(color="black", face="bold", size=rel(1)),
# axis.text.x=element_text(angle=45, hjust=1, color="black", size=rel(0.8)),
# # axis.text.x=element_text(color="black", face="bold", size=rel(1)),
# # axis.title.x=element_text(color="black", face="bold", size=rel(1)),
# axis.title.x=element_blank(),
# axis.title.y=element_text(color="black", face="bold", size=rel(1)),
# panel.border=element_rect(color="black", linewidth=1, fill=NA),
# strip.background=element_rect(color="black", linewidth=1),
# strip.text=element_text(color="black", face="bold", size=rel(0.8)),
# axis.line=element_blank(),
# axis.ticks=element_line(linewidth=1),
# legend.title=element_text(color="black", face="bold", size=rel(1)),
# # legend.title=element_blank(),
# legend.text=element_text(color="black", face="bold", size=rel(0.8)),
# legend.position="bottom",
# # plot.title=element_text(color="black", face="bold", size=rel(1)))+
# plot.title=element_blank())+
# xlab("Genes")+
# ylab("log10 Normalized Counts")+
# labs(title="Significant Genes in Brain THC Positive vs Negative", fill="THC")
# PL
## Order results by padj values.
top20_sig_genes_PL=res_table_tb_PL %>%
arrange(padj) %>% ## Arrange rows by padj values.
pull(gene) %>% ## Extract character vector of ordered genes.
head(n=22) ## Extract the top 20 genes.
## Normalized counts for top 20 significant genes.
top20_sig_norm_PL=normalized_counts_PL %>%
dplyr::filter(gene %in% top20_sig_genes_PL)
## Gathering the columns to have normalized to a single column.
gathered_top20_sig_PL=top20_sig_norm_PL %>%
gather(colnames(top20_sig_norm_PL)[2:53], key="samplename", value="normalized_counts")
## Combine metadata with melted normalized counts data.
gathered_top20_sig_PL=inner_join(meta_PL, gathered_top20_sig_PL)
Joining with `by = join_by(samplename)`
## Plot using ggplot2.
plot1=ggplot()+
geom_point(data=gathered_top20_sig_PL, aes(x=gene, y=normalized_counts, fill=THC, shape=THC),
alpha=0.75, color="black", size=2)+
scale_fill_manual(limits=c("positive", "negative"),
labels=c("Positive", "Negative"),
values=as.vector(polychrome(2)))+
scale_shape_manual(limits=c("positive", "negative"),
labels=c("Positive", "Negative"),
values=c(21, 24))+
scale_y_log10()+
theme_minimal()+
theme(axis.text.y=element_text(color="black", face="bold", size=rel(1)),
axis.text.x=element_text(angle=45, hjust=1, color="black", size=rel(0.8)),
# axis.text.x=element_text(color="black", face="bold", size=rel(1)),
# axis.title.x=element_text(color="black", face="bold", size=rel(1)),
axis.title.x=element_blank(),
axis.title.y=element_text(color="black", face="bold", size=rel(1)),
panel.border=element_rect(color="black", linewidth=1, fill=NA),
strip.background=element_rect(color="black", linewidth=1),
strip.text=element_text(color="black", face="bold", size=rel(0.8)),
axis.line=element_blank(),
axis.ticks=element_line(linewidth=1),
legend.title=element_text(color="black", face="bold", size=rel(1)),
# legend.title=element_blank(),
legend.text=element_text(color="black", face="bold", size=rel(0.8)),
legend.position="bottom",
# plot.title=element_text(color="black", face="bold", size=rel(1)))+
plot.title=element_blank())+
xlab("Genes")+
ylab("log10 Normalized Counts")+
labs(title="Significant Genes in Lung THC Positive vs Negative", fill="THC")
# PM
## Order results by padj values.
top20_sig_genes_PM=res_table_tb_PM %>%
arrange(padj) %>% ## Arrange rows by padj values.
pull(gene) %>% ## Extract character vector of ordered genes.
head(n=4) ## Extract the top 20 genes.
## Normalized counts for top 20 significant genes.
top20_sig_norm_PM=normalized_counts_PM %>%
dplyr::filter(gene %in% top20_sig_genes_PM)
## Gathering the columns to have normalized to a single column.
gathered_top20_sig_PM=top20_sig_norm_PM %>%
gather(colnames(top20_sig_norm_PM)[2:54], key="samplename", value="normalized_counts")
## Combine metadata with melted normalized counts data.
gathered_top20_sig_PM=inner_join(meta_PM, gathered_top20_sig_PM)
Joining with `by = join_by(samplename)`
## Plot using ggplot2.
plot2=ggplot()+
geom_point(data=gathered_top20_sig_PM, aes(x=gene, y=normalized_counts, fill=THC, shape=THC),
alpha=0.75, color="black", size=2)+
scale_fill_manual(limits=c("positive", "negative"),
labels=c("Positive", "Negative"),
values=as.vector(polychrome(2)))+
scale_shape_manual(limits=c("positive", "negative"),
labels=c("Positive", "Negative"),
values=c(21, 24))+
scale_y_log10()+
theme_minimal()+
theme(axis.text.y=element_text(color="black", face="bold", size=rel(1)),
axis.text.x=element_text(angle=45, hjust=1, color="black", size=rel(0.8)),
# axis.text.x=element_text(color="black", face="bold", size=rel(1)),
# axis.title.x=element_text(color="black", face="bold", size=rel(1)),
axis.title.x=element_blank(),
axis.title.y=element_text(color="black", face="bold", size=rel(1)),
panel.border=element_rect(color="black", linewidth=1, fill=NA),
strip.background=element_rect(color="black", linewidth=1),
strip.text=element_text(color="black", face="bold", size=rel(0.8)),
axis.line=element_blank(),
axis.ticks=element_line(linewidth=1),
legend.title=element_text(color="black", face="bold", size=rel(1)),
# legend.title=element_blank(),
legend.text=element_text(color="black", face="bold", size=rel(0.8)),
legend.position="bottom",
# plot.title=element_text(color="black", face="bold", size=rel(1)))+
plot.title=element_blank())+
xlab("Genes")+
ylab("log10 Normalized Counts")+
labs(title="Significant Genes in Muscle THC Positive vs Negative", fill="THC")
# # PS
# ## Order results by padj values.
# top20_sig_genes_PS=res_table_tb_PS %>%
# arrange(padj) %>% ## Arrange rows by padj values.
# pull(gene) %>% ## Extract character vector of ordered genes.
# head(n=20) ## Extract the top 20 genes.
# ## Normalized counts for top 20 significant genes.
# top20_sig_norm_PS=normalized_counts_PS %>%
# dplyr::filter(gene %in% top20_sig_genes_PS)
# ## Gathering the columns to have normalized to a single column.
# gathered_top20_sig_PS=top20_sig_norm_PS %>%
# gather(colnames(top20_sig_norm_PS)[2:29], key="samplename", value="normalized_counts")
# ## Combine metadata with melted normalized counts data.
# gathered_top20_sig_PS=inner_join(meta_PS, gathered_top20_sig_PS)
# ## Plot using ggplot2.
# ggplot()+
# geom_point(data=gathered_top20_sig_PS, aes(x=gene, y=normalized_counts, fill=THC),
# alpha=0.75, color="black", pch=21, size=2)+
# scale_fill_manual(limits=c("positive", "negative"),
# labels=c("Positive", "Negative"),
# values=as.vector(polychrome(2)))+
# scale_y_log10()+
# theme_minimal()+
# theme(axis.text.y=element_text(color="black", face="bold", size=rel(1)),
# axis.text.x=element_text(angle=45, hjust=1, color="black", size=rel(0.8)),
# # axis.text.x=element_text(color="black", face="bold", size=rel(1)),
# # axis.title.x=element_text(color="black", face="bold", size=rel(1)),
# axis.title.x=element_blank(),
# axis.title.y=element_text(color="black", face="bold", size=rel(1)),
# panel.border=element_rect(color="black", linewidth=1, fill=NA),
# strip.background=element_rect(color="black", linewidth=1),
# strip.text=element_text(color="black", face="bold", size=rel(0.8)),
# axis.line=element_blank(),
# axis.ticks=element_line(linewidth=1),
# legend.title=element_text(color="black", face="bold", size=rel(1)),
# # legend.title=element_blank(),
# legend.text=element_text(color="black", face="bold", size=rel(0.8)),
# legend.position="bottom",
# # plot.title=element_text(color="black", face="bold", size=rel(1)))+
# plot.title=element_blank())+
# xlab("Genes")+
# ylab("log10 Normalized Counts")+
# labs(title="Significant Genes in Blood THC Positive vs Negative", fill="THC")
Creating a figure of all significant genes plots.
top_sig_patchwork=plot1 + plot2 + plot_layout(nrow=2)
top_sig_patchwork + plot_annotation(tag_levels="A") + plot_layout(guides="collect") &
theme(legend.position='bottom')
ggsave("PM_top_sig_THC_plots.png", width=7, height=10, unit="in", dpi=320)
Generating a heatmap.
# # PB
# ## Extract normalized expression for significant genes from the samples and set the gene column (1) to row names.
# norm_sig_PB=normalized_counts_PB[,c(1,2:54)] %>%
# dplyr::filter(gene %in% sig_PB$gene) %>%
# data.frame(check.names=FALSE) %>%
# column_to_rownames(var="gene")
# ## Annotate the heatmap.
# annotation_PB=meta_PB %>%
# dplyr::select(samplename, THC) %>%
# data.frame(row.names="samplename")
# ## Set a color palette.
# heat_colors=brewer.pal(6, "YlOrRd")
# ## Run pheatmap.
# # pheatmap::pheatmap(norm_sig_PB,
# # color=heat_colors,
# # cluster_rows=TRUE,
# # show_rownames=FALSE,
# # annotation_col=annotation_PB,
# # border_color=NA,
# # fontsize=9,
# # scale="row", ## Z-scores are plotted rather than actual normalized count values.
# # fontsize_row=10,
# # height=20,
# # # annotation_colors=list(group=c(BP="", RN="", RP="")), ## The annotation colors.
# # annotation_legend=TRUE)
#
# # PL
# ## Extract normalized expression for significant genes from the samples and set the gene column (1) to row names.
# norm_sig_PL=normalized_counts_PL[,c(1,2:52)] %>%
# dplyr::filter(gene %in% sig_PL$gene) %>%
# data.frame(check.names=FALSE) %>%
# column_to_rownames(var="gene")
# ## Annotate the heatmap.
# annotation_PL=meta_PL %>%
# dplyr::select(samplename, THC) %>%
# data.frame(row.names="samplename")
# ## Set a color palette.
# heat_colors=brewer.pal(6, "YlOrRd")
# ## Run pheatmap.
# pheatmap::pheatmap(norm_sig_PL,
# color=heat_colors,
# cluster_rows=TRUE,
# show_rownames=FALSE,
# annotation_col=annotation_PL,
# border_color=NA,
# fontsize=9,
# scale="row", ## Z-scores are plotted rather than actual normalized count values.
# fontsize_row=10,
# height=20,
# # annotation_colors=list(group=c(BP="", RN="", RP="")), ## The annotation colors.
# annotation_legend=TRUE)
#
# # PM
# ## Extract normalized expression for significant genes from the samples and set the gene column (1) to row names.
# norm_sig_PM=normalized_counts_PM[,c(1,2:54)] %>%
# dplyr::filter(gene %in% sig_PM$gene) %>%
# data.frame(check.names=FALSE) %>%
# column_to_rownames(var="gene")
# ## Annotate the heatmap.
# annotation_PM=meta_PM %>%
# dplyr::select(samplename, THC) %>%
# data.frame(row.names="samplename")
# ## Set a color palette.
# heat_colors=brewer.pal(6, "YlOrRd")
# ## Run pheatmap.
# pheatmap::pheatmap(norm_sig_PM,
# color=heat_colors,
# cluster_rows=TRUE,
# show_rownames=FALSE,
# annotation_col=annotation_PM,
# border_color=NA,
# fontsize=9,
# scale="row", ## Z-scores are plotted rather than actual normalized count values.
# fontsize_row=10,
# height=20,
# # annotation_colors=list(group=c(BP="", RN="", RP="")), ## The annotation colors.
# annotation_legend=TRUE)
#
# # PS
# ## Extract normalized expression for significant genes from the samples and set the gene column (1) to row names.
# norm_sig_PS=normalized_counts_PS[,c(1,2:44)] %>%
# dplyr::filter(gene %in% sig_PS$gene) %>%
# data.frame(check.names=FALSE) %>%
# column_to_rownames(var="gene")
# ## Annotate the heatmap.
# annotation_PS=meta_PS %>%
# dplyr::select(samplename, THC) %>%
# data.frame(row.names="samplename")
# ## Set a color palette.
# heat_colors=brewer.pal(6, "YlOrRd")
# ## Run pheatmap.
# # pheatmap::pheatmap(norm_sig_PS,
# # color=heat_colors,
# # cluster_rows=TRUE,
# # show_rownames=FALSE,
# # annotation_col=annotation_PS,
# # border_color=NA,
# # fontsize=9,
# # scale="row", ## Z-scores are plotted rather than actual normalized count values.
# # fontsize_row=10,
# # height=20,
# # # annotation_colors=list(group=c(BP="", RN="", RP="")), ## The annotation colors.
# # annotation_legend=TRUE)
Generating a volcano plot.
# # # PB
# # ## Obtain logical vector where TRUE values denote padj values < 0.1 and fold change > 1.5 in either direction.
# # res_table_tb_vp_PB=res_table_tb_PB %>%
# # dplyr::mutate(threshold=padj<0.1 & abs(log2FoldChange)>=0.58)
# # ## Plots the volcano plot.
# # ggplot()+
# # geom_point(data=res_table_tb_vp_PB, aes(x=log2FoldChange, y=-log10(padj), color=threshold))+
# # ggtitle("THC:negative versus THC:positive Expression")+
# # xlab("log2 fold change")+
# # ylab("-log10 adjusted p-value")+
# # # scale_color_viridis(option="turbo", discrete=TRUE)+
# # theme_minimal()+
# # theme(axis.text.y=element_text(color="black", face="bold", size=rel(1.2)),
# # axis.text.x=element_text(color="black", face="bold", size=rel(1.2)),
# # axis.title.x=element_text(color="black", face="bold", size=rel(1)),
# # axis.title.y=element_text(color="black", face="bold", size=rel(1)),
# # panel.border=element_rect(color="black", linewidth=1, fill=NA),
# # strip.background=element_rect(color="black", linewidth=1),
# # strip.text=element_text(color="black", face="bold", size=rel(0.8)),
# # axis.line=element_blank(),
# # axis.ticks=element_line(linewidth=1),
# # legend.position="none",
# # legend.title=element_text(color="black", face="bold", size=rel(1)),
# # legend.text=element_text(color="black", face="bold", size=rel(0.8)),
# # plot.title=element_text(color="black", face="bold", size=rel(1)))
# #
# # ## Label the top 10 genes with the lowest padj (conditional on max.overlaps value).
# # ## Create a column to indicate which genes to label.
# # res_table_tb_vp_PB=res_table_tb_vp_PB %>% dplyr::arrange(padj) %>% dplyr::mutate(genelabels="")
# # res_table_tb_vp_PB$genelabels[1:10]=res_table_tb_vp_PB$gene[1:10]
# # ## Plots the volcano plot with an additional layer for geom_text_repel().
# # ggplot(res_table_tb_vp_PB, aes(x=log2FoldChange, y=-log10(padj)))+
# # geom_point(aes(color=threshold))+
# # geom_text_repel(aes(label=genelabels), force=3, min.segment.length=0.5, box.padding=0.5)+
# # ggtitle("THC:negative versus THC:positive Expression")+
# # xlab("log2 fold change")+
# # ylab("-log10 adjusted p-value")+
# # # scale_color_viridis(option="turbo", discrete=TRUE)+
# # theme_minimal()+
# # theme(axis.text.y=element_text(color="black", face="bold", size=rel(1.2)),
# # axis.text.x=element_text(color="black", face="bold", size=rel(1.2)),
# # axis.title.x=element_text(color="black", face="bold", size=rel(1)),
# # axis.title.y=element_text(color="black", face="bold", size=rel(1)),
# # panel.border=element_rect(color="black", linewidth=1, fill=NA),
# # strip.background=element_rect(color="black", linewidth=1),
# # strip.text=element_text(color="black", face="bold", size=rel(0.8)),
# # axis.line=element_blank(),
# # axis.ticks=element_line(linewidth=1),
# # legend.position="none",
# # legend.title=element_text(color="black", face="bold", size=rel(1)),
# # legend.text=element_text(color="black", face="bold", size=rel(0.8)),
# # plot.title=element_text(color="black", face="bold", size=rel(1)))
#
# # PL
# ## Obtain logical vector where TRUE values denote padj values < 0.1 and fold change > 1.5 in either direction.
# res_table_tb_vp_PL=res_table_tb_PL %>%
# dplyr::mutate(threshold=padj<0.1 & abs(log2FoldChange)>=0.58)
# ## Plots the volcano plot.
# ggplot()+
# geom_point(data=res_table_tb_vp_PL, aes(x=log2FoldChange, y=-log10(padj), color=threshold))+
# ggtitle("THC:negative versus THC:positive Expression")+
# xlab("log2 fold change")+
# ylab("-log10 adjusted p-value")+
# # scale_color_viridis(option="turbo", discrete=TRUE)+
# theme_minimal()+
# theme(axis.text.y=element_text(color="black", face="bold", size=rel(1.2)),
# axis.text.x=element_text(color="black", face="bold", size=rel(1.2)),
# axis.title.x=element_text(color="black", face="bold", size=rel(1)),
# axis.title.y=element_text(color="black", face="bold", size=rel(1)),
# panel.border=element_rect(color="black", linewidth=1, fill=NA),
# strip.background=element_rect(color="black", linewidth=1),
# strip.text=element_text(color="black", face="bold", size=rel(0.8)),
# axis.line=element_blank(),
# axis.ticks=element_line(linewidth=1),
# legend.position="none",
# legend.title=element_text(color="black", face="bold", size=rel(1)),
# legend.text=element_text(color="black", face="bold", size=rel(0.8)),
# plot.title=element_text(color="black", face="bold", size=rel(1)))
#
# ## Label the top 10 genes with the lowest padj (conditional on max.overlaps value).
# ## Create a column to indicate which genes to label.
# res_table_tb_vp_PL=res_table_tb_vp_PL %>% dplyr::arrange(padj) %>% dplyr::mutate(genelabels="")
# res_table_tb_vp_PL$genelabels[1:29]=res_table_tb_vp_PL$gene[1:29]
# ## Plots the volcano plot with an additional layer for geom_text_repel().
# ggplot(res_table_tb_vp_PL, aes(x=log2FoldChange, y=-log10(padj)))+
# geom_point(aes(color=threshold))+
# geom_text_repel(aes(label=genelabels), force=3, min.segment.length=0.5, box.padding=0.5)+
# ggtitle("THC:negative versus THC:positive Expression")+
# xlab("log2 fold change")+
# ylab("-log10 adjusted p-value")+
# # scale_color_viridis(option="turbo", discrete=TRUE)+
# theme_minimal()+
# theme(axis.text.y=element_text(color="black", face="bold", size=rel(1.2)),
# axis.text.x=element_text(color="black", face="bold", size=rel(1.2)),
# axis.title.x=element_text(color="black", face="bold", size=rel(1)),
# axis.title.y=element_text(color="black", face="bold", size=rel(1)),
# panel.border=element_rect(color="black", linewidth=1, fill=NA),
# strip.background=element_rect(color="black", linewidth=1),
# strip.text=element_text(color="black", face="bold", size=rel(0.8)),
# axis.line=element_blank(),
# axis.ticks=element_line(linewidth=1),
# legend.position="none",
# legend.title=element_text(color="black", face="bold", size=rel(1)),
# legend.text=element_text(color="black", face="bold", size=rel(0.8)),
# plot.title=element_text(color="black", face="bold", size=rel(1)))
#
# # PM
# ## Obtain logical vector where TRUE values denote padj values < 0.1 and fold change > 1.5 in either direction.
# res_table_tb_vp_PM=res_table_tb_PM %>%
# dplyr::mutate(threshold=padj<0.1 & abs(log2FoldChange)>=0.58)
# ## Plots the volcano plot.
# ggplot()+
# geom_point(data=res_table_tb_vp_PM, aes(x=log2FoldChange, y=-log10(padj), color=threshold))+
# ggtitle("THC:negative versus THC:positive Expression")+
# xlab("log2 fold change")+
# ylab("-log10 adjusted p-value")+
# # scale_color_viridis(option="turbo", discrete=TRUE)+
# theme_minimal()+
# theme(axis.text.y=element_text(color="black", face="bold", size=rel(1.2)),
# axis.text.x=element_text(color="black", face="bold", size=rel(1.2)),
# axis.title.x=element_text(color="black", face="bold", size=rel(1)),
# axis.title.y=element_text(color="black", face="bold", size=rel(1)),
# panel.border=element_rect(color="black", linewidth=1, fill=NA),
# strip.background=element_rect(color="black", linewidth=1),
# strip.text=element_text(color="black", face="bold", size=rel(0.8)),
# axis.line=element_blank(),
# axis.ticks=element_line(linewidth=1),
# legend.position="none",
# legend.title=element_text(color="black", face="bold", size=rel(1)),
# legend.text=element_text(color="black", face="bold", size=rel(0.8)),
# plot.title=element_text(color="black", face="bold", size=rel(1)))
#
# ## Label the top 10 genes with the lowest padj (conditional on max.overlaps value).
# ## Create a column to indicate which genes to label.
# res_table_tb_vp_PM=res_table_tb_vp_PM %>% dplyr::arrange(padj) %>% dplyr::mutate(genelabels="")
# res_table_tb_vp_PM$genelabels[1:4]=res_table_tb_vp_PM$gene[1:4]
# ## Plots the volcano plot with an additional layer for geom_text_repel().
# ggplot(res_table_tb_vp_PM, aes(x=log2FoldChange, y=-log10(padj)))+
# geom_point(aes(color=threshold))+
# geom_text_repel(aes(label=genelabels), force=3, min.segment.length=0.5, box.padding=0.5)+
# ggtitle("THC:negative versus THC:positive Expression")+
# xlab("log2 fold change")+
# ylab("-log10 adjusted p-value")+
# # scale_color_viridis(option="turbo", discrete=TRUE)+
# theme_minimal()+
# theme(axis.text.y=element_text(color="black", face="bold", size=rel(1.2)),
# axis.text.x=element_text(color="black", face="bold", size=rel(1.2)),
# axis.title.x=element_text(color="black", face="bold", size=rel(1)),
# axis.title.y=element_text(color="black", face="bold", size=rel(1)),
# panel.border=element_rect(color="black", linewidth=1, fill=NA),
# strip.background=element_rect(color="black", linewidth=1),
# strip.text=element_text(color="black", face="bold", size=rel(0.8)),
# axis.line=element_blank(),
# axis.ticks=element_line(linewidth=1),
# legend.position="none",
# legend.title=element_text(color="black", face="bold", size=rel(1)),
# legend.text=element_text(color="black", face="bold", size=rel(0.8)),
# plot.title=element_text(color="black", face="bold", size=rel(1)))
#
# # # PS
# # ## Obtain logical vector where TRUE values denote padj values < 0.1 and fold change > 1.5 in either direction.
# # res_table_tb_vp_PS=res_table_tb_PS %>%
# # dplyr::mutate(threshold=padj<0.1 & abs(log2FoldChange)>=0.58)
# # ## Plots the volcano plot.
# # ggplot()+
# # geom_point(data=res_table_tb_vp_PS, aes(x=log2FoldChange, y=-log10(padj), color=threshold))+
# # ggtitle("THC:negative versus THC:positive Expression")+
# # xlab("log2 fold change")+
# # ylab("-log10 adjusted p-value")+
# # # scale_color_viridis(option="turbo", discrete=TRUE)+
# # theme_minimal()+
# # theme(axis.text.y=element_text(color="black", face="bold", size=rel(1.2)),
# # axis.text.x=element_text(color="black", face="bold", size=rel(1.2)),
# # axis.title.x=element_text(color="black", face="bold", size=rel(1)),
# # axis.title.y=element_text(color="black", face="bold", size=rel(1)),
# # panel.border=element_rect(color="black", linewidth=1, fill=NA),
# # strip.background=element_rect(color="black", linewidth=1),
# # strip.text=element_text(color="black", face="bold", size=rel(0.8)),
# # axis.line=element_blank(),
# # axis.ticks=element_line(linewidth=1),
# # legend.position="none",
# # legend.title=element_text(color="black", face="bold", size=rel(1)),
# # legend.text=element_text(color="black", face="bold", size=rel(0.8)),
# # plot.title=element_text(color="black", face="bold", size=rel(1)))
# #
# # ## Label the top 10 genes with the lowest padj (conditional on max.overlaps value).
# # ## Create a column to indicate which genes to label.
# # res_table_tb_vp_PS=res_table_tb_vp_PS %>% dplyr::arrange(padj) %>% dplyr::mutate(genelabels="")
# # res_table_tb_vp_PS$genelabels[1:10]=res_table_tb_vp_PS$gene[1:10]
# # ## Plots the volcano plot with an additional layer for geom_text_repel().
# # ggplot(res_table_tb_vp_PS, aes(x=log2FoldChange, y=-log10(padj)))+
# # geom_point(aes(color=threshold))+
# # geom_text_repel(aes(label=genelabels), force=3, min.segment.length=0.5, box.padding=0.5)+
# # ggtitle("THC:negative versus THC:positive Expression")+
# # xlab("log2 fold change")+
# # ylab("-log10 adjusted p-value")+
# # # scale_color_viridis(option="turbo", discrete=TRUE)+
# # theme_minimal()+
# # theme(axis.text.y=element_text(color="black", face="bold", size=rel(1.2)),
# # axis.text.x=element_text(color="black", face="bold", size=rel(1.2)),
# # axis.title.x=element_text(color="black", face="bold", size=rel(1)),
# # axis.title.y=element_text(color="black", face="bold", size=rel(1)),
# # panel.border=element_rect(color="black", linewidth=1, fill=NA),
# # strip.background=element_rect(color="black", linewidth=1),
# # strip.text=element_text(color="black", face="bold", size=rel(0.8)),
# # axis.line=element_blank(),
# # axis.ticks=element_line(linewidth=1),
# # legend.position="none",
# # legend.title=element_text(color="black", face="bold", size=rel(1)),
# # legend.text=element_text(color="black", face="bold", size=rel(0.8)),
# # plot.title=element_text(color="black", face="bold", size=rel(1)))
Only perform one of the following conversions before continuing on with the rest of the code block. Uncomment/comment out the appropriate conversion code as appropriate. The de-duplication process can result in a considerable number of entries in the results table being removed, so only do this if you absolutely need to add annotation information for downstream processes.
# ## Convert from gene symbol to Ensembl IDs (only do this if the dds results has gene symbols).
# ## Return the IDs for the gene symbols in the DE results.
# idx=grch38$symbol %in% rownames(res_table)
# ids=grch38[idx, ]
# ## The gene names can map to more than one Ensembl ID (some genes change ID over time), so duplicate IDs need to be removed.
# non_duplicates=which(duplicated(ids$symbol) == FALSE)
# ids=ids[non_duplicates, ]
# ## Merge the IDs with the results.
# res_ids=inner_join(sig, ids, by=c("gene"="symbol"))
# PB
## Convert from Ensembl IDs to Entrez IDs (only do this if the dds results has Ensembl IDs).
## Return the IDs for the gene symbols in the DE results
idx_PB=grch38$ensgene %in% rownames(res_table_PB)
ids_PB=grch38[idx_PB, ]
## The Entrez IDs can map to more than one Ensembl ID (some genes change ID over time), so duplicate IDs need to be removed.
non_duplicates_PB=which(duplicated(ids_PB$entrez) == FALSE)
ids_PB=ids_PB[non_duplicates_PB, ]
## Merge the IDs with the results.
# Significant DEGs
res_ids_PB=inner_join(sig_PB, ids_PB, by=c("gene"="ensgene"))
# All genes
all_ids_PB=inner_join(res_table_tb_PB, ids_PB, by=c("gene"="ensgene"))
## Move the Entrez ID and gene symbol columns to after the Ensembl ID column, then show the merged table.
# Significant DEGs
res_ids_PB=res_ids_PB %>%
relocate(entrez, .after=gene) %>%
relocate(symbol, .after=entrez)
res_ids_PB
write.csv(res_ids_PB, "gencode_Output_Files/PB_THC_sig_DEGs_annotated.csv", quote=FALSE, row.names=FALSE)
# All genes
all_ids_PB=all_ids_PB %>%
relocate(entrez, .after=gene) %>%
relocate(symbol, .after=entrez)
all_ids_PB
write.csv(all_ids_PB, "gencode_Output_Files/PB_THC_all_genes_annotated.csv", quote=FALSE, row.names=FALSE)
# PL
## Convert from Ensembl IDs to Entrez IDs (only do this if the dds results has Ensembl IDs).
## Return the IDs for the gene symbols in the DE results
idx_PL=grch38$ensgene %in% rownames(res_table_PL)
ids_PL=grch38[idx_PL, ]
## The Entrez IDs can map to more than one Ensembl ID (some genes change ID over time), so duplicate IDs need to be removed.
non_duplicates_PL=which(duplicated(ids_PL$entrez) == FALSE)
ids_PL=ids_PL[non_duplicates_PL, ]
## Merge the IDs with the results.
# Significant DEGs
res_ids_PL=inner_join(sig_PL, ids_PL, by=c("gene"="ensgene"))
# All genes
all_ids_PL=inner_join(res_table_tb_PL, ids_PL, by=c("gene"="ensgene"))
## Move the Entrez ID and gene symbol columns to after the Ensembl ID column, then show the merged table.
# Significant DEGs
res_ids_PL=res_ids_PL %>%
relocate(entrez, .after=gene) %>%
relocate(symbol, .after=entrez)
res_ids_PL
write.csv(res_ids_PL, "gencode_Output_Files/PL_THC_sig_DEGs_annotated.csv", quote=FALSE, row.names=FALSE)
# All genes
all_ids_PL=all_ids_PL %>%
relocate(entrez, .after=gene) %>%
relocate(symbol, .after=entrez)
all_ids_PL
write.csv(all_ids_PL, "gencode_Output_Files/PL_THC_all_genes_annotated.csv", quote=FALSE, row.names=FALSE)
# PM
## Convert from Ensembl IDs to Entrez IDs (only do this if the dds results has Ensembl IDs).
## Return the IDs for the gene symbols in the DE results
idx_PM=grch38$ensgene %in% rownames(res_table_PM)
ids_PM=grch38[idx_PM, ]
## The Entrez IDs can map to more than one Ensembl ID (some genes change ID over time), so duplicate IDs need to be removed.
non_duplicates_PM=which(duplicated(ids_PM$entrez) == FALSE)
ids_PM=ids_PM[non_duplicates_PM, ]
## Merge the IDs with the results.
# Significant DEGs
res_ids_PM=inner_join(sig_PM, ids_PM, by=c("gene"="ensgene"))
# All genes
all_ids_PM=inner_join(res_table_tb_PM, ids_PM, by=c("gene"="ensgene"))
## Move the Entrez ID and gene symbol columns to after the Ensembl ID column, then show the merged table.
# Significant DEGs
res_ids_PM=res_ids_PM %>%
relocate(entrez, .after=gene) %>%
relocate(symbol, .after=entrez)
res_ids_PM
write.csv(res_ids_PM, "gencode_Output_Files/PM_THC_sig_DEGs_annotated.csv", quote=FALSE, row.names=FALSE)
# All genes
all_ids_PM=all_ids_PM %>%
relocate(entrez, .after=gene) %>%
relocate(symbol, .after=entrez)
all_ids_PM
write.csv(all_ids_PM, "gencode_Output_Files/PM_THC_all_genes_annotated.csv", quote=FALSE, row.names=FALSE)
# PS
## Convert from Ensembl IDs to Entrez IDs (only do this if the dds results has Ensembl IDs).
## Return the IDs for the gene symbols in the DE results
idx_PS=grch38$ensgene %in% rownames(res_table_PS)
ids_PS=grch38[idx_PS, ]
## The Entrez IDs can map to more than one Ensembl ID (some genes change ID over time), so duplicate IDs need to be removed.
non_duplicates_PS=which(duplicated(ids_PS$entrez) == FALSE)
ids_PS=ids_PS[non_duplicates_PS, ]
## Merge the IDs with the results.
# Significant DEGs
res_ids_PS=inner_join(sig_PS, ids_PS, by=c("gene"="ensgene"))
# All genes
all_ids_PS=inner_join(res_table_tb_PS, ids_PS, by=c("gene"="ensgene"))
## Move the Entrez ID and gene symbol columns to after the Ensembl ID column, then show the merged table.
# Significant DEGs
res_ids_PS=res_ids_PS %>%
relocate(entrez, .after=gene) %>%
relocate(symbol, .after=entrez)
res_ids_PS
write.csv(res_ids_PS, "gencode_Output_Files/PS_THC_sig_DEGs_annotated.csv", quote=FALSE, row.names=FALSE)
# All genes
all_ids_PS=all_ids_PS %>%
relocate(entrez, .after=gene) %>%
relocate(symbol, .after=entrez)
all_ids_PS
write.csv(all_ids_PS, "gencode_Output_Files/PS_THC_all_genes_annotated.csv", quote=FALSE, row.names=FALSE)
Annotated Significant DEG Tables
# PL
final_res_ids_PL=as.data.frame(res_ids_PL)
final_res_ids_PL[,4:6]=as.data.frame(round(final_res_ids_PL[,4:6], digits=3))
final_res_ids_PL[,7:8]=format(final_res_ids_PL[,7:8], scientific=TRUE, digits=3)
table_PL=flextable(data=final_res_ids_PL)
table_PL=table_PL %>%
bold(bold=TRUE, part="header") %>%
set_table_properties(layout="autofit") %>%
fit_to_width(10)
table_PL
gene | entrez | symbol | baseMean | log2FoldChange | lfcSE | pvalue | padj | chr | start | end | strand | biotype | description |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
ENSG00000063438 | 57,491 | AHRR | 39.680 | 2.281 | 0.485 | 4.47e-08 | 9.72e-04 | 5 | 321,714 | 438,291 | 1 | protein_coding | aryl hydrocarbon receptor repressor |
ENSG00000158164 | 11,013 | TMSB15A | 21.441 | -1.531 | 0.320 | 1.08e-07 | 1.18e-03 | X | 102,513,682 | 102,516,739 | -1 | protein_coding | thymosin beta 15A |
ENSG00000154165 | 2,838 | GPR15 | 17.183 | 2.215 | 0.547 | 9.34e-07 | 6.78e-03 | 3 | 98,531,978 | 98,534,681 | 1 | protein_coding | G protein-coupled receptor 15 |
ENSG00000159339 | 23,569 | PADI4 | 55.639 | 2.098 | 0.542 | 2.59e-06 | 1.12e-02 | 1 | 17,308,195 | 17,364,004 | 1 | protein_coding | peptidyl arginine deiminase 4 |
ENSG00000173535 | 8,794 | TNFRSF10C | 37.443 | 1.353 | 0.340 | 2.76e-06 | 1.12e-02 | 8 | 23,102,921 | 23,117,445 | 1 | protein_coding | TNF receptor superfamily member 10c |
ENSG00000185669 | 333,929 | SNAI3 | 19.078 | 0.825 | 0.213 | 3.59e-06 | 1.12e-02 | 16 | 88,677,688 | 88,686,507 | -1 | protein_coding | snail family transcriptional repressor 3 |
ENSG00000205038 | 93,035 | PKHD1L1 | 681.956 | 1.212 | 0.310 | 3.60e-06 | 1.12e-02 | 8 | 109,362,461 | 109,537,207 | 1 | protein_coding | PKHD1 like 1 |
ENSG00000158517 | 653,361 | NCF1 | 71.617 | 1.249 | 0.355 | 1.48e-05 | 4.03e-02 | 7 | 74,774,011 | 74,789,315 | 1 | protein_coding | neutrophil cytosolic factor 1 |
ENSG00000165178 | 654,817 | NCF1C | 33.028 | 1.238 | 0.365 | 2.08e-05 | 4.52e-02 | 7 | 75,156,639 | 75,172,044 | -1 | unprocessed_pseudogene | neutrophil cytosolic factor 1C pseudogene |
ENSG00000116701 | 4,688 | NCF2 | 532.308 | 0.762 | 0.226 | 2.93e-05 | 5.11e-02 | 1 | 183,554,461 | 183,590,905 | -1 | protein_coding | neutrophil cytosolic factor 2 |
ENSG00000158683 | 168,507 | PKD1L1 | 200.743 | 1.035 | 0.314 | 3.05e-05 | 5.11e-02 | 7 | 47,740,202 | 47,948,466 | -1 | protein_coding | polycystin 1 like 1, transient receptor potential channel interacting |
ENSG00000185640 | 338,785 | KRT79 | 16.521 | 1.525 | 0.484 | 3.87e-05 | 6.01e-02 | 12 | 52,821,408 | 52,834,311 | -1 | protein_coding | keratin 79 |
ENSG00000134955 | 219,855 | SLC37A2 | 113.411 | 0.694 | 0.215 | 4.73e-05 | 6.86e-02 | 11 | 125,063,302 | 125,090,516 | 1 | protein_coding | solute carrier family 37 member 2 |
ENSG00000102575 | 54 | ACP5 | 265.403 | 0.674 | 0.224 | 9.04e-05 | 9.62e-02 | 19 | 11,574,653 | 11,579,993 | -1 | protein_coding | acid phosphatase 5, tartrate resistant |
ENSG00000117115 | 11,240 | PADI2 | 100.079 | 0.916 | 0.310 | 1.03e-04 | 9.62e-02 | 1 | 17,066,761 | 17,119,451 | -1 | protein_coding | peptidyl arginine deiminase 2 |
ENSG00000140465 | 1,543 | CYP1A1 | 308.918 | 2.439 | 0.797 | 7.14e-05 | 9.62e-02 | 15 | 74,719,542 | 74,725,536 | -1 | protein_coding | cytochrome P450 family 1 subfamily A member 1 |
ENSG00000159399 | 3,099 | HK2 | 499.815 | 0.722 | 0.238 | 8.89e-05 | 9.62e-02 | 2 | 74,834,127 | 74,893,359 | 1 | protein_coding | hexokinase 2 |
ENSG00000167680 | 10,501 | SEMA6B | 212.123 | 0.938 | 0.312 | 8.34e-05 | 9.62e-02 | 19 | 4,542,593 | 4,581,776 | -1 | protein_coding | semaphorin 6B |
# PM
final_res_ids_PM=as.data.frame(res_ids_PM)
final_res_ids_PM[,4:6]=as.data.frame(round(final_res_ids_PM[,4:6], digits=3))
final_res_ids_PM[,7:8]=format(final_res_ids_PM[,7:8], scientific=TRUE, digits=3)
table_PM=flextable(data=final_res_ids_PM)
table_PM=table_PM %>%
bold(bold=TRUE, part="header") %>%
set_table_properties(layout="autofit") %>%
fit_to_width(9.19)
table_PM
gene | entrez | symbol | baseMean | log2FoldChange | lfcSE | pvalue | padj | chr | start | end | strand | biotype | description |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
ENSG00000140465 | 1,543 | CYP1A1 | 115.075 | 2.153 | 0.644 | 2.50e-20 | 1.40e-15 | 15 | 74,719,542 | 74,725,536 | -1 | protein_coding | cytochrome P450 family 1 subfamily A member 1 |
ENSG00000124818 | 221,391 | OPN5 | 12.025 | -1.740 | 0.329 | 4.50e-09 | 1.26e-04 | 6 | 47,781,982 | 47,832,780 | 1 | protein_coding | opsin 5 |
ENSG00000074410 | 771 | CA12 | 10.452 | 1.967 | 0.506 | 2.97e-06 | 4.14e-02 | 15 | 63,321,378 | 63,381,846 | -1 | protein_coding | carbonic anhydrase 12 |
ENSG00000163435 | 1,999 | ELF3 | 7.176 | 2.085 | 0.533 | 2.58e-06 | 4.14e-02 | 1 | 202,007,945 | 202,017,183 | 1 | protein_coding | E74 like ETS transcription factor 3 |
sect_properties=prop_section(
page_size=page_size(
orient="landscape",
width=8.5, height=11),
type="continuous",
page_margins=page_mar())
save_as_docx(table_PL, table_PM,
path="postmortem_tables.docx", pr_section=sect_properties)
## Summarize data: omit rows containing NAs, group data by "tissue", then calculate mean and standard deviation for "concentration," "A260_280," "A260_230," and "RIN" and generate a table with the output
tissue_data_na=na.omit(metadata_all)
tissue_data.summary=tissue_data_na %>%
group_by(tissue) %>%
summarise(
mean_conc=mean(concentration),
sd_conc=sd(concentration),
mean_RIN=mean(RIN),
sd_RIN=sd(RIN),
mean_A260_280=mean(A260_280),
sd_A260_280=sd(A260_280),
mean_A260_230=mean(A260_230),
sd_A260_230=sd(A260_230)
)
tissue_data.summary
tissue_data_na %>%
summarise(mean(RIN), sd(RIN))
RNA_stats_PM=as.data.frame(tissue_data.summary)
RNA_stats_PM[,2:9]=as.data.frame(round(RNA_stats_PM[,2:9], digits=3))
table_RNA_stats=flextable(data=RNA_stats_PM)
table_RNA_stats=table_RNA_stats %>%
bold(bold=TRUE, part="header") %>%
set_table_properties(layout="autofit") %>%
fit_to_width(6)
table_RNA_stats
tissue | mean_conc | sd_conc | mean_RIN | sd_RIN | mean_A260_280 | sd_A260_280 | mean_A260_230 | sd_A260_230 |
---|---|---|---|---|---|---|---|---|
blood | 242.279 | 346.551 | 3.094 | 1.236 | 1.992 | 0.087 | 2.100 | 0.218 |
brain | 506.787 | 457.304 | 4.196 | 0.978 | 2.070 | 0.027 | 2.156 | 0.143 |
lung | 164.391 | 154.572 | 3.927 | 1.417 | 2.086 | 0.025 | 2.132 | 0.119 |
muscle | 319.542 | 178.790 | 5.027 | 1.467 | 2.079 | 0.018 | 2.134 | 0.098 |
save_as_docx(table_RNA_stats, path="postmortem_RNA_tables.docx", pr_section=sect_properties)
# Concentration
RNA1=ggplot()+
# create boxplot layer
geom_boxplot(data=tissue_data_na, aes(x=tissue, y=concentration, fill=tissue),
width=0.75, lwd=0.75, show.legend=FALSE)+
# set the x-axis scale to the limits and labels defined below
scale_x_discrete(limits=c("brain", "lung", "muscle", "blood"),
labels=c("Brain", "Lung", "Muscle", "Blood"))+
# set fill color scale
scale_fill_manual(values=alpha(c("black", "black", "black", "black", "black", "black"), 0.25), guide="none",
limits=c("brain", "lung", "muscle", "blood"),
labels=c("Brain", "Lung", "Muscle", "Blood"))+
# change plot scaling ratio
coord_fixed(ratio=1/305)+
# define theme variables
theme_minimal()+
theme(axis.text.y=element_text(color="black", face="bold", size=rel(1)),
axis.text.x=element_text(angle=30, hjust=1, color="black", face="bold", size=rel(1)),
# axis.title.x=element_text(color="black", face="bold", size=rel(0.80)),
axis.title.x=element_blank(),
axis.title.y=element_text(color="black", face="bold", size=rel(1)),
panel.border=element_rect(color="black", linewidth=1, fill=NA),
panel.grid.minor.y=element_blank(),
strip.background=element_rect(color="black", linewidth=1),
strip.text=element_text(color="black", face="bold", size=rel(0.8)),
axis.line=element_blank(),
axis.ticks=element_line(linewidth=1),
legend.position="none",
legend.key.size=unit(rel(0.50), "cm"),
legend.title=element_text(color="black", face="bold", size=rel(0.60)),
legend.text=element_text(color="black", face="bold", size=rel(0.60)),
plot.title=element_text(color="black", face="bold", size=rel(1)),
plot.margin=unit(c(0.10, 0.10, 0.10, 0.10), "cm"))+
# x-axis label
xlab("Tissue")+
# y-axis label
ylab("Concentration (ng/\u00B5L)")
# A260_280
RNA2=ggplot()+
# create boxplot layer
geom_boxplot(data=tissue_data_na, aes(x=tissue, y=A260_280, fill=tissue),
width=0.75, lwd=0.75, show.legend=FALSE)+
# set the x-axis scale to the limits and labels defined below
scale_x_discrete(limits=c("brain", "lung", "muscle", "blood"),
labels=c("Brain", "Lung", "Muscle", "Blood"))+
# set fill color scale
scale_fill_manual(values=alpha(c("black", "black", "black", "black", "black", "black"), 0.25), guide="none",
limits=c("brain", "lung", "muscle", "blood"),
labels=c("Brain", "Lung", "Muscle", "Blood"))+
# change plot scaling ratio
coord_fixed(ratio=20.5/1)+
# define theme variables
theme_minimal()+
theme(axis.text.y=element_text(color="black", face="bold", size=rel(1)),
axis.text.x=element_text(angle=30, hjust=1, color="black", face="bold", size=rel(1)),
# axis.title.x=element_text(color="black", face="bold", size=rel(0.80)),
axis.title.x=element_blank(),
axis.title.y=element_text(color="black", face="bold", size=rel(1)),
panel.border=element_rect(color="black", linewidth=1, fill=NA),
panel.grid.minor.y=element_blank(),
strip.background=element_rect(color="black", linewidth=1),
strip.text=element_text(color="black", face="bold", size=rel(0.8)),
axis.line=element_blank(),
axis.ticks=element_line(linewidth=1),
legend.position="none",
legend.key.size=unit(rel(0.50), "cm"),
legend.title=element_text(color="black", face="bold", size=rel(0.60)),
legend.text=element_text(color="black", face="bold", size=rel(0.60)),
plot.title=element_text(color="black", face="bold", size=rel(1)),
plot.margin=unit(c(0.10, 0.10, 0.10, 0.10), "cm"))+
# x-axis label
xlab("Tissue")+
# y-axis label
ylab("A260/280 Ratio")
# A260_230
RNA3=ggplot()+
# create boxplot layer
geom_boxplot(data=tissue_data_na, aes(x=tissue, y=A260_230, fill=tissue),
width=0.75, lwd=0.75, show.legend=FALSE)+
# set the x-axis scale to the limits and labels defined below
scale_x_discrete(limits=c("brain", "lung", "muscle", "blood"),
labels=c("Brain", "Lung", "Muscle", "Blood"))+
# set fill color scale
scale_fill_manual(values=alpha(c("black", "black", "black", "black", "black", "black"), 0.25), guide="none",
limits=c("brain", "lung", "muscle", "blood"),
labels=c("Brain", "Lung", "Muscle", "Blood"))+
# change plot scaling ratio
coord_fixed(ratio=8.2/1)+
# define theme variables
theme_minimal()+
theme(axis.text.y=element_text(color="black", face="bold", size=rel(1)),
axis.text.x=element_text(angle=30, hjust=1, color="black", face="bold", size=rel(1)),
# axis.title.x=element_text(color="black", face="bold", size=rel(0.80)),
axis.title.x=element_blank(),
axis.title.y=element_text(color="black", face="bold", size=rel(1)),
panel.border=element_rect(color="black", linewidth=1, fill=NA),
panel.grid.minor.y=element_blank(),
strip.background=element_rect(color="black", linewidth=1),
strip.text=element_text(color="black", face="bold", size=rel(0.8)),
axis.line=element_blank(),
axis.ticks=element_line(linewidth=1),
legend.position="none",
legend.key.size=unit(rel(0.50), "cm"),
legend.title=element_text(color="black", face="bold", size=rel(0.60)),
legend.text=element_text(color="black", face="bold", size=rel(0.60)),
plot.title=element_text(color="black", face="bold", size=rel(1)),
plot.margin=unit(c(0.10, 0.10, 0.10, 0.10), "cm"))+
# x-axis label
xlab("Tissue")+
# y-axis label
ylab("A260/230 Ratio")
# RIN
RNA4=ggplot()+
# create boxplot layer
geom_boxplot(data=tissue_data_na, aes(x=tissue, y=RIN, fill=tissue),
width=0.75, lwd=0.75, show.legend=FALSE)+
# set the x-axis scale to the limits and labels defined below
scale_x_discrete(limits=c("brain", "lung", "muscle", "blood"),
labels=c("Brain", "Lung", "Muscle", "Blood"))+
# set fill color scale
scale_fill_manual(values=alpha(c("black", "black", "black", "black", "black", "black"), 0.25), guide="none",
limits=c("brain", "lung", "muscle", "blood"),
labels=c("Brain", "Lung", "Muscle", "Blood"))+
# change plot scaling ratio
coord_fixed(ratio=1.25/1)+
# define theme variables
theme_minimal()+
theme(axis.text.y=element_text(color="black", face="bold", size=rel(1)),
axis.text.x=element_text(angle=30, hjust=1, color="black", face="bold", size=rel(1)),
# axis.title.x=element_text(color="black", face="bold", size=rel(0.80)),
axis.title.x=element_blank(),
axis.title.y=element_text(color="black", face="bold", size=rel(1)),
panel.border=element_rect(color="black", linewidth=1, fill=NA),
panel.grid.minor.y=element_blank(),
strip.background=element_rect(color="black", linewidth=1),
strip.text=element_text(color="black", face="bold", size=rel(0.8)),
axis.line=element_blank(),
axis.ticks=element_line(linewidth=1),
legend.position="none",
legend.key.size=unit(rel(0.50), "cm"),
legend.title=element_text(color="black", face="bold", size=rel(0.60)),
legend.text=element_text(color="black", face="bold", size=rel(0.60)),
plot.title=element_text(color="black", face="bold", size=rel(1)),
plot.margin=unit(c(0.10, 0.10, 0.10, 0.10), "cm"))+
# x-axis label
xlab("Tissue")+
# y-axis label
ylab("RNA Integrity Number")
RNA_patchwork=RNA1 + RNA4 + RNA2 + RNA3 + plot_layout(ncol=2)
RNA_patchwork + plot_annotation(tag_levels="A")
ggsave("PM_RNA_plots.png", width=4.845, height=8, unit="in", dpi=320)