6_Variance_Partition_Analysis
Yang & Steve
07/10/2020
Last updated: 2020-10-28
Checks: 7 0
Knit directory: 20190717_Lardelli_RNASeq_Larvae/
This reproducible R Markdown analysis was created with workflowr (version 1.6.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20200227) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: 1_DE-gene-analysis_cache/
Ignored: 1_DE-gene-analysis_files/
Ignored: analysis/.DS_Store
Ignored: analysis/.Rhistory
Ignored: analysis/.Rproj.user/
Ignored: analysis/5_WGCNA_cache/
Ignored: data/.DS_Store
Ignored: data/0_rawData/.DS_Store
Ignored: data/1_trimmedData/.DS_Store
Ignored: data/2_alignedData/.DS_Store
Ignored: files/
Ignored: keggdiagram/.DS_Store
Ignored: output/.DS_Store
Unstaged changes:
Modified: analysis/index.Rmd
Staged changes:
Modified: analysis/index.Rmd
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 417413a | yangdongau | 2020-10-28 | Add in WGCNA and variance partitioning analysis |
| html | 417413a | yangdongau | 2020-10-28 | Add in WGCNA and variance partitioning analysis |
Setup
library(limma)
library(edgeR)
library(tidyverse)
library(magrittr)
library(pander)
library(ggrepel)
library(scales)
library(variancePartition)
library(lme4)
theme_set(theme_bw())
panderOptions("big.mark", ",")
panderOptions("table.split.table", Inf)
panderOptions("table.style", "rmarkdown")
if (interactive()) setwd(here::here("analysis"))Data load
dgeList <- read_rds(here::here("data","dgeList.rds"))
entrezGenes <- dgeList$genes %>%
dplyr::filter(!is.na(entrez_gene)) %>%
unnest(entrez_gene) %>%
dplyr::rename(entrez_gene = entrez_gene)
topTable <- file.path(here::here("output", "topTable.csv")) %>%
read_csv()
topTableDE <- file.path(here::here("output", "DEgenes.csv")) %>%
read_csv()Application of variance partition analysis
Pair and genotype are used as random factors.
geneExpr <- dgeList$counts
form <- ~ (1|Genotype) + (1|pair)
varPar <- fitExtractVarPartModel(geneExpr, form, dgeList$samples)Dividing work into 100 chunks...
Total: 210 svp <- sortCols(varPar)
plotPercentBars(vp[1:10,] )
| Version | Author | Date |
|---|---|---|
| 417413a | yangdongau | 2020-10-28 |
plotVarPart(vp)
| Version | Author | Date |
|---|---|---|
| 417413a | yangdongau | 2020-10-28 |
head(varPar) Genotype pair Residuals
ENSDARG00000000001 3.313457e-01 1.524640e-10 0.6686543
ENSDARG00000000002 0.000000e+00 0.000000e+00 1.0000000
ENSDARG00000000018 1.230204e-09 2.919213e-01 0.7080787
ENSDARG00000000019 0.000000e+00 0.000000e+00 1.0000000
ENSDARG00000000068 2.184783e-02 4.562347e-01 0.5219175
ENSDARG00000000069 0.000000e+00 0.000000e+00 1.0000000Plot expression stratified by variables
# Get the gene with the highest variation between genotypes
i_genotype <- which.max(varPar$Genotype)
GE_genotype <- data.frame(Expression = geneExpr[i_genotype,], Genotype = dgeList$samples$Genotype)
# Plot expression stratified by Genotype
label_genotype <- paste("Genotype:", format(varPar$Genotype[i_genotype]*100,
digits=3), "%")
plotStratify(Expression ~ Genotype, GE_genotype, text=label_genotype, main=rownames(geneExpr)[i_genotype])
| Version | Author | Date |
|---|---|---|
| 417413a | yangdongau | 2020-10-28 |
# Get the gene with the highest variation among pairs
i_pair <- which.max(varPar$pair)
GE_pair <- data.frame(Expression = geneExpr[i_pair,], pair = dgeList$samples$pair)
# Plot expression stratified by pair
label_pair <- paste("Pair:", format(varPar$pair[i_pair]*100,
digits=3), "%")
plotStratify(Expression ~ pair, GE_pair, text=label_pair, main=rownames(geneExpr)[i_pair])
| Version | Author | Date |
|---|---|---|
| 417413a | yangdongau | 2020-10-28 |
# no colour
plotStratify(Expression ~ pair, GE_pair, colorBy = NULL, text=label_pair, main=rownames(geneExpr)[i_pair])
| Version | Author | Date |
|---|---|---|
| 417413a | yangdongau | 2020-10-28 |
i_genotype_list <- varPar[order(-varPar$Genotype),]
i_genotype_list <- cbind(rownames(i_genotype_list), i_genotype_list)
rownames(i_genotype_list) <- NULL
colnames(i_genotype_list) <- c("gene_id","Genotype","pair","Residuals")
gene_names <- dgeList$genes %>%
dplyr::select(external_gene_name)
gene_names <- cbind(rownames(gene_names),gene_names)
rownames(gene_names) <- NULL
colnames(gene_names) <- c("gene_id", "gene_name")
i_genotype_list <- i_genotype_list %>%
left_join(gene_names)
i_genotype_list[1:10,] gene_id Genotype pair Residuals gene_name
1 ENSDARG00000070140 0.7987162 7.650881e-02 0.12477504 RETSAT
2 ENSDARG00000054510 0.7922638 5.726178e-10 0.20773616 adrb2b
3 ENSDARG00000076870 0.7809935 2.414483e-02 0.19486172 piezo1
4 ENSDARG00000101331 0.7795316 0.000000e+00 0.22046840 tekt1
5 ENSDARG00000037613 0.7475010 1.683271e-01 0.08417193 lgals8b
6 ENSDARG00000102658 0.7451735 2.510753e-09 0.25482651 zgc:174624
7 ENSDARG00000076805 0.7435091 1.836092e-02 0.23812994 leng8
8 ENSDARG00000099511 0.7434086 6.280672e-08 0.25659136 CABZ01034698.2
9 ENSDARG00000032430 0.7407945 8.080070e-02 0.17840480 ppp2r1bb
10 ENSDARG00000030110 0.7393170 5.459288e-02 0.20609012 myod1Add Genotype column to DE gene list
i_genotype_list_genotype <- i_genotype_list %>%
dplyr::select(gene_id, Genotype)
colnames(i_genotype_list_genotype) <- c("ensembl_gene_id","Genotype")
topTableDE_genotype <- topTableDE %>%
left_join(i_genotype_list_genotype)
# Save results
write.csv(topTableDE_genotype, here::here("output", "DEgenes_with_genotype.csv"))
devtools::session_info()─ Session info ──────────────────────────────────────────────────────────
setting value
version R version 3.6.0 (2019-04-26)
os macOS Mojave 10.14.6
system x86_64, darwin15.6.0
ui X11
language (EN)
collate en_AU.UTF-8
ctype en_AU.UTF-8
tz Australia/Adelaide
date 2020-10-28
─ Packages ──────────────────────────────────────────────────────────────
package * version date lib source
assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.0)
backports 1.1.4 2019-04-10 [1] CRAN (R 3.6.0)
Biobase * 2.44.0 2019-05-02 [1] Bioconductor
BiocGenerics * 0.30.0 2019-05-02 [1] Bioconductor
BiocParallel 1.18.1 2019-08-06 [1] Bioconductor
bitops 1.0-6 2013-08-17 [1] CRAN (R 3.6.0)
boot 1.3-23 2019-07-05 [1] CRAN (R 3.6.0)
broom 0.5.2 2019-04-07 [1] CRAN (R 3.6.0)
callr 3.3.1 2019-07-18 [1] CRAN (R 3.6.0)
caTools 1.17.1.2 2019-03-06 [1] CRAN (R 3.6.0)
cellranger 1.1.0 2016-07-27 [1] CRAN (R 3.6.0)
cli 1.1.0 2019-03-19 [1] CRAN (R 3.6.0)
codetools 0.2-16 2018-12-24 [1] CRAN (R 3.6.0)
colorRamps 2.3 2012-10-29 [1] CRAN (R 3.6.0)
colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.6.0)
crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.0)
desc 1.2.0 2018-05-01 [1] CRAN (R 3.6.0)
devtools 2.2.2 2020-02-17 [1] CRAN (R 3.6.0)
digest 0.6.20 2019-07-04 [1] CRAN (R 3.6.0)
doParallel 1.0.15 2019-08-02 [1] CRAN (R 3.6.0)
dplyr * 0.8.3 2019-07-04 [1] CRAN (R 3.6.0)
edgeR * 3.26.7 2019-08-13 [1] Bioconductor
ellipsis 0.3.0 2019-09-20 [1] CRAN (R 3.6.0)
evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.0)
forcats * 0.4.0 2019-02-17 [1] CRAN (R 3.6.0)
foreach * 1.4.7 2019-07-27 [1] CRAN (R 3.6.0)
fs 1.3.1 2019-05-06 [1] CRAN (R 3.6.0)
gdata 2.18.0 2017-06-06 [1] CRAN (R 3.6.0)
generics 0.0.2 2018-11-29 [1] CRAN (R 3.6.0)
ggplot2 * 3.2.1 2019-08-10 [1] CRAN (R 3.6.0)
ggrepel * 0.8.1 2019-05-07 [1] CRAN (R 3.6.0)
git2r 0.26.1 2019-06-29 [1] CRAN (R 3.6.0)
glue 1.3.1 2019-03-12 [1] CRAN (R 3.6.0)
gplots 3.0.1.1 2019-01-27 [1] CRAN (R 3.6.0)
gtable 0.3.0 2019-03-25 [1] CRAN (R 3.6.0)
gtools 3.8.1 2018-06-26 [1] CRAN (R 3.6.0)
haven 2.1.1 2019-07-04 [1] CRAN (R 3.6.0)
here 0.1 2017-05-28 [1] CRAN (R 3.6.0)
hms 0.5.1 2019-08-23 [1] CRAN (R 3.6.0)
htmltools 0.3.6 2017-04-28 [1] CRAN (R 3.6.0)
httpuv 1.5.1 2019-04-05 [1] CRAN (R 3.6.0)
httr 1.4.1 2019-08-05 [1] CRAN (R 3.6.0)
iterators 1.0.12 2019-07-26 [1] CRAN (R 3.6.0)
jsonlite 1.6 2018-12-07 [1] CRAN (R 3.6.0)
KernSmooth 2.23-15 2015-06-29 [1] CRAN (R 3.6.0)
knitr 1.24 2019-08-08 [1] CRAN (R 3.6.0)
labeling 0.3 2014-08-23 [1] CRAN (R 3.6.0)
later 0.8.0 2019-02-11 [1] CRAN (R 3.6.0)
lattice 0.20-38 2018-11-04 [1] CRAN (R 3.6.0)
lazyeval 0.2.2 2019-03-15 [1] CRAN (R 3.6.0)
limma * 3.40.6 2019-07-26 [1] Bioconductor
lme4 * 1.1-21 2019-03-05 [1] CRAN (R 3.6.0)
locfit 1.5-9.1 2013-04-20 [1] CRAN (R 3.6.0)
lubridate 1.7.4 2018-04-11 [1] CRAN (R 3.6.0)
magrittr * 1.5 2014-11-22 [1] CRAN (R 3.6.0)
MASS 7.3-51.4 2019-03-31 [1] CRAN (R 3.6.0)
Matrix * 1.2-17 2019-03-22 [1] CRAN (R 3.6.0)
memoise 1.1.0 2017-04-21 [1] CRAN (R 3.6.0)
minqa 1.2.4 2014-10-09 [1] CRAN (R 3.6.0)
modelr 0.1.5 2019-08-08 [1] CRAN (R 3.6.0)
munsell 0.5.0 2018-06-12 [1] CRAN (R 3.6.0)
nlme 3.1-141 2019-08-01 [1] CRAN (R 3.6.0)
nloptr 1.2.1 2018-10-03 [1] CRAN (R 3.6.0)
pander * 0.6.3 2018-11-06 [1] CRAN (R 3.6.0)
pbkrtest 0.4-7 2017-03-15 [1] CRAN (R 3.6.0)
pillar 1.4.2 2019-06-29 [1] CRAN (R 3.6.0)
pkgbuild 1.0.6 2019-10-09 [1] CRAN (R 3.6.0)
pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.6.0)
pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.6.0)
plyr 1.8.4 2016-06-08 [1] CRAN (R 3.6.0)
prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.6.0)
processx 3.4.1 2019-07-18 [1] CRAN (R 3.6.0)
progress 1.2.2 2019-05-16 [1] CRAN (R 3.6.0)
promises 1.0.1 2018-04-13 [1] CRAN (R 3.6.0)
ps 1.3.0 2018-12-21 [1] CRAN (R 3.6.0)
purrr * 0.3.3 2019-10-18 [1] CRAN (R 3.6.0)
R6 2.4.0 2019-02-14 [1] CRAN (R 3.6.0)
Rcpp 1.0.2 2019-07-25 [1] CRAN (R 3.6.0)
readr * 1.3.1 2018-12-21 [1] CRAN (R 3.6.0)
readxl 1.3.1 2019-03-13 [1] CRAN (R 3.6.0)
remotes 2.1.1 2020-02-15 [1] CRAN (R 3.6.0)
reshape2 1.4.3 2017-12-11 [1] CRAN (R 3.6.0)
rlang 0.4.4 2020-01-28 [1] CRAN (R 3.6.0)
rmarkdown 1.15 2019-08-21 [1] CRAN (R 3.6.0)
rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.6.0)
rstudioapi 0.10 2019-03-19 [1] CRAN (R 3.6.0)
rvest 0.3.4 2019-05-15 [1] CRAN (R 3.6.0)
scales * 1.0.0 2018-08-09 [1] CRAN (R 3.6.0)
sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.0)
stringi 1.4.3 2019-03-12 [1] CRAN (R 3.6.0)
stringr * 1.4.0 2019-02-10 [1] CRAN (R 3.6.0)
testthat 2.3.1 2019-12-01 [1] CRAN (R 3.6.0)
tibble * 2.1.3 2019-06-06 [1] CRAN (R 3.6.0)
tidyr * 0.8.3 2019-03-01 [1] CRAN (R 3.6.0)
tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.6.0)
tidyverse * 1.2.1 2017-11-14 [1] CRAN (R 3.6.0)
usethis 1.5.1 2019-07-04 [1] CRAN (R 3.6.0)
variancePartition * 1.14.1 2019-10-01 [1] Bioconductor
vctrs 0.2.0 2019-07-05 [1] CRAN (R 3.6.0)
whisker 0.4 2019-08-28 [1] CRAN (R 3.6.0)
withr 2.1.2 2018-03-15 [1] CRAN (R 3.6.0)
workflowr 1.6.0 2019-12-19 [1] CRAN (R 3.6.0)
xfun 0.9 2019-08-21 [1] CRAN (R 3.6.0)
xml2 1.2.2 2019-08-09 [1] CRAN (R 3.6.0)
yaml 2.2.0 2018-07-25 [1] CRAN (R 3.6.0)
zeallot 0.1.0 2018-01-28 [1] CRAN (R 3.6.0)
[1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library