Limma-voom paired analysis of ALL samples
Last updated: 2022-04-27
Checks: 7 0
Knit directory: amnio-cell-free-RNA/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). 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(20200224) 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 results in this page were generated with repository version 43f8f24. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
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: .Rhistory
Ignored: .Rproj.user/
Ignored: .bpipe/
Ignored: analysis/obsolete_analysis/
Ignored: code/.bpipe/
Ignored: code/.rnaseq-test.groovy.swp
Ignored: code/obsolete_analysis/
Ignored: data/.bpipe/
Ignored: data/190717_A00692_0021_AHLLHFDSXX/
Ignored: data/190729_A00692_0022_AHLLHFDSXX/
Ignored: data/190802_A00692_0023_AHLLHFDSXX/
Ignored: data/200612_A00692_0107_AHN3YCDMXX.tar
Ignored: data/200612_A00692_0107_AHN3YCDMXX/
Ignored: data/200626_A00692_0111_AHNJH7DMXX.tar
Ignored: data/200626_A00692_0111_AHNJH7DMXX/
Ignored: data/CMV-AF-database-final-included-samples.csv
Ignored: data/GONE4.10.13.txt
Ignored: data/HK_exons.csv
Ignored: data/HK_exons.txt
Ignored: data/IPA molecule summary.xls
Ignored: data/IPA-molecule-summary.csv
Ignored: data/deduped_rRNA_coverage.txt
Ignored: data/gene-transcriptome-analysis/
Ignored: data/hg38_rRNA.bed
Ignored: data/hg38_rRNA.saf
Ignored: data/ignore-overlap-mapping/
Ignored: data/ignore/
Ignored: data/rds/
Ignored: data/salmon-pilot-analysis/
Ignored: data/star-genome-analysis/
Ignored: output/c2Ens.RData
Ignored: output/c5Ens.RData
Ignored: output/hEns.RData
Ignored: output/keggEns.RData
Ignored: output/obsolete_output/
Unstaged changes:
Modified: .gitignore
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 repository in which changes were made to the R Markdown (analysis/STAR-FC-all.Rmd) and HTML (docs/STAR-FC-all.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 43f8f24 | Jovana Maksimovic | 2022-04-27 | wflow_publish(c("analysis/STAR-FC-all.Rmd", "analysis/STAR-FC-exclude-US-ab.Rmd", |
| html | 3b51ada | Jovana Maksimovic | 2021-12-03 | Build site. |
| Rmd | 85a6971 | Jovana Maksimovic | 2021-12-03 | wflow_publish(c("analysis/STAR-FC-all.Rmd", "analysis/STAR-FC-exclude-US-ab.Rmd", |
| html | 41ad3da | Jovana Maksimovic | 2021-11-24 | Build site. |
| Rmd | d2811a1 | Jovana Maksimovic | 2021-11-24 | wflow_publish(files = paste0("analysis/", list.files("analysis/", |
| html | ccd6c80 | Jovana Maksimovic | 2021-10-29 | Build site. |
| Rmd | 02eba43 | Jovana Maksimovic | 2021-10-29 | wflow_publish(c("analysis/index.Rmd", "analysis/STAR-FC-all.Rmd", |
| html | a3deccf | Jovana Maksimovic | 2021-09-27 | Build site. |
| Rmd | af2af97 | Jovana Maksimovic | 2021-09-27 | wflow_publish(c("analysis/index.Rmd", "analysis/STAR-FC-all.Rmd", |
| Rmd | 2b6479a | Jovana Maksimovic | 2021-07-26 | Move/remove olds files. |
| Rmd | 76d6992 | Jovana Maksimovic | 2020-12-04 | Explore rRNA counts |
| html | d96660b | Jovana Maksimovic | 2020-11-09 | Build site. |
| Rmd | d32e63b | Jovana Maksimovic | 2020-11-09 | wflow_publish(c("analysis/STAR-FC-all.Rmd", "analysis/STAR-FC-exclude-US-ab.Rmd", |
| html | 10dcedf | Jovana Maksimovic | 2020-11-06 | Build site. |
| Rmd | 787e0a1 | Jovana Maksimovic | 2020-11-06 | wflow_publish(c("analysis/STAR-FC-all.Rmd", "analysis/STAR-FC-exclude-US-ab.Rmd", |
| html | a91102d | Jovana Maksimovic | 2020-10-13 | Build site. |
| Rmd | 2ee5593 | Jovana Maksimovic | 2020-10-13 | wflow_publish(c("analysis/index.Rmd", "analysis/STAR-FC-all.Rmd", |
library(here)
library(tidyverse)
library(EnsDb.Hsapiens.v86)
library(readr)
library(limma)
library(edgeR)
library(NMF)
library(patchwork)
library(EGSEA)
source(here("code/output.R"))The data showed some adapter contamination and sequence duplication issues. Adapters were removed using Trimmomatic and both paired and unpaired reads were retained. Only paired reads were initially mapped with Star in conjunction with GRCh38 and gencode_v34 to detect all junctions, across all samples. Paired and unpaired reads were then mapped to GRCh38 separately using Star. Duplicates were removed from paired and unpaired mapped data using Picard MarkDuplicates. Reads were then counted across features from gencode_v34 using featureCounts.
Data import
Set up DGElist object for downstream analysis. Sum paired and unpaired counts prior to downstream analysis.
rawPE <- read_delim(here("data/star-genome-analysis/counts-pe/counts.txt"), delim = "\t", skip = 1)
rawSE <- read_delim(here("data/star-genome-analysis/counts-se/counts.txt"), delim = "\t", skip = 1)
samps <- strsplit2(colnames(rawPE)[c(7:ncol(rawPE))], "_")[,5]
batch <- factor(strsplit2(colnames(rawPE)[c(7:ncol(rawPE))],
"_")[,1], labels = 1:2)
batch <- tibble(batch = batch, id = samps)
colnames(rawPE)[7:ncol(rawPE)] <- samps
colnames(rawSE)[7:ncol(rawSE)] <- samps
counts <- rawPE[, 7:ncol(rawPE)] + rawSE[, 7:ncol(rawSE)]
dge <- DGEList(counts = counts,
genes = rawPE[,c(1,6)])
dgeAn object of class "DGEList"
$counts
CMV30 CMV31 CMV8 CMV9 CMV26 CMV27 CMV14 CMV15 CMV20 CMV21 CMV1 CMV2 CMV3 CMV4
1 0 0 0 0 2 2 0 1 0 1 1 0 1 0
2 58 95 58 59 113 101 60 48 79 71 54 63 39 46
3 1 0 0 0 0 0 0 0 0 0 0 0 0 0
4 0 0 0 0 0 0 0 0 0 0 1 0 1 0
5 0 0 0 0 0 0 0 0 0 0 0 0 0 0
CMV10 CMV11 CMV18 CMV19 CMV35 Corriel NTC-2 CMV51 CMV52 CMV53 CMV54 CMV56
1 0 0 0 0 1 1 0 0 0 0 0 0
2 62 35 51 45 59 84 0 63 28 49 46 37
3 0 0 0 0 0 1 0 0 0 0 0 0
4 0 0 0 0 0 0 0 0 0 0 0 0
5 0 0 0 0 0 0 0 0 0 0 0 0
CMV57 CMV58 CMV60 CMV61
1 0 2 1 2
2 59 82 44 36
3 0 0 0 0
4 0 0 1 0
5 0 0 0 0
60664 more rows ...
$samples
group lib.size norm.factors
CMV30 1 4673630 1
CMV31 1 5232010 1
CMV8 1 3594801 1
CMV9 1 3425478 1
CMV26 1 4892776 1
25 more rows ...
$genes
Geneid Length
1 ENSG00000223972.5 1735
2 ENSG00000227232.5 1351
3 ENSG00000278267.1 68
4 ENSG00000243485.5 1021
5 ENSG00000284332.1 138
60664 more rows ...Load sample information and file names.
samps1 <- read_csv(here("data/CMV-AF-database-corrected-oct-2020.csv"))
samps2 <- read_csv(here("data/samples.csv"))
samps1 %>% full_join(samps2, by = c("sequencing_ID" = "SampleId")) %>%
mutate(pair = ifelse(!is.na(matched_pair), matched_pair,
ifelse(!is.na(MatchedPair), MatchedPair, NA)),
CMV_status = ifelse(!is.na(CMV_status), CMV_status,
ifelse(!is.na(TestResult), TestResult, NA)),
Sex = toupper(Sex),
Indication = tolower(Indication)) %>%
dplyr::rename(sex = Sex,
id = sequencing_ID,
indication = Indication,
GA_at_amnio = `GA_at_amnio-completed_weeks`) -> samps
read_csv(file = here("data/metadata.csv")) %>%
inner_join(read_csv(file = here("data/joindata.csv")),
by = c("Record.ID" = "UR")) %>%
right_join(samps, by = c("ID post-extraction" = "id")) %>%
na_if("NA") %>%
mutate(sex = ifelse(!is.na(sex), sex,
ifelse(!is.na(Sex), toupper(Sex), NA)),
GA_at_amnio = ifelse(!is.na(GA_at_amnio), GA_at_amnio,
ifelse(!is.na(GA.at.amnio), GA.at.amnio, NA))) %>%
dplyr::rename(id = `ID post-extraction`) %>%
dplyr::select(id,
CMV_status,
pair,
sex,
GA_at_amnio,
indication) %>%
left_join(batch) %>%
dplyr::filter(id %in% colnames(dge)) %>%
drop_na() -> targets
m <- match(colnames(dge), targets$id)
targets <- targets[m[!is.na(m)], ]
targets# A tibble: 26 x 7
id CMV_status pair sex GA_at_amnio indication batch
<chr> <chr> <chr> <chr> <chr> <chr> <fct>
1 CMV30 pos L1 F 21 no_us_ab 1
2 CMV31 neg L1 F 21 no_us_ab 1
3 CMV8 neg L2 F 23 no_us_ab 1
4 CMV9 pos L2 F 23 no_us_ab 1
5 CMV26 pos L3 F 22 no_us_ab 1
6 CMV14 neg L4 F 21 no_us_ab 1
7 CMV15 pos L4 F 22 no_us_ab 1
8 CMV20 pos L5 M 21 no_us_ab 1
9 CMV21 neg NC1 F 21 no_us_ab 1
10 CMV1 pos M1 F 21 no_us_ab 1
# … with 16 more rowsQuality control
Genes that do not have an adequate number of reads in any sample should be filtered out prior to downstream analyses. From a biological perspective, genes that are not expressed at a biologically meaningful level in any condition are not of interest. Statistically, we get a better estimate of the mean-variance relationship in the data and reduce the number of statistical tests that are performed during differential expression analyses.
Filter out lowly expressed genes and genes without Entrez IDs and calculate TMM normalization factors.
z <- dge[, colnames(dge) %in% targets$id] # retain only relevant samples
z$genes$Ensembl <- strsplit2(z$genes$Geneid, ".", fixed = TRUE)[,1]
z$group <- targets$CMV_status
edb <- EnsDb.Hsapiens.v86 # add Gene Symbols and Entrez IDs
z$genes <- left_join(z$genes, ensembldb::genes(edb,
filter = GeneIdFilter(z$genes$Ensembl),
columns = c("gene_id",
"symbol",
"entrezid"),
return.type = "data.frame"),
by = c("Ensembl" = "gene_id"))
z$genes$entrezid <- unlist(sapply(z$genes$entrezid, function(x) {
if(is.null(x)) NA else x[length(x)]
}), use.names = FALSE)
keep <- !is.na(z$genes$entrezid) & !is.null(z$genes$entrezid)
x <- z[keep, ] # remove genes without Entrez IDs
keep <- filterByExpr(x, group = z$group)
x <- x[keep, ] # remove lowly expressed genes
y <- calcNormFactors(x)
yAn object of class "DGEList"
$counts
CMV30 CMV31 CMV8 CMV9 CMV26 CMV14 CMV15 CMV20 CMV21 CMV1 CMV2 CMV3 CMV4
32 20 36 28 42 28 25 19 26 18 32 25 55 31
52 88 73 55 43 55 75 53 64 55 61 68 36 62
55 6 15 17 15 15 9 14 13 12 7 10 16 12
63 148 172 148 126 175 179 176 141 179 155 194 121 123
64 14 15 14 11 20 14 9 13 14 13 13 19 10
CMV10 CMV11 CMV19 CMV35 CMV51 CMV52 CMV53 CMV54 CMV56 CMV57 CMV58 CMV60
32 15 11 22 24 25 19 23 20 18 23 24 8
52 69 11 35 47 59 49 42 37 38 49 65 49
55 7 2 11 9 19 9 6 12 7 7 7 5
63 156 44 107 164 144 71 71 137 136 122 131 88
64 9 3 6 9 14 3 9 12 5 7 12 5
CMV61
32 20
52 49
55 6
63 108
64 3
12727 more rows ...
$samples
group lib.size norm.factors
CMV30 1 4673630 1.017563
CMV31 1 5232010 1.052277
CMV8 1 3594801 1.052059
CMV9 1 3425478 1.026378
CMV26 1 4892776 1.068518
21 more rows ...
$genes
Geneid Length Ensembl symbol entrezid
32 ENSG00000230021.10 5495 ENSG00000230021 RP5-857K21.4 101928626
52 ENSG00000228794.10 15682 ENSG00000228794 LINC01128 643837
55 ENSG00000230368.2 1971 ENSG00000230368 FAM41C 284593
63 ENSG00000188976.11 5540 ENSG00000188976 NOC2L 26155
64 ENSG00000187961.14 3402 ENSG00000187961 KLHL17 339451
12727 more rows ...
$group
[1] "pos" "neg" "neg" "pos" "pos"
21 more elements ...Plotting the distribution log-CPM values shows that a majority of genes within each sample are either not expressed or lowly-expressed with log-CPM values that are small or negative.
L <- mean(z$samples$lib.size) * 1e-6
M <- median(z$samples$lib.size) * 1e-6
par(mfrow = c(1,2))
lcpmz <- cpm(z, log = TRUE)
lcpm.cutoff <- log2(10/M + 2/L)
nsamples <- ncol(z)
col <- scales::hue_pal()(nsamples)
plot(density(lcpmz[,1]), col = col[1], lwd = 1, ylim = c(0, 2), las = 2,
main = "", xlab = "")
title(main = "Unfiltered data", xlab = "Log-cpm")
abline(v = lcpm.cutoff, lty = 3)
for (i in 2:nsamples){
den <- density(lcpmz[,i])
lines(den$x, den$y, col = col[i], lwd = 1)
}
lcpmy <- cpm(y, log=TRUE)
plot(density(lcpmy[,1]), col = col[1], lwd = 1, ylim = c(0, 0.25), las = 2,
main = "", xlab = "")
title(main = "Filtered data", xlab = "Log-cpm")
abline(v = lcpm.cutoff, lty = 3)
for (i in 2:nsamples){
den <- density(lcpmy[,i])
lines(den$x, den$y, col = col[i], lwd = 1)
}
| Version | Author | Date |
|---|---|---|
| 10dcedf | Jovana Maksimovic | 2020-11-06 |
Although in excess of 30 million reads were obtained per sample, we can see that after mapping, duplicate removal and quantification of gene expression the median library size is just under than 4 million reads. This suggests that we are likely to only be capturing the most abundant cfRNAs.
It is assumed that all samples should have a similar range and distribution of expression values. The raw data looks fairly uniform between samples, although TMM normalization further improves this.
dat <- data.frame(lib = y$samples$lib.size,
status = y$group,
sample = colnames(y))
p1 <- ggplot(dat, aes(x = sample, y = lib, fill = status)) +
geom_bar(stat = "identity") +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(x = "Sample", y = "Library size",
fill = "CMV Status") +
geom_hline(yintercept = median(dat$lib), linetype = "dashed") +
scale_x_discrete(limits = dat$sample)
dat <- reshape2::melt(cpm(y, normalized.lib.sizes = FALSE, log = TRUE),
value.name = "cpm")
dat$status <- rep(y$group, each = nrow(y))
colnames(dat)[2] <- "sample"
p2 <- ggplot(dat, aes(x = sample, y = cpm, fill = status)) +
geom_boxplot(show.legend = FALSE, outlier.size = 0.75) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7)) +
labs(x = "Sample", y = "Library size",
fill = "CMV Status") +
geom_hline(yintercept = median(dat$lib), linetype = "dashed")
dat <- reshape2::melt(cpm(y, normalized.lib.sizes = TRUE, log = TRUE),
value.name = "cpm")
dat$status <- rep(y$group, each = nrow(y))
colnames(dat)[2] <- "sample"
p3 <- ggplot(dat, aes(x = sample, y = cpm, fill = status)) +
geom_boxplot(show.legend = FALSE, outlier.size = 0.75) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7)) +
labs(x = "Sample", y = "Library size",
fill = "CMV Status") +
geom_hline(yintercept = median(dat$lib), linetype = "dashed")
p1 / (p2 + p3) + plot_layout(guides = "collect")
Multi-dimensional scaling (MDS) plots show the largest sources of variation in the data. They are a good way of exploring the relationships between the samples and identifying structure in the data. The following series of MDS plots examines the first four principal components. The samples are coloured by various known features of the samples such as CMV Status and foetal sex. The MDS plots do not show the samples strongly clustering by any of the known features of the data, although there does seem to be some separation between the CMV positive and negative samples in the 1st and 2nd principal components. This indicates that there are possibly some differentially expressed genes between CMV positive and negative samples.
A weak batch effect is also evident in the 3rd principal component, when we examine the plots coloured by batch.
dims <- list(c(1,2), c(1,3), c(2,3), c(3,4))
vars <- c("CMV_status", "pair", "sex", "GA_at_amnio", "indication", "batch")
patches <- vector("list", length(vars))
for(i in 1:length(vars)){
p <- vector("list", length(dims))
for(j in 1:length(dims)){
mds <- plotMDS(cpm(y, log = TRUE), top = 1000, gene.selection="common",
plot = FALSE, dim.plot = dims[[j]])
dat <- tibble::tibble(x = mds$x, y = mds$y,
sample = targets$id,
variable = pull(targets, vars[i]))
p[[j]] <- ggplot(dat, aes(x = x, y = y, colour = variable)) +
geom_text(aes(label = sample), size = 2.5) +
labs(x = glue::glue("Principal component {dims[[j]][1]}"),
y = glue::glue("Principal component {dims[[j]][2]}"),
colour = vars[i])
}
patches[[i]] <- wrap_elements(wrap_plots(p, ncol = 2, guides = "collect") +
plot_annotation(title = glue::glue("Coloured by: {vars[i]}")) &
theme(legend.position = "bottom"))
}
wrap_plots(patches, ncol = 1)
Differential expression analysis
Due to the variability in the data, the TMM normalised data was transformed using voomWithQualityWeights. This takes into account the differing library sizes and the mean variance relationship in the data as well as calculating sample-specific quality weights. Linear models were fit in limma, taking into account the voom weights. The CMV positive samples were compared to the CMV negative samples, taking into account the sample pairs. A summary of the number of differentially expressed genes is shown below.
design <- model.matrix(~0 + y$group + targets$pair, data = targets)
colnames(design) <- c(levels(factor(y$group)), unique(targets$pair)[-1])
v <- voomWithQualityWeights(y, design, plot = TRUE)
cont <- makeContrasts(contrasts = "pos - neg",
levels = design)
fit <- lmFit(v, design)
cfit <- contrasts.fit(fit, cont)
fit2 <- eBayes(cfit, robust = TRUE)
fitSum <- summary(decideTests(fit2, p.value = 0.05))
fitSum pos - neg
Down 17
NotSig 12505
Up 210There were 17 down-regulated and 210 up-regulated genes between CMV positive and CMV negative samples at FDR < 0.05.
These are the top 10 differentially expressed genes.
top <- topTable(fit2, n = 500)
readr::write_csv(top, path = here("output/star-fc-limma-voom-all.csv"))Warning: The `path` argument of `write_csv()` is deprecated as of readr 1.4.0.
Please use the `file` argument instead.head(top, n=10) Geneid Length Ensembl symbol entrezid logFC
45751 ENSG00000140853.15 12386 ENSG00000140853 NLRC5 84166 2.797155
8687 ENSG00000115415.20 9770 ENSG00000115415 STAT1 6772 1.199505
19329 ENSG00000206337.12 11058 ENSG00000206337 HCP5 10866 4.440444
15204 ENSG00000137628.17 6746 ENSG00000137628 DDX60 55601 1.590297
8282 ENSG00000115267.8 5094 ENSG00000115267 IFIH1 64135 1.069293
14201 ENSG00000138646.9 4764 ENSG00000138646 HERC5 51191 2.943603
2095 ENSG00000137965.11 2038 ENSG00000137965 IFI44 10561 1.954652
46903 ENSG00000132530.17 6615 ENSG00000132530 XAF1 54739 2.858189
27479 ENSG00000107201.10 4640 ENSG00000107201 DDX58 23586 1.110450
19309 ENSG00000234745.11 3051 ENSG00000234745 HLA-B 3106 3.138597
AveExpr t P.Value adj.P.Val B
45751 2.163443 9.405942 4.679903e-08 0.0005193123 7.426915
8687 7.926887 8.204584 3.789399e-07 0.0012061657 6.855860
19329 1.845420 9.507407 8.157593e-08 0.0005193123 6.512669
15204 5.968763 7.878387 9.119445e-07 0.0020294388 6.033060
8282 5.614202 7.524000 9.868764e-07 0.0020294388 5.912610
14201 1.990993 8.409535 2.397649e-07 0.0010175624 5.881695
2095 4.394309 7.738868 1.138190e-06 0.0020294388 5.756492
46903 4.797937 7.488143 1.705713e-06 0.0024130150 5.359573
27479 5.639234 7.054829 2.755727e-06 0.0035085914 4.933196
19309 4.348936 7.058890 3.844793e-06 0.0041278571 4.573300The following plots show the expression of the top 12 ranked differentially expressed genes for CMV positive and CMV negative samples. Although there is significant variability within the groups and the log2 fold changes are not large, there are obvious differences in expression for the top ranked genes.
dat <- reshape2::melt(cpm(y, log = TRUE),
value.name = "cpm")
dat$status <- rep(y$group, each = nrow(y))
dat$gene <- rep(y$genes$Geneid, ncol(y))
p <- vector("list", 12)
for(i in 1:length(p)){
p[[i]] <- ggplot(data = subset(dat, dat$gene == top$Geneid[i]),
aes(x = status, y = cpm, colour = status)) +
geom_jitter(width = 0.25) +
stat_summary(fun = "mean", geom = "crossbar") +
labs(x = "Status", y = "log2 CPM", colour = "Status") +
ggtitle(top$symbol[i]) +
theme(plot.title = element_text(size = 8),
plot.subtitle = element_text(size = 7),
axis.title = element_text(size = 8),
axis.text.x = element_text(size = 7))
}
wrap_plots(p, guides = "collect", ncol = 3) &
theme(legend.position = "bottom")
topTable(fit2, num = Inf) %>%
mutate(sig = ifelse(adj.P.Val <= 0.05, "<= 0.05", "> 0.05")) -> dat
ggplot(dat, aes(x = logFC, y = -log10(P.Value), color = sig)) +
geom_point(alpha = 0.75) +
ggrepel::geom_text_repel(data = subset(dat, adj.P.Val < 0.05),
aes(x = logFC, y = -log10(P.Value),
label = symbol),
size = 2, colour = "black", max.overlaps = 15) +
labs(x = expression(~log[2]~"(Fold Change)"),
y = expression(~-log[10]~"(P-value)"),
colour = "FDR") +
scale_colour_brewer(palette = "Set1")Warning: ggrepel: 150 unlabeled data points (too many overlaps). Consider
increasing max.overlaps
cut <- 0.05
o <- order(targets$CMV_status)
aheatmap(v$E[v$genes$Geneid %in% top$Geneid[top$adj.P.Val < cut], o],
annCol = list(CMV_result = targets$CMV_status[o]),
Colv = NA,
labRow = v$genes$symbol[v$genes$Geneid %in%
top$Geneid[top$adj.P.Val < cut]],
main = glue::glue("DEGs at FDR < {cut}"))
Gene set enrichment analysis (GSEA)
Testing for enrichment of Gene Ontology (GO) categories among statistically significant differentially expressed genes.
go <- goana(top$entrezid[top$adj.P.Val < 0.05], universe = v$genes$entrezid,
trend = v$genes$Length)
topGO(go, number = Inf) %>%
mutate(FDR = p.adjust(P.DE)) %>%
dplyr::filter(FDR < 0.05) %>%
knitr::kable(format.args = list(scientific = -1), digits = 50)| Term | Ont | N | DE | P.DE | FDR | |
|---|---|---|---|---|---|---|
| GO:0009607 | response to biotic stimulus | BP | 876 | 62 | 9.627568e-22 | 1.987900e-17 |
| GO:0051607 | defense response to virus | BP | 178 | 30 | 5.069312e-21 | 1.046661e-16 |
| GO:0043207 | response to external biotic stimulus | BP | 849 | 60 | 5.634435e-21 | 1.163285e-16 |
| GO:0051707 | response to other organism | BP | 849 | 60 | 5.634435e-21 | 1.163285e-16 |
| GO:0045087 | innate immune response | BP | 516 | 47 | 8.260606e-21 | 1.705320e-16 |
| GO:0009615 | response to virus | BP | 231 | 33 | 1.026368e-20 | 2.118731e-16 |
| GO:0071357 | cellular response to type I interferon | BP | 70 | 21 | 1.410360e-20 | 2.911266e-16 |
| GO:0060337 | type I interferon signaling pathway | BP | 70 | 21 | 1.410360e-20 | 2.911266e-16 |
| GO:0034340 | response to type I interferon | BP | 74 | 21 | 5.205041e-20 | 1.074320e-15 |
| GO:0098542 | defense response to other organism | BP | 628 | 49 | 7.916235e-19 | 1.633832e-14 |
| GO:0006955 | immune response | BP | 1177 | 67 | 2.011052e-18 | 4.150408e-14 |
| GO:0006952 | defense response | BP | 913 | 57 | 2.286585e-17 | 4.718825e-13 |
| GO:0009605 | response to external stimulus | BP | 1674 | 76 | 1.969940e-15 | 4.065169e-11 |
| GO:0002376 | immune system process | BP | 1887 | 80 | 1.245515e-14 | 2.570121e-10 |
| GO:0044419 | interspecies interaction between organisms | BP | 1439 | 68 | 1.548258e-14 | 3.194675e-10 |
| GO:0034341 | response to interferon-gamma | BP | 121 | 20 | 3.483005e-14 | 7.186484e-10 |
| GO:0071346 | cellular response to interferon-gamma | BP | 109 | 19 | 5.693487e-14 | 1.174680e-09 |
| GO:0002252 | immune effector process | BP | 776 | 47 | 8.171695e-14 | 1.685902e-09 |
| GO:0060333 | interferon-gamma-mediated signaling pathway | BP | 61 | 15 | 1.316373e-13 | 2.715678e-09 |
| GO:0034097 | response to cytokine | BP | 802 | 46 | 1.113023e-12 | 2.296054e-08 |
| GO:0071345 | cellular response to cytokine stimulus | BP | 740 | 41 | 6.912916e-11 | 1.425996e-06 |
| GO:0019221 | cytokine-mediated signaling pathway | BP | 505 | 33 | 9.514700e-11 | 1.962597e-06 |
| GO:0050776 | regulation of immune response | BP | 513 | 33 | 1.436948e-10 | 2.963849e-06 |
| GO:0001817 | regulation of cytokine production | BP | 440 | 30 | 2.673853e-10 | 5.514823e-06 |
| GO:0001816 | cytokine production | BP | 483 | 31 | 5.835633e-10 | 1.203541e-05 |
| GO:0045071 | negative regulation of viral genome replication | BP | 48 | 11 | 5.949646e-10 | 1.226995e-05 |
| GO:0007166 | cell surface receptor signaling pathway | BP | 1939 | 71 | 7.588803e-10 | 1.564963e-05 |
| GO:1903900 | regulation of viral life cycle | BP | 126 | 16 | 8.102998e-10 | 1.670919e-05 |
| GO:0048525 | negative regulation of viral process | BP | 77 | 13 | 9.103396e-10 | 1.877120e-05 |
| GO:0031224 | intrinsic component of membrane | CC | 2637 | 86 | 1.895169e-09 | 3.907649e-05 |
| GO:0045824 | negative regulation of innate immune response | BP | 41 | 10 | 1.899346e-09 | 3.916071e-05 |
| GO:0045069 | regulation of viral genome replication | BP | 85 | 13 | 3.226057e-09 | 6.651162e-05 |
| GO:0050777 | negative regulation of immune response | BP | 87 | 13 | 4.329391e-09 | 8.925472e-05 |
| GO:0071310 | cellular response to organic substance | BP | 1817 | 66 | 5.472067e-09 | 1.128067e-04 |
| GO:0016021 | integral component of membrane | CC | 2563 | 83 | 6.204369e-09 | 1.278969e-04 |
| GO:0009617 | response to bacterium | BP | 335 | 24 | 7.175135e-09 | 1.479011e-04 |
| GO:0019079 | viral genome replication | BP | 109 | 14 | 8.208250e-09 | 1.691885e-04 |
| GO:0043903 | regulation of symbiotic process | BP | 190 | 18 | 8.356821e-09 | 1.722424e-04 |
| GO:0002682 | regulation of immune system process | BP | 909 | 42 | 9.784380e-09 | 2.016561e-04 |
| GO:1903901 | negative regulation of viral life cycle | BP | 63 | 11 | 1.281074e-08 | 2.640166e-04 |
| GO:0031347 | regulation of defense response | BP | 432 | 27 | 1.458622e-08 | 3.005929e-04 |
| GO:0050792 | regulation of viral process | BP | 178 | 17 | 1.950967e-08 | 4.020358e-04 |
| GO:0010033 | response to organic substance | BP | 2227 | 74 | 2.204467e-08 | 4.542524e-04 |
| GO:0002831 | regulation of response to biotic stimulus | BP | 277 | 21 | 2.482535e-08 | 5.115264e-04 |
| GO:0001818 | negative regulation of cytokine production | BP | 160 | 16 | 2.691509e-08 | 5.545585e-04 |
| GO:0070887 | cellular response to chemical stimulus | BP | 2197 | 73 | 2.896441e-08 | 5.967537e-04 |
| GO:0002697 | regulation of immune effector process | BP | 231 | 19 | 3.231642e-08 | 6.657828e-04 |
| GO:0046977 | TAP binding | MF | 7 | 5 | 3.575728e-08 | 7.366358e-04 |
| GO:0048584 | positive regulation of response to stimulus | BP | 1500 | 56 | 4.649513e-08 | 9.577996e-04 |
| GO:0044403 | symbiotic process | BP | 828 | 38 | 6.670711e-08 | 1.374100e-03 |
| GO:0005887 | integral component of plasma membrane | CC | 667 | 33 | 9.779821e-08 | 2.014448e-03 |
| GO:0031226 | intrinsic component of plasma membrane | CC | 705 | 34 | 1.126292e-07 | 2.319824e-03 |
| GO:0005903 | brush border | CC | 82 | 11 | 2.161947e-07 | 4.452745e-03 |
| GO:0042221 | response to chemical | BP | 2816 | 84 | 2.305232e-07 | 4.747625e-03 |
| GO:0003725 | double-stranded RNA binding | MF | 66 | 10 | 2.417977e-07 | 4.979581e-03 |
| GO:0002832 | negative regulation of response to biotic stimulus | BP | 68 | 10 | 3.229737e-07 | 6.650997e-03 |
| GO:1904970 | brush border assembly | BP | 5 | 4 | 4.917484e-07 | 1.012608e-02 |
| GO:0016032 | viral process | BP | 793 | 35 | 6.002323e-07 | 1.235938e-02 |
| GO:0019058 | viral life cycle | BP | 279 | 19 | 6.247614e-07 | 1.286384e-02 |
| GO:0035455 | response to interferon-alpha | BP | 19 | 6 | 6.848944e-07 | 1.410129e-02 |
| GO:0050896 | response to stimulus | BP | 5571 | 136 | 7.291092e-07 | 1.501090e-02 |
| GO:0007154 | cell communication | BP | 3936 | 105 | 8.632247e-07 | 1.777121e-02 |
| GO:0032101 | regulation of response to external stimulus | BP | 631 | 30 | 9.222986e-07 | 1.898644e-02 |
| GO:0051241 | negative regulation of multicellular organismal process | BP | 772 | 34 | 9.334281e-07 | 1.921462e-02 |
| GO:0002703 | regulation of leukocyte mediated immunity | BP | 115 | 12 | 9.769944e-07 | 2.011045e-02 |
| GO:0045088 | regulation of innate immune response | BP | 208 | 16 | 1.007157e-06 | 2.073031e-02 |
| GO:0023052 | signaling | BP | 3907 | 104 | 1.154867e-06 | 2.376947e-02 |
| GO:0007165 | signal transduction | BP | 3630 | 98 | 1.628915e-06 | 3.352471e-02 |
| GO:0048583 | regulation of response to stimulus | BP | 2814 | 81 | 1.992047e-06 | 4.099633e-02 |
| GO:0002483 | antigen processing and presentation of endogenous peptide antigen | BP | 13 | 5 | 2.006796e-06 | 4.129785e-02 |
| GO:0019885 | antigen processing and presentation of endogenous peptide antigen via MHC class I | BP | 13 | 5 | 2.006796e-06 | 4.129785e-02 |
GSEA helps us to interpret the results of a differential expression analysis. The camera function performs a competitive test to assess whether the genes in a given set are highly ranked in terms of differential expression relative to genes that are not in the set. We have tested several collections of gene sets from the Broad Institute’s Molecular Signatures Database MSigDB.
Build gene set indexes.
gsAnnots <- buildIdx(entrezIDs = v$genes$entrezid, species = "human",
msigdb.gsets = c("h", "c2", "c5"))[1] "Loading MSigDB Gene Sets ... "
[1] "Loaded gene sets for the collection h ..."
[1] "Indexed the collection h ..."
[1] "Created annotation for the collection h ..."
[1] "Loaded gene sets for the collection c2 ..."
[1] "Indexed the collection c2 ..."
[1] "Created annotation for the collection c2 ..."
[1] "Loaded gene sets for the collection c5 ..."
[1] "Indexed the collection c5 ..."
[1] "Created annotation for the collection c5 ..."
[1] "Building KEGG pathways annotation object ... "The GO gene sets consist of genes annotated by the same GO terms.
c5Cam <- camera(v, gsAnnots$c5@idx, design, contrast = cont, trend.var = TRUE)
write.csv(c5Cam[c5Cam$FDR < 0.05,],
file = here("output/star-fc-limma-voom-all-gsea-c5.csv"))
head(c5Cam, n = 20) NGenes
GO_RESPONSE_TO_TYPE_I_INTERFERON 48
GO_DEFENSE_RESPONSE_TO_VIRUS 111
GO_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_ENDOGENOUS_PEPTIDE_ANTIGEN 11
GO_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_ENDOGENOUS_ANTIGEN 13
GO_INTERFERON_GAMMA_MEDIATED_SIGNALING_PATHWAY 43
GO_DIGESTION 56
GO_RESPONSE_TO_INTERFERON_GAMMA 71
GO_CELLULAR_RESPONSE_TO_INTERFERON_GAMMA 59
GO_REGULATION_OF_CELL_KILLING 32
GO_REGULATION_OF_LEUKOCYTE_MEDIATED_CYTOTOXICITY 29
GO_RESPONSE_TO_VIRUS 166
GO_SOLUTE_PROTON_SYMPORTER_ACTIVITY 18
GO_CYTOKINE_MEDIATED_SIGNALING_PATHWAY 245
GO_NEGATIVE_REGULATION_OF_CELL_KILLING 12
GO_MHC_PROTEIN_COMPLEX 11
GO_REGULATION_OF_T_CELL_MEDIATED_CYTOTOXICITY 12
GO_RESPONSE_TO_INTERFERON_ALPHA 18
GO_POSITIVE_REGULATION_OF_T_CELL_MEDIATED_CYTOTOXICITY 8
GO_REGULATION_OF_PLASMA_LIPOPROTEIN_PARTICLE_LEVELS 29
GO_MHC_CLASS_I_PROTEIN_COMPLEX 7
Direction
GO_RESPONSE_TO_TYPE_I_INTERFERON Up
GO_DEFENSE_RESPONSE_TO_VIRUS Up
GO_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_ENDOGENOUS_PEPTIDE_ANTIGEN Up
GO_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_ENDOGENOUS_ANTIGEN Up
GO_INTERFERON_GAMMA_MEDIATED_SIGNALING_PATHWAY Up
GO_DIGESTION Up
GO_RESPONSE_TO_INTERFERON_GAMMA Up
GO_CELLULAR_RESPONSE_TO_INTERFERON_GAMMA Up
GO_REGULATION_OF_CELL_KILLING Up
GO_REGULATION_OF_LEUKOCYTE_MEDIATED_CYTOTOXICITY Up
GO_RESPONSE_TO_VIRUS Up
GO_SOLUTE_PROTON_SYMPORTER_ACTIVITY Up
GO_CYTOKINE_MEDIATED_SIGNALING_PATHWAY Up
GO_NEGATIVE_REGULATION_OF_CELL_KILLING Up
GO_MHC_PROTEIN_COMPLEX Up
GO_REGULATION_OF_T_CELL_MEDIATED_CYTOTOXICITY Up
GO_RESPONSE_TO_INTERFERON_ALPHA Up
GO_POSITIVE_REGULATION_OF_T_CELL_MEDIATED_CYTOTOXICITY Up
GO_REGULATION_OF_PLASMA_LIPOPROTEIN_PARTICLE_LEVELS Up
GO_MHC_CLASS_I_PROTEIN_COMPLEX Up
PValue
GO_RESPONSE_TO_TYPE_I_INTERFERON 5.918561e-26
GO_DEFENSE_RESPONSE_TO_VIRUS 5.511815e-14
GO_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_ENDOGENOUS_PEPTIDE_ANTIGEN 4.049629e-13
GO_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_ENDOGENOUS_ANTIGEN 4.365482e-13
GO_INTERFERON_GAMMA_MEDIATED_SIGNALING_PATHWAY 6.081715e-13
GO_DIGESTION 4.430786e-11
GO_RESPONSE_TO_INTERFERON_GAMMA 6.699823e-11
GO_CELLULAR_RESPONSE_TO_INTERFERON_GAMMA 6.410894e-10
GO_REGULATION_OF_CELL_KILLING 1.371315e-08
GO_REGULATION_OF_LEUKOCYTE_MEDIATED_CYTOTOXICITY 2.794535e-08
GO_RESPONSE_TO_VIRUS 2.890393e-08
GO_SOLUTE_PROTON_SYMPORTER_ACTIVITY 5.791605e-08
GO_CYTOKINE_MEDIATED_SIGNALING_PATHWAY 9.989903e-08
GO_NEGATIVE_REGULATION_OF_CELL_KILLING 1.023995e-07
GO_MHC_PROTEIN_COMPLEX 2.027142e-07
GO_REGULATION_OF_T_CELL_MEDIATED_CYTOTOXICITY 2.805316e-07
GO_RESPONSE_TO_INTERFERON_ALPHA 3.328300e-07
GO_POSITIVE_REGULATION_OF_T_CELL_MEDIATED_CYTOTOXICITY 3.395518e-07
GO_REGULATION_OF_PLASMA_LIPOPROTEIN_PARTICLE_LEVELS 3.468113e-07
GO_MHC_CLASS_I_PROTEIN_COMPLEX 3.535987e-07
FDR
GO_RESPONSE_TO_TYPE_I_INTERFERON 3.646426e-22
GO_DEFENSE_RESPONSE_TO_VIRUS 1.697915e-10
GO_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_ENDOGENOUS_PEPTIDE_ANTIGEN 6.723934e-10
GO_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_ENDOGENOUS_ANTIGEN 6.723934e-10
GO_INTERFERON_GAMMA_MEDIATED_SIGNALING_PATHWAY 7.493889e-10
GO_DIGESTION 4.549679e-08
GO_RESPONSE_TO_INTERFERON_GAMMA 5.896801e-08
GO_CELLULAR_RESPONSE_TO_INTERFERON_GAMMA 4.937189e-07
GO_REGULATION_OF_CELL_KILLING 9.387410e-06
GO_REGULATION_OF_LEUKOCYTE_MEDIATED_CYTOTOXICITY 1.618883e-05
GO_RESPONSE_TO_VIRUS 1.618883e-05
GO_SOLUTE_PROTON_SYMPORTER_ACTIVITY 2.973507e-05
GO_CYTOKINE_MEDIATED_SIGNALING_PATHWAY 4.506309e-05
GO_NEGATIVE_REGULATION_OF_CELL_KILLING 4.506309e-05
GO_MHC_PROTEIN_COMPLEX 8.326149e-05
GO_REGULATION_OF_T_CELL_MEDIATED_CYTOTOXICITY 1.080222e-04
GO_RESPONSE_TO_INTERFERON_ALPHA 1.089261e-04
GO_POSITIVE_REGULATION_OF_T_CELL_MEDIATED_CYTOTOXICITY 1.089261e-04
GO_REGULATION_OF_PLASMA_LIPOPROTEIN_PARTICLE_LEVELS 1.089261e-04
GO_MHC_CLASS_I_PROTEIN_COMPLEX 1.089261e-04The Hallmark gene sets are coherently expressed signatures derived by aggregating many MSigDB gene sets to represent well-defined biological states or processes.
hCam <- camera(v, gsAnnots$h@idx, design, contrast = cont, trend.var = TRUE)
head(hCam, n = 20) NGenes Direction PValue FDR
HALLMARK_INTERFERON_ALPHA_RESPONSE 81 Up 4.065265e-46 2.032632e-44
HALLMARK_INTERFERON_GAMMA_RESPONSE 147 Up 5.643640e-33 1.410910e-31
HALLMARK_E2F_TARGETS 196 Down 5.883620e-08 9.806033e-07
HALLMARK_INFLAMMATORY_RESPONSE 116 Up 8.645708e-08 1.080713e-06
HALLMARK_MYC_TARGETS_V1 198 Down 1.853964e-07 1.853964e-06
HALLMARK_TNFA_SIGNALING_VIA_NFKB 162 Up 2.346065e-06 1.955054e-05
HALLMARK_KRAS_SIGNALING_UP 133 Up 3.327963e-06 2.377116e-05
HALLMARK_IL6_JAK_STAT3_SIGNALING 47 Up 4.590363e-05 2.868977e-04
HALLMARK_COMPLEMENT 140 Up 1.187809e-04 6.598941e-04
HALLMARK_G2M_CHECKPOINT 194 Down 1.898061e-04 9.490303e-04
HALLMARK_ALLOGRAFT_REJECTION 102 Up 3.379307e-04 1.536049e-03
HALLMARK_MITOTIC_SPINDLE 192 Down 1.686661e-03 7.027753e-03
HALLMARK_MYC_TARGETS_V2 54 Down 9.271710e-03 3.566042e-02
HALLMARK_DNA_REPAIR 145 Down 1.106814e-02 3.952908e-02
HALLMARK_XENOBIOTIC_METABOLISM 139 Up 2.211981e-02 7.373271e-02
HALLMARK_COAGULATION 84 Up 2.375926e-02 7.424769e-02
HALLMARK_IL2_STAT5_SIGNALING 142 Up 3.101322e-02 9.121536e-02
HALLMARK_ADIPOGENESIS 182 Up 4.278309e-02 1.188419e-01
HALLMARK_MYOGENESIS 122 Down 5.072696e-02 1.334920e-01
HALLMARK_ANDROGEN_RESPONSE 96 Up 5.671960e-02 1.417990e-01Barcode plots show the enrichment of gene sets among up or down-regulated genes. The following barcode plots show the enrichment of the top 4 hallmark gene sets among the genes differentially expressed between CMV positive and CMV negative samples.
par(mfrow=c(2,2))
sapply(1:4, function(i){
barcodeplot(fit2$t[,1], gsAnnots$h@idx[[rownames(hCam)[i]]],
main = rownames(hCam)[i], cex.main = 0.75)
})
[[1]]
NULL
[[2]]
NULL
[[3]]
NULL
[[4]]
NULLThe curated gene sets are compiled from online pathway databases, publications in PubMed, and knowledge of domain experts.
c2Cam <- camera(v, gsAnnots$c2@idx, design, contrast = cont, trend.var = TRUE)
write.csv(c2Cam[c2Cam$FDR < 0.05,],
file = here("output/star-fc-limma-voom-all-gsea-c2.csv"))
head(c2Cam, n = 20) NGenes Direction PValue
MOSERLE_IFNA_RESPONSE 27 Up 4.452715e-42
BROWNE_INTERFERON_RESPONSIVE_GENES 55 Up 2.853789e-40
HECKER_IFNB1_TARGETS 59 Up 1.074420e-36
SANA_RESPONSE_TO_IFNG_UP 48 Up 4.805722e-36
FARMER_BREAST_CANCER_CLUSTER_1 15 Up 2.184050e-27
DAUER_STAT3_TARGETS_DN 49 Up 8.940232e-27
BOSCO_INTERFERON_INDUCED_ANTIVIRAL_MODULE 59 Up 4.536072e-26
SANA_TNF_SIGNALING_UP 63 Up 1.091778e-24
EINAV_INTERFERON_SIGNATURE_IN_CANCER 24 Up 5.226263e-22
REACTOME_INTERFERON_ALPHA_BETA_SIGNALING 44 Up 7.206877e-22
BOWIE_RESPONSE_TO_TAMOXIFEN 16 Up 7.637944e-21
RADAEVA_RESPONSE_TO_IFNA1_UP 45 Up 1.858830e-20
BENNETT_SYSTEMIC_LUPUS_ERYTHEMATOSUS 24 Up 2.002087e-20
BOWIE_RESPONSE_TO_EXTRACELLULAR_MATRIX 16 Up 1.543851e-19
ZHANG_INTERFERON_RESPONSE 18 Up 1.924959e-18
DER_IFN_ALPHA_RESPONSE_UP 67 Up 9.637826e-18
SEITZ_NEOPLASTIC_TRANSFORMATION_BY_8P_DELETION_UP 56 Up 1.992852e-17
KRASNOSELSKAYA_ILF3_TARGETS_UP 30 Up 4.035282e-17
TAKEDA_TARGETS_OF_NUP98_HOXA9_FUSION_3D_UP 129 Up 1.819293e-16
ROETH_TERT_TARGETS_UP 13 Up 5.131494e-16
FDR
MOSERLE_IFNA_RESPONSE 1.665761e-38
BROWNE_INTERFERON_RESPONSIVE_GENES 5.338013e-37
HECKER_IFNB1_TARGETS 1.339802e-33
SANA_RESPONSE_TO_IFNG_UP 4.494551e-33
FARMER_BREAST_CANCER_CLUSTER_1 1.634106e-24
DAUER_STAT3_TARGETS_DN 5.574235e-24
BOSCO_INTERFERON_INDUCED_ANTIVIRAL_MODULE 2.424206e-23
SANA_TNF_SIGNALING_UP 5.105427e-22
EINAV_INTERFERON_SIGNATURE_IN_CANCER 2.172383e-19
REACTOME_INTERFERON_ALPHA_BETA_SIGNALING 2.696093e-19
BOWIE_RESPONSE_TO_TAMOXIFEN 2.597595e-18
RADAEVA_RESPONSE_TO_IFNA1_UP 5.761389e-18
BENNETT_SYSTEMIC_LUPUS_ERYTHEMATOSUS 5.761389e-18
BOWIE_RESPONSE_TO_EXTRACELLULAR_MATRIX 4.125391e-17
ZHANG_INTERFERON_RESPONSE 4.800848e-16
DER_IFN_ALPHA_RESPONSE_UP 2.253444e-15
SEITZ_NEOPLASTIC_TRANSFORMATION_BY_8P_DELETION_UP 4.385447e-15
KRASNOSELSKAYA_ILF3_TARGETS_UP 8.386661e-15
TAKEDA_TARGETS_OF_NUP98_HOXA9_FUSION_3D_UP 3.582092e-14
ROETH_TERT_TARGETS_UP 9.598459e-14The following barcode plots show the enrichment of the top 4 curated gene sets among the genes differentially expressed between CMV positive and CMV negative samples.
par(mfrow=c(2,2))
sapply(1:4, function(i){
barcodeplot(fit2$t[,1], gsAnnots$c2@idx[[rownames(c2Cam)[i]]],
main = rownames(c2Cam)[i], cex.main = 0.75)
})
[[1]]
NULL
[[2]]
NULL
[[3]]
NULL
[[4]]
NULLThe KEGG gene sets encompass all of the pathways defined in the Kegg pathway database.
keggCam <- camera(v, gsAnnots$kegg@idx, design, contrast = cont,
trend.var = TRUE)
head(keggCam, n = 20) NGenes Direction PValue
Mineral absorption 36 Up 1.375770e-10
Fat digestion and absorption 21 Up 4.933923e-08
Type I diabetes mellitus 12 Up 4.088334e-07
Antigen processing and presentation 39 Up 9.194603e-07
Graft-versus-host disease 8 Up 1.542188e-06
Autoimmune thyroid disease 10 Up 2.363013e-06
Cytokine-cytokine receptor interaction 83 Up 2.773616e-06
Allograft rejection 9 Up 7.401576e-06
Ribosome 125 Down 1.150378e-05
Spliceosome 128 Down 1.501836e-05
Vitamin digestion and absorption 19 Up 1.620295e-05
NOD-like receptor signaling pathway 117 Up 1.756741e-05
Bile secretion 34 Up 2.845918e-05
Hepatitis C 98 Up 4.236596e-05
Influenza A 119 Up 1.399519e-04
Measles 89 Up 2.186454e-04
Maturity onset diabetes of the young 9 Up 3.538335e-04
Cell adhesion molecules (CAMs) 54 Up 3.720338e-04
Mismatch repair 22 Down 4.258843e-04
Leishmaniasis 39 Up 4.565431e-04
FDR
Mineral absorption 4.003492e-08
Fat digestion and absorption 7.178858e-06
Type I diabetes mellitus 3.965684e-05
Antigen processing and presentation 6.689074e-05
Graft-versus-host disease 8.975534e-05
Autoimmune thyroid disease 1.146062e-04
Cytokine-cytokine receptor interaction 1.153032e-04
Allograft rejection 2.692323e-04
Ribosome 3.719557e-04
Spliceosome 4.260096e-04
Vitamin digestion and absorption 4.260096e-04
NOD-like receptor signaling pathway 4.260096e-04
Bile secretion 6.370478e-04
Hepatitis C 8.806067e-04
Influenza A 2.715068e-03
Measles 3.976612e-03
Maturity onset diabetes of the young 6.014546e-03
Cell adhesion molecules (CAMs) 6.014546e-03
Mismatch repair 6.522754e-03
Leishmaniasis 6.642703e-03par(mfrow=c(2,2))
sapply(1:4, function(i){
barcodeplot(fit2$t[,1], gsAnnots$kegg@idx[[rownames(keggCam)[i]]],
main = rownames(keggCam)[i], cex.main = 0.75)
})
[[1]]
NULL
[[2]]
NULL
[[3]]
NULL
[[4]]
NULLBrain development genes
Test only the specialised brain development gene set for differential expression between CMV positive and CMV negative samples. This reduces the multiple testing burden and can identify DEGs from a particular set of interest.
brainSet <- readr::read_delim(file = here("data/brain-development-geneset.txt"),
delim = "\t", skip = 2, col_names = "BRAIN_DEV")
brainSet# A tibble: 51 x 1
BRAIN_DEV
<chr>
1 AC139768.1
2 ADGRG1
3 AFF2
4 ALK
5 ALX1
6 BPTF
7 CDK5R1
8 CEP290
9 CLN5
10 CNTN4
# … with 41 more rowsfit <- lmFit(v[v$genes$symbol %in% brainSet$BRAIN_DEV, ], design)
cfit <- contrasts.fit(fit, cont)
fit2 <- eBayes(cfit, robust = TRUE)
fitSum <- summary(decideTests(fit2, p.value = 0.05))
fitSum pos - neg
Down 1
NotSig 26
Up 0topBrain <- topTable(fit2, n = Inf)
topBrain Geneid Length Ensembl symbol entrezid logFC
54718 ENSG00000053438.11 1329 ENSG00000053438 NNAT 4826 -0.658397816
46775 ENSG00000040531.16 4949 ENSG00000040531 CTNS 1497 0.207087793
47853 ENSG00000176749.9 3948 ENSG00000176749 CDK5R1 8851 0.197638751
59499 ENSG00000102038.15 4363 ENSG00000102038 SMARCA1 6594 -0.169096017
28963 ENSG00000167081.18 4904 ENSG00000167081 PBX3 5090 0.162490669
7681 ENSG00000074047.21 7341 ENSG00000074047 GLI2 2736 -0.136586759
59824 ENSG00000155966.14 14241 ENSG00000155966 AFF2 2334 -0.114952294
13911 ENSG00000132467.4 2020 ENSG00000132467 UTP3 57050 0.103859874
3760 ENSG00000185630.19 21370 ENSG00000185630 PBX1 5087 -0.099918011
51917 ENSG00000130479.11 4933 ENSG00000130479 MAP1S 55201 0.089549092
9697 ENSG00000144619.15 9123 ENSG00000144619 CNTN4 152330 0.015179716
45778 ENSG00000205336.13 8659 ENSG00000205336 ADGRG1 9289 0.131955954
58386 ENSG00000158352.15 10333 ENSG00000158352 SHROOM4 57477 0.080489267
24536 ENSG00000164690.8 5234 ENSG00000164690 SHH 6469 0.069843390
39080 ENSG00000102805.16 22814 ENSG00000102805 CLN5 1203 -0.063771691
42132 ENSG00000114062.21 12808 ENSG00000114062 UBE3A 7337 0.050810995
10970 ENSG00000185008.17 16663 ENSG00000185008 ROBO2 6092 -0.036031512
33637 ENSG00000110697.13 6442 ENSG00000110697 PITPNM1 9600 0.089169404
2823 ENSG00000092621.12 9217 ENSG00000092621 PHGDH 26227 0.030791118
1293 ENSG00000131238.17 5088 ENSG00000131238 PPT1 5538 0.034873502
16275 ENSG00000164258.12 1174 ENSG00000164258 NDUFS4 4724 0.019761214
49155 ENSG00000171634.18 15119 ENSG00000171634 BPTF 2186 0.022251393
57747 ENSG00000146950.13 8218 ENSG00000146950 SHROOM2 357 -0.025131407
23659 ENSG00000128573.26 16334 ENSG00000128573 FOXP2 93986 0.040384825
47791 ENSG00000196712.18 27130 ENSG00000196712 NF1 4763 0.001432021
37089 ENSG00000198707.16 10442 ENSG00000198707 CEP290 80184 -0.009571616
54624 ENSG00000198646.14 11133 ENSG00000198646 NCOA6 23054 -0.006603637
AveExpr t P.Value adj.P.Val B
54718 6.704590 -4.03247492 0.0008889658 0.02400208 -0.5600968
46775 2.783450 2.12484119 0.0481112739 0.64950220 -3.8586145
47853 3.817792 1.70114074 0.1074223665 0.81797536 -4.7434003
59499 7.397718 -1.61608178 0.1247714223 0.81797536 -5.0555150
28963 4.749400 1.42912478 0.1713709361 0.81797536 -5.1513516
7681 2.125700 -0.92352083 0.3682863535 0.85268245 -5.1713022
59824 1.620371 -0.51390733 0.6137403741 0.95281611 -5.2117635
13911 4.929279 1.33379242 0.1993605308 0.81797536 -5.2983565
3760 8.662670 -1.28478125 0.2156034668 0.81797536 -5.4118594
51917 4.779260 1.16607815 0.2592187731 0.81797536 -5.5049297
9697 2.365792 0.10097458 0.9207177180 0.97006631 -5.5409852
45778 4.413092 0.99341767 0.3346552258 0.85268245 -5.6086049
58386 5.385050 1.13258409 0.2726584540 0.81797536 -5.6225705
24536 3.394605 0.51655473 0.6119286850 0.95281611 -5.6786584
39080 3.405381 -0.73383954 0.4727578109 0.91174721 -5.7417837
42132 7.708283 0.90263101 0.3789699778 0.85268245 -5.8709223
10970 3.886193 -0.24404785 0.8101607935 0.95281611 -5.9016656
33637 5.297035 0.75961566 0.4580529238 0.91174721 -5.9082439
2823 3.987221 0.27427429 0.7870829789 0.95281611 -5.9084511
1293 4.185204 0.38936485 0.7017004855 0.95281611 -5.9459601
16275 4.834618 0.24191244 0.8116581647 0.95281611 -6.0689442
49155 9.098445 0.33664492 0.7403851592 0.95281611 -6.1661418
57747 5.863156 -0.40069222 0.6934930242 0.95281611 -6.2123748
23659 5.960003 0.28640956 0.7780801393 0.95281611 -6.2473473
47791 7.809457 0.01968990 0.9845133510 0.98451335 -6.2771041
37089 6.668566 -0.08389176 0.9341379292 0.97006631 -6.2818256
54624 7.381450 -0.13412714 0.8948314579 0.97006631 -6.2842430The following plots show the expression of the top 9 genes from the brain development set as ranked by their differential expression with regard to CMV positive and CMV negative status.
b <- y[y$genes$entrezid %in% topBrain$entrezid[1:9], ]
dat <- reshape2::melt(cpm(b, log = TRUE),
value.name = "cpm")
dat$status <- rep(b$group, each = nrow(b))
dat$gene <- rep(b$genes$Geneid, ncol(b))
p <- vector("list", 9)
for(i in 1:length(p)){
p[[i]] <- ggplot(data = subset(dat, dat$gene == topBrain$Geneid[i]),
aes(x = status, y = cpm, colour = status)) +
geom_jitter(width = 0.25) +
stat_summary(fun = "mean", geom = "crossbar") +
labs(x = "Status", y = "log2 CPM", colour = "Status") +
ggtitle(topBrain$symbol[i]) +
theme(plot.title = element_text(size = 8),
plot.subtitle = element_text(size = 7),
axis.title = element_text(size = 8),
axis.text.x = element_text(size = 7))
}
wrap_plots(p, guides = "collect", ncol = 3) &
theme(legend.position = "bottom")
Summary
Although the effective library sizes were low, the data is generally of good quality. We found at total of 0 differentially expressed genes at FDR < 0.05. The significant genes were enriched for GO terms associated with interferon response and similar. Further gene set testing results indicate an upregulation of interferon response genes in the CMV positive samples, relative to the CMV negative samples, which is consistent with the top genes from the DE analysis.
sessionInfo()R version 4.0.2 (2020-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS: /config/binaries/R/4.0.2/lib64/R/lib/libRblas.so
LAPACK: /config/binaries/R/4.0.2/lib64/R/lib/libRlapack.so
locale:
[1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_AU.UTF-8 LC_COLLATE=en_AU.UTF-8
[5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8
[7] LC_PAPER=en_AU.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] RColorBrewer_1.1-2 EGSEA_1.18.1
[3] pathview_1.30.1 topGO_2.42.0
[5] SparseM_1.78 GO.db_3.12.1
[7] graph_1.68.0 gage_2.40.1
[9] patchwork_1.1.1 NMF_0.23.0
[11] cluster_2.1.0 rngtools_1.5
[13] pkgmaker_0.32.2 registry_0.5-1
[15] edgeR_3.32.1 limma_3.46.0
[17] EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.14.0
[19] AnnotationFilter_1.14.0 GenomicFeatures_1.42.1
[21] AnnotationDbi_1.52.0 Biobase_2.50.0
[23] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7
[25] IRanges_2.24.1 S4Vectors_0.28.1
[27] BiocGenerics_0.36.1 forcats_0.5.1
[29] stringr_1.4.0 dplyr_1.0.4
[31] purrr_0.3.4 readr_1.4.0
[33] tidyr_1.1.2 tibble_3.1.2
[35] ggplot2_3.3.5 tidyverse_1.3.0
[37] here_1.0.1 workflowr_1.6.2
loaded via a namespace (and not attached):
[1] rappdirs_0.3.3 rtracklayer_1.50.0
[3] Glimma_2.0.0 bit64_4.0.5
[5] knitr_1.31 multcomp_1.4-16
[7] DelayedArray_0.16.3 data.table_1.13.6
[9] hwriter_1.3.2 KEGGREST_1.30.1
[11] RCurl_1.98-1.3 doParallel_1.0.16
[13] generics_0.1.0 metap_1.4
[15] org.Mm.eg.db_3.12.0 TH.data_1.0-10
[17] RSQLite_2.2.5 bit_4.0.4
[19] mutoss_0.1-12 xml2_1.3.2
[21] lubridate_1.7.9.2 httpuv_1.5.5
[23] SummarizedExperiment_1.20.0 assertthat_0.2.1
[25] xfun_0.23 hms_1.0.0
[27] evaluate_0.14 promises_1.2.0.1
[29] fansi_0.5.0 progress_1.2.2
[31] caTools_1.18.1 dbplyr_2.1.0
[33] readxl_1.3.1 Rgraphviz_2.34.0
[35] DBI_1.1.1 geneplotter_1.68.0
[37] tmvnsim_1.0-2 htmlwidgets_1.5.3
[39] ellipsis_0.3.2 backports_1.2.1
[41] annotate_1.68.0 PADOG_1.32.0
[43] gbRd_0.4-11 gridBase_0.4-7
[45] biomaRt_2.46.3 MatrixGenerics_1.2.1
[47] HTMLUtils_0.1.7 vctrs_0.3.8
[49] cachem_1.0.4 withr_2.4.2
[51] globaltest_5.44.0 GenomicAlignments_1.26.0
[53] prettyunits_1.1.1 mnormt_2.0.2
[55] lazyeval_0.2.2 crayon_1.4.1
[57] genefilter_1.72.1 labeling_0.4.2
[59] pkgconfig_2.0.3 nlme_3.1-152
[61] ProtGenerics_1.22.0 GSA_1.03.1
[63] rlang_0.4.11 lifecycle_1.0.0
[65] sandwich_3.0-0 BiocFileCache_1.14.0
[67] mathjaxr_1.2-0 modelr_0.1.8
[69] cellranger_1.1.0 rprojroot_2.0.2
[71] GSVA_1.38.2 matrixStats_0.59.0
[73] Matrix_1.3-2 zoo_1.8-9
[75] reprex_1.0.0 whisker_0.4
[77] png_0.1-7 viridisLite_0.4.0
[79] bitops_1.0-7 KernSmooth_2.23-18
[81] Biostrings_2.58.0 blob_1.2.1
[83] R2HTML_2.3.2 doRNG_1.8.2
[85] scales_1.1.1 memoise_2.0.0.9000
[87] GSEABase_1.52.1 magrittr_2.0.1
[89] plyr_1.8.6 safe_3.30.0
[91] gplots_3.1.1 zlibbioc_1.36.0
[93] compiler_4.0.2 plotrix_3.8-1
[95] KEGGgraph_1.50.0 DESeq2_1.30.1
[97] Rsamtools_2.6.0 cli_3.0.0
[99] XVector_0.30.0 EGSEAdata_1.18.0
[101] MASS_7.3-53.1 tidyselect_1.1.0
[103] stringi_1.5.3 highr_0.8
[105] yaml_2.2.1 askpass_1.1
[107] locfit_1.5-9.4 ggrepel_0.9.1
[109] grid_4.0.2 tools_4.0.2
[111] rstudioapi_0.13 foreach_1.5.1
[113] git2r_0.28.0 farver_2.1.0
[115] digest_0.6.27 Rcpp_1.0.6
[117] broom_0.7.4 later_1.1.0.1
[119] org.Hs.eg.db_3.12.0 httr_1.4.2
[121] Rdpack_2.1 colorspace_2.0-2
[123] rvest_0.3.6 XML_3.99-0.5
[125] fs_1.5.0 splines_4.0.2
[127] statmod_1.4.35 sn_1.6-2
[129] multtest_2.46.0 plotly_4.9.3
[131] xtable_1.8-4 jsonlite_1.7.2
[133] R6_2.5.0 TFisher_0.2.0
[135] KEGGdzPathwaysGEO_1.28.0 pillar_1.6.1
[137] htmltools_0.5.1.1 hgu133plus2.db_3.2.3
[139] glue_1.4.2 fastmap_1.1.0
[141] DT_0.17 BiocParallel_1.24.1
[143] codetools_0.2-18 mvtnorm_1.1-1
[145] utf8_1.2.1 lattice_0.20-41
[147] numDeriv_2016.8-1.1 hgu133a.db_3.2.3
[149] curl_4.3 gtools_3.8.2
[151] openssl_1.4.3 survival_3.2-7
[153] rmarkdown_2.6 org.Rn.eg.db_3.12.0
[155] munsell_0.5.0 GenomeInfoDbData_1.2.4
[157] iterators_1.0.13 haven_2.3.1
[159] reshape2_1.4.4 gtable_0.3.0
[161] rbibutils_2.0
sessionInfo()R version 4.0.2 (2020-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS: /config/binaries/R/4.0.2/lib64/R/lib/libRblas.so
LAPACK: /config/binaries/R/4.0.2/lib64/R/lib/libRlapack.so
locale:
[1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_AU.UTF-8 LC_COLLATE=en_AU.UTF-8
[5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8
[7] LC_PAPER=en_AU.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] RColorBrewer_1.1-2 EGSEA_1.18.1
[3] pathview_1.30.1 topGO_2.42.0
[5] SparseM_1.78 GO.db_3.12.1
[7] graph_1.68.0 gage_2.40.1
[9] patchwork_1.1.1 NMF_0.23.0
[11] cluster_2.1.0 rngtools_1.5
[13] pkgmaker_0.32.2 registry_0.5-1
[15] edgeR_3.32.1 limma_3.46.0
[17] EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.14.0
[19] AnnotationFilter_1.14.0 GenomicFeatures_1.42.1
[21] AnnotationDbi_1.52.0 Biobase_2.50.0
[23] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7
[25] IRanges_2.24.1 S4Vectors_0.28.1
[27] BiocGenerics_0.36.1 forcats_0.5.1
[29] stringr_1.4.0 dplyr_1.0.4
[31] purrr_0.3.4 readr_1.4.0
[33] tidyr_1.1.2 tibble_3.1.2
[35] ggplot2_3.3.5 tidyverse_1.3.0
[37] here_1.0.1 workflowr_1.6.2
loaded via a namespace (and not attached):
[1] rappdirs_0.3.3 rtracklayer_1.50.0
[3] Glimma_2.0.0 bit64_4.0.5
[5] knitr_1.31 multcomp_1.4-16
[7] DelayedArray_0.16.3 data.table_1.13.6
[9] hwriter_1.3.2 KEGGREST_1.30.1
[11] RCurl_1.98-1.3 doParallel_1.0.16
[13] generics_0.1.0 metap_1.4
[15] org.Mm.eg.db_3.12.0 TH.data_1.0-10
[17] RSQLite_2.2.5 bit_4.0.4
[19] mutoss_0.1-12 xml2_1.3.2
[21] lubridate_1.7.9.2 httpuv_1.5.5
[23] SummarizedExperiment_1.20.0 assertthat_0.2.1
[25] xfun_0.23 hms_1.0.0
[27] evaluate_0.14 promises_1.2.0.1
[29] fansi_0.5.0 progress_1.2.2
[31] caTools_1.18.1 dbplyr_2.1.0
[33] readxl_1.3.1 Rgraphviz_2.34.0
[35] DBI_1.1.1 geneplotter_1.68.0
[37] tmvnsim_1.0-2 htmlwidgets_1.5.3
[39] ellipsis_0.3.2 backports_1.2.1
[41] annotate_1.68.0 PADOG_1.32.0
[43] gbRd_0.4-11 gridBase_0.4-7
[45] biomaRt_2.46.3 MatrixGenerics_1.2.1
[47] HTMLUtils_0.1.7 vctrs_0.3.8
[49] cachem_1.0.4 withr_2.4.2
[51] globaltest_5.44.0 GenomicAlignments_1.26.0
[53] prettyunits_1.1.1 mnormt_2.0.2
[55] lazyeval_0.2.2 crayon_1.4.1
[57] genefilter_1.72.1 labeling_0.4.2
[59] pkgconfig_2.0.3 nlme_3.1-152
[61] ProtGenerics_1.22.0 GSA_1.03.1
[63] rlang_0.4.11 lifecycle_1.0.0
[65] sandwich_3.0-0 BiocFileCache_1.14.0
[67] mathjaxr_1.2-0 modelr_0.1.8
[69] cellranger_1.1.0 rprojroot_2.0.2
[71] GSVA_1.38.2 matrixStats_0.59.0
[73] Matrix_1.3-2 zoo_1.8-9
[75] reprex_1.0.0 whisker_0.4
[77] png_0.1-7 viridisLite_0.4.0
[79] bitops_1.0-7 KernSmooth_2.23-18
[81] Biostrings_2.58.0 blob_1.2.1
[83] R2HTML_2.3.2 doRNG_1.8.2
[85] scales_1.1.1 memoise_2.0.0.9000
[87] GSEABase_1.52.1 magrittr_2.0.1
[89] plyr_1.8.6 safe_3.30.0
[91] gplots_3.1.1 zlibbioc_1.36.0
[93] compiler_4.0.2 plotrix_3.8-1
[95] KEGGgraph_1.50.0 DESeq2_1.30.1
[97] Rsamtools_2.6.0 cli_3.0.0
[99] XVector_0.30.0 EGSEAdata_1.18.0
[101] MASS_7.3-53.1 tidyselect_1.1.0
[103] stringi_1.5.3 highr_0.8
[105] yaml_2.2.1 askpass_1.1
[107] locfit_1.5-9.4 ggrepel_0.9.1
[109] grid_4.0.2 tools_4.0.2
[111] rstudioapi_0.13 foreach_1.5.1
[113] git2r_0.28.0 farver_2.1.0
[115] digest_0.6.27 Rcpp_1.0.6
[117] broom_0.7.4 later_1.1.0.1
[119] org.Hs.eg.db_3.12.0 httr_1.4.2
[121] Rdpack_2.1 colorspace_2.0-2
[123] rvest_0.3.6 XML_3.99-0.5
[125] fs_1.5.0 splines_4.0.2
[127] statmod_1.4.35 sn_1.6-2
[129] multtest_2.46.0 plotly_4.9.3
[131] xtable_1.8-4 jsonlite_1.7.2
[133] R6_2.5.0 TFisher_0.2.0
[135] KEGGdzPathwaysGEO_1.28.0 pillar_1.6.1
[137] htmltools_0.5.1.1 hgu133plus2.db_3.2.3
[139] glue_1.4.2 fastmap_1.1.0
[141] DT_0.17 BiocParallel_1.24.1
[143] codetools_0.2-18 mvtnorm_1.1-1
[145] utf8_1.2.1 lattice_0.20-41
[147] numDeriv_2016.8-1.1 hgu133a.db_3.2.3
[149] curl_4.3 gtools_3.8.2
[151] openssl_1.4.3 survival_3.2-7
[153] rmarkdown_2.6 org.Rn.eg.db_3.12.0
[155] munsell_0.5.0 GenomeInfoDbData_1.2.4
[157] iterators_1.0.13 haven_2.3.1
[159] reshape2_1.4.4 gtable_0.3.0
[161] rbibutils_2.0