library(readr)
library(limma)
library(BiocParallel)
## Warning: package 'BiocParallel' was built under R version 4.3.1
library(ComplexHeatmap)
## Warning: package 'ComplexHeatmap' was built under R version 4.3.1
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.18.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
library(clusterProfiler)
## Warning: package 'clusterProfiler' was built under R version 4.3.1
##
## Registered S3 methods overwritten by 'treeio':
## method from
## MRCA.phylo tidytree
## MRCA.treedata tidytree
## Nnode.treedata tidytree
## Ntip.treedata tidytree
## ancestor.phylo tidytree
## ancestor.treedata tidytree
## child.phylo tidytree
## child.treedata tidytree
## full_join.phylo tidytree
## full_join.treedata tidytree
## groupClade.phylo tidytree
## groupClade.treedata tidytree
## groupOTU.phylo tidytree
## groupOTU.treedata tidytree
## inner_join.phylo tidytree
## inner_join.treedata tidytree
## is.rooted.treedata tidytree
## nodeid.phylo tidytree
## nodeid.treedata tidytree
## nodelab.phylo tidytree
## nodelab.treedata tidytree
## offspring.phylo tidytree
## offspring.treedata tidytree
## parent.phylo tidytree
## parent.treedata tidytree
## root.treedata tidytree
## rootnode.phylo tidytree
## sibling.phylo tidytree
## clusterProfiler v4.8.3 For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
##
## If you use clusterProfiler in published research, please cite:
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141
##
## Attaching package: 'clusterProfiler'
## The following object is masked from 'package:stats':
##
## filter
library(EnhancedVolcano)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.3.3
## Loading required package: ggrepel
## Warning: package 'ggrepel' was built under R version 4.3.3
library(tidyverse)
## Warning: package 'tidyr' was built under R version 4.3.2
## Warning: package 'purrr' was built under R version 4.3.3
## Warning: package 'lubridate' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ stringr 1.5.1
## ✔ forcats 1.0.0.9000 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks clusterProfiler::filter(), stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ purrr::simplify() masks clusterProfiler::simplify()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(clusterProfiler)
library(org.Hs.eg.db)
## Loading required package: AnnotationDbi
## Warning: package 'AnnotationDbi' was built under R version 4.3.2
## Loading required package: stats4
## Loading required package: BiocGenerics
## Warning: package 'BiocGenerics' was built under R version 4.3.1
##
## Attaching package: 'BiocGenerics'
##
## The following objects are masked from 'package:lubridate':
##
## intersect, setdiff, union
##
## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union
##
## The following object is masked from 'package:limma':
##
## plotMA
##
## 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
##
## Loading required package: Biobase
## Warning: package 'Biobase' was built under R version 4.3.1
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Loading required package: IRanges
## Warning: package 'IRanges' was built under R version 4.3.1
## Loading required package: S4Vectors
## Warning: package 'S4Vectors' was built under R version 4.3.2
##
## Attaching package: 'S4Vectors'
##
## The following objects are masked from 'package:lubridate':
##
## second, second<-
##
## The following objects are masked from 'package:dplyr':
##
## first, rename
##
## The following object is masked from 'package:tidyr':
##
## expand
##
## The following object is masked from 'package:clusterProfiler':
##
## rename
##
## The following object is masked from 'package:utils':
##
## findMatches
##
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
##
##
## Attaching package: 'IRanges'
##
## The following object is masked from 'package:lubridate':
##
## %within%
##
## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice
##
## The following object is masked from 'package:purrr':
##
## reduce
##
## The following object is masked from 'package:clusterProfiler':
##
## slice
##
##
## Attaching package: 'AnnotationDbi'
##
## The following object is masked from 'package:dplyr':
##
## select
##
## The following object is masked from 'package:clusterProfiler':
##
## select
library(dplyr)
library(preprocessCore)
library(M3C)
Colsum_median_Norm_Spleen_Proteomics <- read.csv("~/Desktop/Spleen.Analysis.13Aug/Colsum.median.Norm.Spleen.Proteomics.csv", row.names = 1)
Spleen_Metadata <- read_csv("~/Desktop/Spleen.Analysis.13Aug/Spleen.Metadata.csv")
## Rows: 20 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): Sample, ID, Condition, Group
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dim(Colsum_median_Norm_Spleen_Proteomics)
## [1] 7733 20
Spleen_Metadata <- as.data.frame(Spleen_Metadata)
row.names(Spleen_Metadata) <- colnames(Colsum_median_Norm_Spleen_Proteomics)
metadata <- Spleen_Metadata
##PCA Analysis
library(M3C)
#PCA plot after Normalization
# Define a list of 10 distinct colors for the groups
color_list <- c('darkgreen', 'purple', 'brown', 'orange', 'pink', 'brown', 'grey', 'blue')
# Perform PCA with the specified color vector and other parameters
pca(Colsum_median_Norm_Spleen_Proteomics,
colvec = color_list,
labels = Spleen_Metadata$Group,
controlscale = TRUE,
scale = 3,
printres = TRUE,
printwidth = 25,
legendtitle = "Group")
## ***PCA wrapper function***
## running...
## printing PCA to current directory...
## done.
### case1.Group3.Mutated.Saline.vs.Group1.Control.Saline
#Step 1: Subset metadata
# Subset metadata for CASE 1
case1_metadata.Group3.Group1 <- Spleen_Metadata[Spleen_Metadata$Condition %in% c("Group3.Mutated.Saline", "Group1.Control.Saline"), ]
#Step 2: Subset count matrix based on sample IDs
# Extract sample IDs from case1 metadata
case1_ids <- case1_metadata.Group3.Group1$Sample
# Subset counts using sample IDs
case1_counts.Group3.Group1 <- Colsum_median_Norm_Spleen_Proteomics[, case1_ids]
dim(case1_counts.Group3.Group1)
## [1] 7733 9
colnames(case1_counts.Group3.Group1)==row.names(case1_metadata.Group3.Group1)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
saveRDS(case1_counts.Group3.Group1, file = 'cleanDat.case1_counts.Group3.Group1.Spleen.rds')
saveRDS(case1_metadata.Group3.Group1, file = 'numericMeta.case1_counts.Group3.Group1.Spleen.rds')
#DEP Analysis
cleanDat <- case1_counts.Group3.Group1
numericMeta <- case1_metadata.Group3.Group1
condition <- numericMeta$Condition
Grouping <- numericMeta$Condition
source("/Users/usri/Desktop/Screen/parANOVA.dex.R")
ANOVAout <- parANOVA.dex()
## - parallelThreads variable not set. Running with 2 threads only.
## - Network color assignment vector not supplied or not of length in rows of cleanDat; will not be included in output table and data frame.
## - twoGroupCorrMethod variable not set to a correction method for p.adjust when only 2 groups of samples specified in Grouping. Using Benjamini Hochberg 'BH' FDR.
## - fallbackIfSmallTukeyP variable not set. Using recommended Bonferroni t-test FDR for unreliable Tukey p values <10^-8.5
##
## Attaching package: 'foreach'
##
## The following objects are masked from 'package:purrr':
##
## accumulate, when
##
## ...Tukey p<10^-8.5 Fallback calculations using Bonferroni corrected T test: 0 [0%]
labelTop <- 10
FCmin= 1
flip = -3
signifP=0.05
plotVolc()
## - No comparison p value columns selected in selectComps. Using ALL comparisons.
## - selectComps may not reference valid integer p value column indexes of ANOVAout (or CORout).
## Output will be for all comparisons or correlation(s).
## - useNETcolors not set or not TRUE/FALSE. NETcolors not found in ANOVAout, so 3 color volcano(es) will be drawn.
## - Variable outputfigs not specified. Saving volcano plots to /Users/usri/Desktop/Spleen.Analysis.13Aug/Spleen.Analysis .
##
## ...Thresholds:
## [x] Applying a 100% minimum fold change threshold at + and - x=1 .
## [y] with minimum significance cutoff p < 0.05, equivalent to -log10(p) y=1.3 .
##
## Processing ANOVA column 3 (Group3.Mutated.Saline vs Group1.Control.Saline) for volcano ...
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Generating PDF volcano for ANOVAout column 3 (Group3.Mutated.Saline vs Group1.Control.Saline) ...
## Generating HTML volcano for ANOVAout column 3 (Group3.Mutated.Saline vs Group1.Control.Saline) ...
## Warning in geom2trace.default(dots[[1L]][[3L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomLabelRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
## Warning in geom2trace.default(dots[[1L]][[3L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomLabelRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
## Warning in geom2trace.default(dots[[1L]][[3L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomLabelRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
# Step 1: Subset metadata
case2_metadata.Group4.Group3 <- Spleen_Metadata[Spleen_Metadata$Condition %in% c("Group4.Mutated.LPS", "Group3.Mutated.Saline"), ]
# Step 2: Subset count matrix based on sample IDs
case2_ids <- case2_metadata.Group4.Group3$Sample
case2_counts.Group4.Group3 <- Colsum_median_Norm_Spleen_Proteomics[, case2_ids]
# Check dimensions and ID matching
dim(case2_counts.Group4.Group3)
## [1] 7733 13
# Save outputs
saveRDS(case2_counts.Group4.Group3, file = 'cleanDat.case2_counts.Group4.Group3.Spleen.rds')
saveRDS(case2_metadata.Group4.Group3, file = 'numericMeta.case2_counts.Group4.Group3.Spleen.rds')
#DEP Analysis
cleanDat <- case2_counts.Group4.Group3
numericMeta <- case2_metadata.Group4.Group3
condition <- numericMeta$Condition
Grouping <- numericMeta$Condition
source("/Users/usri/Desktop/Screen/parANOVA.dex.R")
ANOVAout <- parANOVA.dex()
## - parallelThreads variable not set. Running with 2 threads only.
## - Network color assignment vector not supplied or not of length in rows of cleanDat; will not be included in output table and data frame.
## - twoGroupCorrMethod variable not set to a correction method for p.adjust when only 2 groups of samples specified in Grouping. Using Benjamini Hochberg 'BH' FDR.
## - fallbackIfSmallTukeyP variable not set. Using recommended Bonferroni t-test FDR for unreliable Tukey p values <10^-8.5
##
## ...Tukey p<10^-8.5 Fallback calculations using Bonferroni corrected T test: 13 [0.17%]
labelTop <- 15
FCmin= 1
flip = -3
signifP=0.05
plotVolc()
## - No comparison p value columns selected in selectComps. Using ALL comparisons.
## - selectComps may not reference valid integer p value column indexes of ANOVAout (or CORout).
## Output will be for all comparisons or correlation(s).
## - useNETcolors not set or not TRUE/FALSE. NETcolors not found in ANOVAout, so 3 color volcano(es) will be drawn.
## - Variable outputfigs not specified. Saving volcano plots to /Users/usri/Desktop/Spleen.Analysis.13Aug/Spleen.Analysis .
##
## ...Thresholds:
## [x] Applying a 100% minimum fold change threshold at + and - x=1 .
## [y] with minimum significance cutoff p < 0.05, equivalent to -log10(p) y=1.3 .
##
## Processing ANOVA column 3 (Group4.Mutated.LPS vs Group3.Mutated.Saline) for volcano ...
## Generating PDF volcano for ANOVAout column 3 (Group4.Mutated.LPS vs Group3.Mutated.Saline) ...
## Generating HTML volcano for ANOVAout column 3 (Group4.Mutated.LPS vs Group3.Mutated.Saline) ...
## Warning in geom2trace.default(dots[[1L]][[3L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomLabelRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
## Warning in geom2trace.default(dots[[1L]][[3L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomLabelRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
## Warning in geom2trace.default(dots[[1L]][[3L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomLabelRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
# Step 1: Subset metadata
case3_metadata.Group4.Group2 <- Spleen_Metadata[Spleen_Metadata$Condition %in% c("Group4.Mutated.LPS", "Group2.Control.LPS"), ]
# Step 2: Subset count matrix based on sample IDs
case3_ids <- case3_metadata.Group4.Group2$Sample
case3_counts.Group4.Group2 <- Colsum_median_Norm_Spleen_Proteomics[, case3_ids]
# Check dimensions and ID matching
dim(case3_counts.Group4.Group2)
## [1] 7733 10
colnames(case3_counts.Group4.Group2) == row.names(case3_metadata.Group4.Group2)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
saveRDS(case3_counts.Group4.Group2, file = 'cleanDat.case3_counts.Group4.Group2.Spleen.rds')
saveRDS(case3_metadata.Group4.Group2, file = 'numericMeta.case3_counts.Group4.Group2.Spleen.rds')
#DEP Analysis
#DEP Analysis
cleanDat <- case3_counts.Group4.Group2
numericMeta <- case3_metadata.Group4.Group2
condition <- numericMeta$Condition
Grouping <- numericMeta$Condition
source("/Users/usri/Desktop/Screen/parANOVA.dex.R")
ANOVAout <- parANOVA.dex()
## - parallelThreads variable not set. Running with 2 threads only.
## - Network color assignment vector not supplied or not of length in rows of cleanDat; will not be included in output table and data frame.
## - twoGroupCorrMethod variable not set to a correction method for p.adjust when only 2 groups of samples specified in Grouping. Using Benjamini Hochberg 'BH' FDR.
## - fallbackIfSmallTukeyP variable not set. Using recommended Bonferroni t-test FDR for unreliable Tukey p values <10^-8.5
##
## ...Tukey p<10^-8.5 Fallback calculations using Bonferroni corrected T test: 0 [0%]
labelTop <- 15
FCmin= 1
flip = -3
signifP=0.05
plotVolc()
## - No comparison p value columns selected in selectComps. Using ALL comparisons.
## - selectComps may not reference valid integer p value column indexes of ANOVAout (or CORout).
## Output will be for all comparisons or correlation(s).
## - useNETcolors not set or not TRUE/FALSE. NETcolors not found in ANOVAout, so 3 color volcano(es) will be drawn.
## - Variable outputfigs not specified. Saving volcano plots to /Users/usri/Desktop/Spleen.Analysis.13Aug/Spleen.Analysis .
##
## ...Thresholds:
## [x] Applying a 100% minimum fold change threshold at + and - x=1 .
## [y] with minimum significance cutoff p < 0.05, equivalent to -log10(p) y=1.3 .
##
## Processing ANOVA column 3 (Group4.Mutated.LPS vs Group2.Control.LPS) for volcano ...
## Generating PDF volcano for ANOVAout column 3 (Group4.Mutated.LPS vs Group2.Control.LPS) ...
## Generating HTML volcano for ANOVAout column 3 (Group4.Mutated.LPS vs Group2.Control.LPS) ...
## Warning in geom2trace.default(dots[[1L]][[3L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomLabelRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
## Warning in geom2trace.default(dots[[1L]][[3L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomLabelRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
## Warning in geom2trace.default(dots[[1L]][[3L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomLabelRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
