The code presented here supports the Mukherjee et al publication entitled “DAB2IP-Low Luminal A Breast Cancer Patients Exhibit Gene Expression Profiles Overlapping with More Aggressive Cancer Partly Related to IKK/NF-kB Activation” (2023), under review. The steps outlined here recapitulate the analysis performed for preparation of the figures, however, the figures will be subtly different given slight changes in software versions that occurred during manuscript preparation.
The data utilized here were obtained from cBioPortal in February 2019 corresponding to the PanCan Atlas 2018 release and Firehouse Legacy data. The latter contains ER/PR/HER2 status for each patient sample.
Note that the GDC portal no longer includes HTseq gene quantifications but rather includes augmented_star_gene_counts
, which should be similar, and RNA-seq data also now utilize RSEM whereas the data used here do not.
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.0 ✔ readr 2.1.2
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.0 ✔ tibble 3.1.8
## ✔ lubridate 1.9.3 ✔ tidyr 1.2.0
## ✔ purrr 0.3.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ComplexHeatmap)
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.10.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:
## 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(circlize)
## ========================================
## circlize version 0.4.14
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
##
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
## in R. Bioinformatics 2014.
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(circlize))
## ========================================
library(tidyHeatmap)
## ========================================
## tidyHeatmap version 1.8.1
## If you use tidyHeatmap in published research, please cite:
## 1) Mangiola et al. tidyHeatmap: an R package for modular heatmap production
## based on tidy principles. JOSS 2020.
## 2) Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
## This message can be suppressed by:
## suppressPackageStartupMessages(library(tidyHeatmap))
## ========================================
##
##
## Attaching package: 'tidyHeatmap'
##
## The following object is masked from 'package:stats':
##
## heatmap
library(readxl)
library(DESeq2)
## Loading required package: S4Vectors
## Loading required package: stats4
## Loading required package: BiocGenerics
##
## 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 objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
##
## The following objects are masked from 'package:base':
##
## anyDuplicated, 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 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 objects are masked from 'package:base':
##
## expand.grid, I, unname
##
## Loading required package: IRanges
##
## 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
##
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
##
## Attaching package: 'matrixStats'
##
## The following object is masked from 'package:dplyr':
##
## count
##
##
## 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(ashr)
library(paletteer)
library(patchwork)
EPHstatus <- read_tsv("brca_tcga_clinical_data_ER_PR_HER2.tsv") %>%
dplyr::select(c(`Sample ID`,`ER Status By IHC`,`IHC-HER2`,`PR status by ihc`)) %>%
mutate(
TNBC = case_when(
`ER Status By IHC` == "Negative" & `IHC-HER2` == "Negative" & `PR status by ihc` == "Negative" ~ "TNBC",
T ~ "nonTNBC"
)
) %>%
dplyr::rename("SampleID" = `Sample ID`,
"ERstatus" = `ER Status By IHC`,
"HER2status" = `IHC-HER2`,
"PRstatus" = `PR status by ihc`) %>%
mutate(SampleID = paste0(SampleID,"A"))
## New names:
## Rows: 1108 Columns: 139
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "\t" chr
## (88): Study ID, Patient ID, Sample ID, American Joint Committee on Cance... dbl
## (20): Diagnosis Age, Days to Sample Collection., Last Alive Less Initial... lgl
## (31): Neoplasm American Joint Committee on Cancer Clinical Group Stage, ...
## ℹ 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.
## • `First Pathologic Diagnosis Biospecimen Acquisition Method Type` -> `First
## Pathologic Diagnosis Biospecimen Acquisition Method Type...68`
## • `First Pathologic Diagnosis Biospecimen Acquisition Method Type` -> `First
## Pathologic Diagnosis Biospecimen Acquisition Method Type...70`
Subtype <- read_tsv("data_clinical_patient.txt",skip=4) %>%
dplyr::rename("SampleID" = PATIENT_ID) %>%
dplyr::select(c(SampleID,SUBTYPE)) %>%
mutate(SampleID = paste0(SampleID,"-01A"))
## Rows: 1084 Columns: 38
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (29): PATIENT_ID, SUBTYPE, CANCER_TYPE_ACRONYM, OTHER_PATIENT_ID, SEX, A...
## dbl (8): AGE, DAYS_LAST_FOLLOWUP, DAYS_TO_BIRTH, DAYS_TO_INITIAL_PATHOLOGIC...
## lgl (1): WEIGHT
##
## ℹ 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.
TCGAexpr <- read_tsv("data_RNA_Seq_v2_mRNA_median_all_sample_Zscores.txt") %>%
dplyr::filter(!is.na(Hugo_Symbol)) %>%
dplyr::select(!starts_with("Entrez")) %>%
pivot_longer(cols=!starts_with("Hugo"),names_to="SampleID",values_to="Expression") %>%
dplyr::rename("Gene" = Hugo_Symbol) %>%
mutate(SampleID = paste0(SampleID,"A")) %>%
dplyr::select(c(SampleID,Gene,Expression))
## Rows: 20531 Columns: 1084
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): Hugo_Symbol
## dbl (1083): Entrez_Gene_Id, TCGA-3C-AAAU-01, TCGA-3C-AALI-01, TCGA-3C-AALJ-0...
##
## ℹ 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.
DAB2IP <- TCGAexpr %>%
dplyr::filter(Gene=="DAB2IP") %>%
mutate(quartile = ntile(Expression, 4)) %>%
mutate(
DAB2IP = case_when(
quartile == 1 ~ "DAB2IP_low",
quartile == 4 ~ "DAB2IP_high",
T ~ "DAB2IP_mid"
)
) %>%
dplyr::select(SampleID,DAB2IP,Expression)
combined <- inner_join(inner_join(DAB2IP,EPHstatus),Subtype)
## Joining with `by = join_by(SampleID)`
## Joining with `by = join_by(SampleID)`
combined <- combined %>%
drop_na(SUBTYPE)
# Get min cutoffs for each DAB2IP expression quartile
DAB2IP.cuts <- TCGAexpr %>%
filter(Gene=="DAB2IP") %>%
mutate(quartile = ntile(Expression, 4)) %>%
group_by(quartile) %>%
summarize(max = max(Expression)) %>%
pull(max)
combined %>%
select(SampleID,Expression,SUBTYPE) %>%
group_by(SUBTYPE) %>%
ggplot(aes(x=SUBTYPE,y=Expression,group=SUBTYPE,fill=SUBTYPE)) +
geom_boxplot(outlier.shape=NA) +
annotate("rect", xmin=5.5, xmax=Inf, ymin=-Inf, ymax=DAB2IP.cuts[1], alpha=0.4, fill="gray90") +
annotate("rect", xmin=5.5, xmax=Inf, ymin=DAB2IP.cuts[1], ymax=DAB2IP.cuts[2], alpha=0.4, fill = "gray50") +
annotate("rect", xmin=5.5, xmax=Inf, ymin=DAB2IP.cuts[2], ymax=DAB2IP.cuts[3], alpha=0.4, fill = "gray20") +
annotate("rect", xmin=5.5, xmax=Inf, ymin=DAB2IP.cuts[3], ymax=Inf, alpha=0.4, fill = "black") +
geom_jitter(width=0.25) +
ylab("DAB2IP expression") +
xlab("Subtype") +
theme_classic() +
geom_text(label="1",x=5.55,y=-2.5) +
geom_text(label="2",x=5.55,y=-0.3) +
geom_text(label="3",x=5.55,y=0.4) +
geom_text(label="4",x=5.55,y=2)
## png
## 2
RORP <- read_tsv("BRCA.1218_pam50scores.FINAL.txt") %>%
dplyr::rename("SampleID" = bcr_patient_barcode, "RORPscores" = `ROR-P (Subtype + Proliferation)`) %>%
dplyr::select(SampleID,RORPscores) %>%
mutate(SampleID = paste0(SampleID,"-01A")) %>%
group_by(SampleID) %>%
summarize(RORPavg = mean(RORPscores))
## New names:
## Rows: 1222 Columns: 20
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "\t" chr
## (9): bcr_patient_barcode, Barcode, Run.Lane.barcode, Call, ROR-S Group ... dbl
## (11): Basal, Her2...5, LumA, LumB, Normal, Confidence, ROR-S (Subtype On...
## ℹ 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.
## • `Her2` -> `Her2...5`
## • `Her2` -> `Her2...17`
inner_join(TCGAexpr,inner_join(RORP,combined,by="SampleID"),by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE,RORPavg) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(Gene=="DAB2IP") %>%
dplyr::filter(DAB2IP != "DAB2IP_mid") %>%
ggplot(aes(x=fct_relevel(DAB2IP,c("DAB2IP_low","DAB2IP_high")),y=RORPavg,fill=SUBTYPE)) +
geom_boxplot(outlier.shape=NA) +
geom_jitter(width=0.25) +
xlab("DAB2IP expression") +
facet_wrap(~SUBTYPE,ncol=5) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
## png
## 2
rorp.d.low.lumA <- inner_join(TCGAexpr,inner_join(RORP,combined,by="SampleID"),by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE,RORPavg) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(Gene=="DAB2IP") %>%
dplyr::filter(SUBTYPE=="BRCA_LumA") %>%
dplyr::filter(DAB2IP=="DAB2IP_low") %>%
pull(RORPavg)
rorp.d.high.lumA <- inner_join(TCGAexpr,inner_join(RORP,combined,by="SampleID"),by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE,RORPavg) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(Gene=="DAB2IP") %>%
dplyr::filter(SUBTYPE=="BRCA_LumA") %>%
dplyr::filter(DAB2IP=="DAB2IP_high") %>%
pull(RORPavg)
t.test(rorp.d.high.lumA,rorp.d.low.lumA)
##
## Welch Two Sample t-test
##
## data: rorp.d.high.lumA and rorp.d.low.lumA
## t = -6.6269, df = 201.56, p-value = 3.064e-10
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -17.189772 -9.306095
## sample estimates:
## mean of x mean of y
## 7.464012 20.711945
rorp.d.low.lumB <- inner_join(TCGAexpr,inner_join(RORP,combined,by="SampleID"),by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE,RORPavg) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(Gene=="DAB2IP") %>%
dplyr::filter(SUBTYPE=="BRCA_LumB") %>%
dplyr::filter(DAB2IP=="DAB2IP_low") %>%
pull(RORPavg)
rorp.d.high.lumB <- inner_join(TCGAexpr,inner_join(RORP,combined,by="SampleID"),by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE,RORPavg) %>%
dplyr:: rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(Gene=="DAB2IP") %>%
dplyr::filter(SUBTYPE=="BRCA_LumB") %>%
dplyr::filter(DAB2IP=="DAB2IP_high") %>%
pull(RORPavg)
t.test(rorp.d.high.lumB,rorp.d.low.lumB)
##
## Welch Two Sample t-test
##
## data: rorp.d.high.lumB and rorp.d.low.lumB
## t = 0.96542, df = 23.001, p-value = 0.3444
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.704416 10.187711
## sample estimates:
## mean of x mean of y
## 49.27263 46.03098
rorp.d.low.Basal <- inner_join(TCGAexpr,inner_join(RORP,combined,by="SampleID"),by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE,RORPavg) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(Gene=="DAB2IP") %>%
dplyr::filter(SUBTYPE=="BRCA_Basal") %>%
dplyr::filter(DAB2IP=="DAB2IP_low") %>%
pull(RORPavg)
rorp.d.high.Basal <- inner_join(TCGAexpr,inner_join(RORP,combined,by="SampleID"),by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE,RORPavg) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(Gene=="DAB2IP") %>%
dplyr::filter(SUBTYPE=="BRCA_Basal") %>%
dplyr::filter(DAB2IP=="DAB2IP_high") %>%
pull(RORPavg)
t.test(rorp.d.high.Basal,rorp.d.low.Basal)
##
## Welch Two Sample t-test
##
## data: rorp.d.high.Basal and rorp.d.low.Basal
## t = -2.0816, df = 59.824, p-value = 0.04167
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -10.8012654 -0.2147586
## sample estimates:
## mean of x mean of y
## 53.06848 58.57649
rorp.d.low.Her2 <- inner_join(TCGAexpr,inner_join(RORP,combined,by="SampleID"),by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE,RORPavg) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(Gene=="DAB2IP") %>%
dplyr::filter(SUBTYPE=="BRCA_Her2") %>%
dplyr::filter(DAB2IP=="DAB2IP_low") %>%
pull(RORPavg)
rorp.d.high.Her2 <- inner_join(TCGAexpr,inner_join(RORP,combined,by="SampleID"),by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE,RORPavg) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(Gene=="DAB2IP") %>%
dplyr::filter(SUBTYPE=="BRCA_Her2") %>%
dplyr::filter(DAB2IP=="DAB2IP_high") %>%
pull(RORPavg)
t.test(rorp.d.high.Her2,rorp.d.low.Her2)
##
## Welch Two Sample t-test
##
## data: rorp.d.high.Her2 and rorp.d.low.Her2
## t = -0.54965, df = 28.418, p-value = 0.5869
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -14.479670 8.349809
## sample estimates:
## mean of x mean of y
## 43.70862 46.77355
tumorsDE <- combined %>%
dplyr::select(SampleID,DAB2IP,ERstatus,SUBTYPE) %>%
dplyr::filter(ERstatus=="Positive") %>%
dplyr::filter(SUBTYPE=="BRCA_LumA") %>%
dplyr::filter(DAB2IP!="DAB2IP_mid")
write_tsv(tumorsDE,"TCGA_BRCA_tumors_DAB2IP_low_high_ERpositiveOnly_LumA_010924.tsv")
tumorsDE
## # A tibble: 206 × 4
## SampleID DAB2IP ERstatus SUBTYPE
## <chr> <chr> <chr> <chr>
## 1 TCGA-3C-AAAU-01A DAB2IP_high Positive BRCA_LumA
## 2 TCGA-4H-AAAK-01A DAB2IP_high Positive BRCA_LumA
## 3 TCGA-A1-A0SE-01A DAB2IP_high Positive BRCA_LumA
## 4 TCGA-A1-A0SQ-01A DAB2IP_high Positive BRCA_LumA
## 5 TCGA-A2-A04V-01A DAB2IP_high Positive BRCA_LumA
## 6 TCGA-A2-A0CK-01A DAB2IP_high Positive BRCA_LumA
## 7 TCGA-A2-A0CP-01A DAB2IP_high Positive BRCA_LumA
## 8 TCGA-A2-A0CR-01A DAB2IP_high Positive BRCA_LumA
## 9 TCGA-A2-A0CT-01A DAB2IP_low Positive BRCA_LumA
## 10 TCGA-A2-A0CV-01A DAB2IP_high Positive BRCA_LumA
## # … with 196 more rows
# Collapse all htseq output files from isoform->gene
tx2gene <- read_tsv("rawCounts_DE/gencode.gene.info.v22.lookup.tsv",col_names=c("tx","gene"))
sampleFiles <- grep("htseq.counts$",list.files("rawCounts_DE"),value=TRUE)
for (i in 1:length(sampleFiles)) {
file <- sampleFiles[i]
tsv <- read_tsv(paste0("rawCounts_DE/",file),col_names=c("tx","count"))
joined <- inner_join(tsv,tx2gene) %>%
group_by(gene) %>%
summarize(sum = sum(count))
newFile <- paste0("rawCounts_DE/",file,".030323.collapsed")
write_tsv(joined,newFile,col_names=F)
}
# Join sample IDs with filenames of new collapsed counts files and compile DESeq2 metadata
sampleFiles <- grep("counts.030323.collapsed",list.files("rawCounts_DE"),value=TRUE)
sampleInfo <- read_tsv("TCGA_BRCA_tumors_DAB2IP_low_high_ERpositiveOnly_LumA_010924.tsv")
## Rows: 206 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): SampleID, DAB2IP, ERstatus, SUBTYPE
##
## ℹ 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.
countFiles <- read_tsv("rawCounts_DE/gdc_sample_sheet.2020-11-25.tsv") %>%
dplyr::select(`Sample ID`,`File Name`) %>%
dplyr::rename("SampleID" = `Sample ID`,"Filename" = `File Name`)
## Rows: 1222 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (8): File ID, File Name, Data Category, Data Type, Project ID, Case ID, ...
##
## ℹ 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.
coldata <- inner_join(sampleInfo,countFiles) %>%
dplyr::select(SampleID,Filename,DAB2IP) %>%
mutate(Filename = str_replace_all(Filename,"counts.gz","counts.030323.collapsed")) %>%
dplyr::filter(Filename != "0c7d119e-22c0-47ed-b736-3404ab2c7065.htseq.counts.030323.collapsed") %>%
dplyr::filter(Filename != "6aa69f7f-28c6-4943-b907-62d8085e9c6c.htseq.counts.030323.collapsed")
## Joining with `by = join_by(SampleID)`
## Warning in inner_join(sampleInfo, countFiles): Each row in `x` is expected to match at most 1 row in `y`.
## ℹ Row 27 of `x` matches multiple rows.
## ℹ If multiple matches are expected, set `multiple = "all"` to silence this
## warning.
# These two samples have two counts files, removing the second instance of both
coldata <- as.data.frame(coldata)
coldata$DAB2IP <- factor(coldata$DAB2IP)
# Set up DESeq2 object
dds <- DESeqDataSetFromHTSeqCount(sampleTable = coldata, directory = "rawCounts_DE", design= ~DAB2IP)
# Run DE test
dds <- DESeq(dds,betaPrior=FALSE)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 4830 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
# Use shrunken log-fold-changes with `ashr`
res <- lfcShrink(dds, contrast=c("DAB2IP","DAB2IP_high","DAB2IP_low"), type="ashr")
## using 'ashr' for LFC shrinkage. If used in published research, please cite:
## Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
## https://doi.org/10.1093/biostatistics/kxw041
# VST normalization for downstream plotting
vsd <- vst(dds)
DESeq2::plotPCA(vsd, intgroup="DAB2IP")
# Prepare output files
res <- res[order(res$padj), ]
resdata <- merge(as.data.frame(res,row.names=rownames(res)), as.data.frame(assay(vsd)), by="row.names", sort=FALSE)
rownames(resdata)=resdata$Row.names
resdata=resdata[,-1]
# Plot histogram of significant adjusted p-values
hist(res$padj, breaks=50, col="grey")
# Filter results
res_tbl <- rownames_to_column(as.data.frame(res), var = "Gene") %>%
as_tibble()
res_filt <- res_tbl %>%
dplyr::select(Gene,log2FoldChange,padj) %>%
dplyr::filter(padj < 1e-5)
# Write output
write_tsv(res_filt,"TCGA_BRCA_tumors_DAB2IP_low_high_ERpositiveOnly_LumA_DESeq2_results_adjP1e5_010924.txt")
DEG <- read_tsv("TCGA_BRCA_tumors_DAB2IP_low_high_ERpositiveOnly_LumA_DESeq2_results_adjP1e5_010924.txt")
## Rows: 2073 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): Gene
## dbl (2): log2FoldChange, padj
##
## ℹ 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.
DEG_any <- DEG %>%
pull(Gene)
# Subset TCGA gene expression data to just the differentially-expressed genes
### Include SFRS1 manually, since we know this to be differentially-expressed but uses an Entrez synonymous gene symbol for SRSF1
TCGAexpr_DEG <- TCGAexpr %>%
dplyr::filter(Gene %in% DEG_any | Gene=="SFRS1")
DAB2IP_order <- inner_join(TCGAexpr_DEG,combined,by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(!is.na(SUBTYPE)) %>%
dplyr::filter(ERstatus=="Positive" | ERstatus=="Negative") %>%
dplyr::filter(DAB2IP=="DAB2IP_high" | DAB2IP=="DAB2IP_low") %>%
# dplyr::filter(!Gene %in% c("SNORD104","SNORD11","SNORD6","SNORD93")) %>% # Remove these due to SD=NA
group_by(SUBTYPE) %>%
dplyr::arrange(DAB2IP_expression,.by_group = T) %>%
dplyr::select(SampleID,DAB2IP_expression) %>%
unique() %>%
pull(SampleID)
## Adding missing grouping variables: `SUBTYPE`
# Now plot heatmap of DAB2IP DEGs for all DAB2IP high/low tumors
inner_join(TCGAexpr_DEG,combined,by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(DAB2IP=="DAB2IP_high" | DAB2IP=="DAB2IP_low") %>%
dplyr::filter(!is.na(SUBTYPE)) %>%
dplyr::filter(ERstatus=="Positive" | ERstatus=="Negative") %>%
# dplyr::filter(!Gene %in% c("SNORD104","SNORD11","SNORD6","SNORD93")) %>% # Remove these due to SD=NA
group_by(SUBTYPE) %>%
tidyHeatmap::heatmap(Gene,
SampleID,
Expression,
scale="none",
palette_value=colorRamp2(c(-3, 0, 3), c("blue", "white", "red")),
cluster_columns=F,
column_order=DAB2IP_order,
show_row_names = FALSE,
show_column_names = FALSE) %>%
add_tile(SUBTYPE) %>%
add_tile(DAB2IP)
## tidyHeatmap says: (once per session) from release 1.7.0 the scaling is set to "none" by default. Please use scale = "row", "column" or "both" to apply scaling
## png
## 2
# Format sample metadata
coldata <- inner_join(TCGAexpr_DEG,combined,by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(DAB2IP=="DAB2IP_high" | DAB2IP=="DAB2IP_low") %>%
dplyr::filter(!is.na(SUBTYPE)) %>%
dplyr::filter(ERstatus=="Positive" | ERstatus=="Negative") %>%
dplyr::filter(Gene=="DAB2IP") %>%
dplyr::select(-Gene,-DAB2IP_expression,-Expression) %>%
mutate(SUBTYPE = str_replace_all(SUBTYPE,"BRCA_","")) %>%
unite(Group,c(SUBTYPE,DAB2IP),sep="_")
coldata.df <- as.data.frame(column_to_rownames(coldata,var="SampleID"))
# Format expression data matrix
degmat <- inner_join(TCGAexpr_DEG,combined,by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(DAB2IP=="DAB2IP_high" | DAB2IP=="DAB2IP_low") %>%
dplyr::filter(!is.na(SUBTYPE)) %>%
dplyr::filter(ERstatus=="Positive" | ERstatus=="Negative") %>%
# dplyr::filter(!Gene %in% c("SNORD104","SNORD11","SNORD6","SNORD93")) %>% # Remove these due to SD=NA
dplyr::select(-DAB2IP,-DAB2IP_expression,-ERstatus,-SUBTYPE) %>%
pivot_wider(names_from=SampleID,values_from=Expression,values_fn = median)
degmat.df <- as.data.frame(column_to_rownames(degmat,var="Gene"))
# Compute Pearson correlations, Luminal A DAB2IP-low vs all
cors <- cor(degmat.df[,coldata$SampleID[coldata$Group=="LumA_DAB2IP_low"]],degmat.df)
# Set up sample groupings and colors
cols <- c(paletteer_d("nationalparkcolors::Redwoods")[1], # Light green
paletteer_d("nationalparkcolors::Everglades")[5], # Dark green
paletteer_d("nationalparkcolors::Everglades")[1], # Light blue
paletteer_d("nationalparkcolors::Everglades")[2], # Dark blue
paletteer_d("nationalparkcolors::Redwoods")[3], # Light yellow
paletteer_d("nationalparkcolors::BryceCanyon")[1], # Dark yellow
paletteer_d("nationalparkcolors::Badlands")[3], # Light red
paletteer_d("nationalparkcolors::DeathValley")[1], # Dark red
paletteer_d("nationalparkcolors::Acadia")[6] # Dark gray
)
names(cols) <- levels(as.factor(coldata$Group))
ha1 <- HeatmapAnnotation(Group = as.factor(coldata$Group), col = list(Group = cols))
ha2 <- rowAnnotation(Group = as.factor(coldata$Group[coldata$Group %in% "LumA_DAB2IP_low"]), col = list(Group = cols), show_legend = F)
# Plot such that DAB2IP_low_LumA (rows) plotted in same order as the clustered columns
# Plot once to get column orders, then plot again with those specified to order the rows
h1 <- ComplexHeatmap::Heatmap(cors,
show_row_names = F,
show_column_names = F,
name = "Expression correlation",
clustering_distance_rows=function(x) as.dist(1-cor(t(x))),
clustering_method_rows = "complete",
clustering_distance_columns=function(x) as.dist(1-cor(t(x))),
clustering_method_columns = "complete",
column_title = NULL,
row_title = NULL,
top_annotation = ha1,
left_annotation = ha2,
cluster_column_slices = F,
column_split = as.factor(coldata$Group),
border = T
)
# Draw unordered version to get order
ht <- draw(h1)
orders = column_order(ht)
row_order <- orders$LumA_DAB2IP_low
# Draw in sorted order
ComplexHeatmap::Heatmap(cors[colnames(cors)[row_order],],
show_row_names = F,
show_column_names = F,
cluster_rows = F,
name = "Expression correlation",
clustering_distance_columns=function(x) as.dist(1-cor(t(x))),
clustering_method_columns = "complete",
column_title = NULL,
row_title = NULL,
top_annotation = ha1,
left_annotation = ha2,
cluster_column_slices = F,
column_split = as.factor(coldata$Group),
border = T
)
## png
## 2
PDscores <- read_tsv("TCGA_Dscore_Pscore.txt")
## Rows: 1198 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): SampleID
## dbl (2): D_Score, P_score_11genesMean
##
## ℹ 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.
left_join(combined,PDscores) %>%
dplyr::select(SampleID,ERstatus,SUBTYPE,D_Score,P_score_11genesMean,Expression,DAB2IP) %>%
dplyr::filter(ERstatus=="Positive") %>%
drop_na(SUBTYPE) %>%
drop_na(P_score_11genesMean) %>%
dplyr::filter(DAB2IP!="DAB2IP_mid") %>%
mutate(DAB2IP = fct_relevel(DAB2IP,c("DAB2IP_low","DAB2IP_high"))) %>%
ggplot(aes(x=DAB2IP,y=P_score_11genesMean,color=SUBTYPE)) +
geom_boxplot() +
geom_jitter() +
facet_wrap(~SUBTYPE,ncol=5)
## Joining with `by = join_by(SampleID)`
## Joining with `by = join_by(SampleID)`
## png
## 2
# LumA
Pscore_LumA_low <- left_join(combined,PDscores) %>%
dplyr::select(SampleID,ERstatus,SUBTYPE,D_Score,P_score_11genesMean,Expression,DAB2IP) %>%
dplyr::filter(ERstatus=="Positive") %>%
drop_na(SUBTYPE) %>%
dplyr::filter(DAB2IP=="DAB2IP_low" & SUBTYPE=="BRCA_LumA") %>%
pull(P_score_11genesMean)
## Joining with `by = join_by(SampleID)`
Pscore_LumA_high <- left_join(combined,PDscores) %>%
dplyr::select(SampleID,ERstatus,SUBTYPE,D_Score,P_score_11genesMean,Expression,DAB2IP) %>%
dplyr::filter(ERstatus=="Positive") %>%
drop_na(SUBTYPE) %>%
dplyr::filter(DAB2IP=="DAB2IP_high" & SUBTYPE=="BRCA_LumA") %>%
pull(P_score_11genesMean)
## Joining with `by = join_by(SampleID)`
t.test(Pscore_LumA_low,Pscore_LumA_high)
##
## Welch Two Sample t-test
##
## data: Pscore_LumA_low and Pscore_LumA_high
## t = 7.1667, df = 194.46, p-value = 1.555e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.5406667 0.9512256
## sample estimates:
## mean of x mean of y
## -0.267900 -1.013846
# LumB
Pscore_LumB_low <- left_join(combined,PDscores) %>%
dplyr::select(SampleID,ERstatus,SUBTYPE,D_Score,P_score_11genesMean,Expression,DAB2IP) %>%
dplyr::filter(ERstatus=="Positive") %>%
drop_na(SUBTYPE) %>%
dplyr::filter(DAB2IP=="DAB2IP_low" & SUBTYPE=="BRCA_LumB") %>%
pull(P_score_11genesMean)
## Joining with `by = join_by(SampleID)`
Pscore_LumB_high <- left_join(combined,PDscores) %>%
dplyr::select(SampleID,ERstatus,SUBTYPE,D_Score,P_score_11genesMean,Expression,DAB2IP) %>%
dplyr::filter(ERstatus=="Positive") %>%
drop_na(SUBTYPE) %>%
dplyr::filter(DAB2IP=="DAB2IP_high" & SUBTYPE=="BRCA_LumB") %>%
pull(P_score_11genesMean)
## Joining with `by = join_by(SampleID)`
t.test(Pscore_LumB_low,Pscore_LumB_high)
##
## Welch Two Sample t-test
##
## data: Pscore_LumB_low and Pscore_LumB_high
## t = 0.017274, df = 15.495, p-value = 0.9864
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2266552 0.2303695
## sample estimates:
## mean of x mean of y
## 0.7344341 0.7325769
# Basal
Pscore_Basal_low <- left_join(combined,PDscores) %>%
dplyr::select(SampleID,ERstatus,SUBTYPE,D_Score,P_score_11genesMean,Expression,DAB2IP) %>%
dplyr::filter(ERstatus=="Positive") %>%
drop_na(SUBTYPE) %>%
dplyr::filter(DAB2IP=="DAB2IP_low" & SUBTYPE=="BRCA_Basal") %>%
pull(P_score_11genesMean)
## Joining with `by = join_by(SampleID)`
Pscore_Basal_high <- left_join(combined,PDscores) %>%
dplyr::select(SampleID,ERstatus,SUBTYPE,D_Score,P_score_11genesMean,Expression,DAB2IP) %>%
dplyr::filter(ERstatus=="Positive") %>%
drop_na(SUBTYPE) %>%
dplyr::filter(DAB2IP=="DAB2IP_high" & SUBTYPE=="BRCA_Basal") %>%
pull(P_score_11genesMean)
## Joining with `by = join_by(SampleID)`
t.test(Pscore_Basal_low,Pscore_Basal_high)
##
## Welch Two Sample t-test
##
## data: Pscore_Basal_low and Pscore_Basal_high
## t = 0.27631, df = 9, p-value = 0.7885
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1950053 0.2492719
## sample estimates:
## mean of x mean of y
## 0.6153000 0.5881667
# Her2
Pscore_Her2_low <- left_join(combined,PDscores) %>%
dplyr::select(SampleID,ERstatus,SUBTYPE,D_Score,P_score_11genesMean,Expression,DAB2IP) %>%
dplyr::filter(ERstatus=="Positive") %>%
drop_na(SUBTYPE) %>%
dplyr::filter(DAB2IP=="DAB2IP_low" & SUBTYPE=="BRCA_Her2") %>%
pull(P_score_11genesMean)
## Joining with `by = join_by(SampleID)`
Pscore_Her2_high <- left_join(combined,PDscores) %>%
dplyr::select(SampleID,ERstatus,SUBTYPE,D_Score,P_score_11genesMean,Expression,DAB2IP) %>%
dplyr::filter(ERstatus=="Positive") %>%
drop_na(SUBTYPE) %>%
dplyr::filter(DAB2IP=="DAB2IP_high" & SUBTYPE=="BRCA_Her2") %>%
pull(P_score_11genesMean)
## Joining with `by = join_by(SampleID)`
t.test(Pscore_Her2_low,Pscore_Her2_high)
##
## Welch Two Sample t-test
##
## data: Pscore_Her2_low and Pscore_Her2_high
## t = -0.27922, df = 7.811, p-value = 0.7873
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.4086571 0.3207142
## sample estimates:
## mean of x mean of y
## 0.3629286 0.4069000
These putative targets were identified using ChIP-Atlas. Signal within 5kb of the transcription start site (TSS) of target genes was averaged for HMEC, MCF7, and MDA-MB-231 cells (n=20 studies total), and filtered to retain promoters with signal greater than or equal to 10
rela_targets <- read_tsv("RELA.5.breast.gt10.tsv",col_names="Gene") %>%
pull(Gene)
## Rows: 2871 Columns: 1
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): Gene
##
## ℹ 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.
inner_join(TCGAexpr_DEG,combined,by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(DAB2IP=="DAB2IP_high" | DAB2IP=="DAB2IP_low") %>%
dplyr::filter(!is.na(SUBTYPE)) %>%
dplyr::filter(ERstatus=="Positive" | ERstatus=="Negative") %>%
dplyr::filter(Gene %in% rela_targets) %>%
# dplyr::filter(!Gene %in% c("SNORD104","SNORD11","SNORD6","SNORD93")) %>% # Remove these due to SD=NA
group_by(SUBTYPE) %>%
tidyHeatmap::heatmap(Gene,
SampleID,
Expression,
scale="none",
palette_value=colorRamp2(c(-3, 0, 3), c("blue", "white", "red")),
cluster_columns=F,
column_order=DAB2IP_order,
show_row_names = FALSE,
show_column_names = FALSE) %>%
add_tile(SUBTYPE) %>%
add_tile(DAB2IP)
## png
## 2
goi <- c("BIRC5", "CDK5", "TBCA", "SFRS1") # Note SFRS1 = SRSF1 as per Entrez gene ID
inner_join(TCGAexpr_DEG,combined,by="SampleID") %>%
dplyr::filter(Gene %in% goi) %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(SUBTYPE != "BRCA_Normal") %>%
dplyr::filter(DAB2IP != "DAB2IP_mid") %>%
ggplot(aes(x=fct_relevel(DAB2IP,c("DAB2IP_low","DAB2IP_high")),y=Expression,fill=SUBTYPE)) +
geom_boxplot(outlier.shape=NA) +
geom_jitter(width=0.25) +
xlab("Expression") +
facet_grid(Gene~SUBTYPE) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
## png
## 2
kegg.splice.genes <- c("SNRPE","SFRS1","SNRPG","SFRS2","SNRPD1","SNRPF","BUD31","MAGOHB","HNRNPU","HNRNPC","SNRPB2","ALYREF","SNRPC","SF3B6","SNRNP27","SNRPA1","USP39","PPIL1","LSM3","SNRPB")
inner_join(TCGAexpr_DEG,combined,by="SampleID") %>%
dplyr::select(SampleID,Gene,Expression.x,DAB2IP,Expression.y,ERstatus,SUBTYPE) %>%
dplyr::filter(Gene %in% kegg.splice.genes) %>%
dplyr::rename("Expression" = Expression.x,"DAB2IP_expression" = Expression.y) %>%
dplyr::filter(DAB2IP=="DAB2IP_high" | DAB2IP=="DAB2IP_low") %>%
dplyr::filter(!is.na(SUBTYPE)) %>%
dplyr::filter(ERstatus=="Positive" | ERstatus=="Negative") %>%
group_by(SUBTYPE) %>%
tidyHeatmap::heatmap(Gene,
SampleID,
Expression,
scale="none",
palette_value=colorRamp2(c(-3, 0, 3), c("blue", "white", "red")),
cluster_columns=F,
column_order=DAB2IP_order,
show_row_names = T,
show_column_names = FALSE) %>%
add_tile(SUBTYPE) %>%
add_tile(DAB2IP)
## png
## 2
save.image("TCGA-BRCA_DAB2IP_010924.RData")
sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux 8.8 (Ootpa)
##
## Matrix products: default
## BLAS/LAPACK: /nas/longleaf/rhel8/apps/r/4.1.0/lib/libopenblas_haswellp-r0.3.5.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 grid stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] patchwork_1.1.3 paletteer_1.4.0
## [3] ashr_2.2-54 DESeq2_1.34.0
## [5] SummarizedExperiment_1.24.0 Biobase_2.54.0
## [7] MatrixGenerics_1.6.0 matrixStats_0.62.0
## [9] GenomicRanges_1.46.1 GenomeInfoDb_1.35.16
## [11] IRanges_2.28.0 S4Vectors_0.32.4
## [13] BiocGenerics_0.40.0 readxl_1.4.3
## [15] tidyHeatmap_1.8.1 circlize_0.4.14
## [17] ComplexHeatmap_2.10.0 lubridate_1.9.3
## [19] forcats_1.0.0 stringr_1.5.1
## [21] dplyr_1.1.0 purrr_0.3.5
## [23] readr_2.1.2 tidyr_1.2.0
## [25] tibble_3.1.8 ggplot2_3.4.0
## [27] tidyverse_2.0.0 rmarkdown_2.12
## [29] knitr_1.45
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-3 rjson_0.2.21 ellipsis_0.3.2
## [4] XVector_0.34.0 GlobalOptions_0.1.2 clue_0.3-60
## [7] farver_2.1.1 bit64_4.0.5 AnnotationDbi_1.56.2
## [10] fansi_1.0.3 codetools_0.2-19 splines_4.1.0
## [13] doParallel_1.0.17 cachem_1.0.8 geneplotter_1.72.0
## [16] jsonlite_1.8.0 annotate_1.72.0 cluster_2.1.4
## [19] png_0.1-8 compiler_4.1.0 httr_1.4.7
## [22] Matrix_1.6-3 fastmap_1.1.1 cli_3.6.0
## [25] htmltools_0.5.7 tools_4.1.0 gtable_0.3.0
## [28] glue_1.6.2 GenomeInfoDbData_1.2.7 Rcpp_1.0.8.3
## [31] cellranger_1.1.0 jquerylib_0.1.4 vctrs_0.5.2
## [34] Biostrings_2.62.0 iterators_1.0.14 xfun_0.41
## [37] irlba_2.3.5.1 timechange_0.2.0 lifecycle_1.0.4
## [40] XML_3.99-0.9 dendextend_1.17.1 zlibbioc_1.40.0
## [43] scales_1.2.0 vroom_1.5.7 hms_1.1.3
## [46] parallel_4.1.0 rematch2_2.1.2 RColorBrewer_1.1-3
## [49] prismatic_1.1.0 yaml_2.3.5 memoise_2.0.1
## [52] gridExtra_2.3 sass_0.4.7 stringi_1.7.6
## [55] RSQLite_2.2.10 SQUAREM_2021.1 highr_0.10
## [58] genefilter_1.76.0 foreach_1.5.2 BiocParallel_1.28.3
## [61] truncnorm_1.0-9 shape_1.4.6 rlang_1.1.2
## [64] pkgconfig_2.0.3 bitops_1.0-7 evaluate_0.15
## [67] lattice_0.20-45 invgamma_1.1 labeling_0.4.3
## [70] bit_4.0.5 tidyselect_1.2.0 magrittr_2.0.3
## [73] R6_2.5.1 magick_2.7.4 generics_0.1.2
## [76] DelayedArray_0.20.0 DBI_1.1.2 pillar_1.7.0
## [79] withr_2.5.0 mixsqp_0.3-48 survival_3.5-7
## [82] KEGGREST_1.34.0 RCurl_1.98-1.6 crayon_1.5.1
## [85] utf8_1.2.2 tzdb_0.2.0 viridis_0.6.2
## [88] GetoptLong_1.0.5 locfit_1.5-9.5 blob_1.2.4
## [91] digest_0.6.29 xtable_1.8-4 munsell_0.5.0
## [94] viridisLite_0.4.2 bslib_0.5.1