Last updated: 2025-12-12

Checks: 4 2

Knit directory: ~/Documents/DEPMAPforBRCA1P1revisions/

This reproducible R Markdown analysis was created with workflowr (version 1.7.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

Environment: empty

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.

Seed: set.seed(12345)

The command set.seed(12345) 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.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: absolute

Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.

absolute	relative
C:/Users/maryamshariff/Documents/DEPMAPforBRCA1P1revisions/HCC1806 (ACH-000624) log2(TPM+1) vs MDAMB231 (ACH-000768) log2(TPM+1) (1) (1).png	HCC1806 (ACH-000624) log2(TPM+1) vs MDAMB231 (ACH-000768) log2(TPM+1) (1) (1).png
~/DEPMAPforBRCA1P1revisions	.

Repository version: no version control

Tracking code development and connecting the code version to the results is critical for reproducibility. To start using Git, open the Terminal and type git init in your project directory.

This project is not being versioned with Git. To obtain the full reproducibility benefits of using workflowr, please see ?wflow_start.

For “Regulation of Antiviral and Antitumor Immunity by the BRCA1 Pseudogene in Human Cancers” we have received reviewer comment 13 as stated “The switch from MDA-MB-231 to HCC-1806 cells is justified by the complete growth inhibition observed in the former; however, a more thorough explanation of how HCC-1806 better models the biology of BRCA1P1 loss would be helpful. Differences in baseline BRCA1P1 expression, immune interactions, and antiviral pathway activation between these cell lines should be addressed.”

We planned to address this comment by checking BRCA1P1 expression, immune interactions, and antiviral pathway activation between MDA-MB-231 and HCC-1806 cells.

This document is split up by the following sections: Section 1 : BRCA1P1 expression Section 2: Antiviral gene expression Section 3: Gene Ontology (GO) expression Section 3a: Gene Enrichment Analysis (GSEA) for GO Section 3b: Over-Representation Analysis (ORA) for GO Section 4: KEGG Pathway Expression Section 4a: GSEA for KEGG Section 4b: ORA for KEGG Section 5: MSigDB Hallmark Pathway Expression Section 5a: GSEA for MSigDB Hallmark Section 5b: ORA for MSigDB Hallmark Section 6: Reactome Pathway Expression Section 6a: GSEA for Reactome Section 6b: ORA for Reactome

Section 1 : BRCA1P1 expression

BRCA1P1 is located at chr17:43168170-43168249 based on UCSC Genome Browser. In this analysis, we examined the copy number status of the BRCA1P1 genomic region in HCC1806 (ACH-000624) and MDA-MB-231 (ACH-000768)—using segment-level copy number data from the DepMap WGS dataset. Because BRCA1P1 is a pseudogene, it is not included in DepMap’s gene-level copy number matrices (which require valid Entrez IDs). Therefore, we extracted the segment-level CNV data (OmicsCNSegmentsWGS.csv), which provides the raw genomic intervals and their associated copy number estimates. For each cell line, we identified the unique CNV segment whose genomic coordinates span the BRCA1P1 locus (chr17:43,168,170–43,168,249). Segment boundaries differ across cell lines because CNV segments are called independently.

Section 2: Antiviral gene expression

DEPMAP + Top 100 Gene Lists

First I wanted to calculate the top 100 up-regulated genes in both MDA-MB-231 and HCC1806 cell lines. This value is calculated by the fold change between the two cell lines (MDA-MB-231 – HCC1806).

Gene expression analyses were performed using publicly available DepMap expression data (DepMap Public 25Q3), comprising log-transformed transcript abundance values [log(TPM + 1)] derived from unstranded RNA-sequencing of human protein-coding genes. Expression profiles corresponding to ACH-000624 (HCC1806) and ACH-000768 (MDA-MB-231) breast cancer cell lines were extracted for downstream analysis. Because only a single expression profile was available for each cell line, statistical differential expression testing was not performed. To enable pathway-level comparison between the two cell lines, genes were ranked using a fold-change metric defined as the difference in log-transformed expression values between MDA-MB-231 and HCC1806 (MDA-MB-231 – HCC1806). This ranking metric is consistent with recommended pre-ranking strategies for Gene Set Enrichment Analysis (GSEA) in settings where only log-scale expression values (and no replicate-based statistics) are available. Because these results are based on single-sample profiles per cell line, enrichment findings should be interpreted as hypothesis-generating rather than formal differential expression statistics.

Warning: The above code chunk cached its results, but it won’t be re-run if previous chunks it depends on are updated. If you need to use caching, it is highly recommended to also set knitr::opts_chunk$set(autodep = TRUE) at the top of the file (in a chunk that is not cached). Alternatively, you can customize the option dependson for each individual chunk that is cached. Using either autodep or dependson will remove this warning. See the knitr cache options for more details.

I validated this ranking using DepMap visualization software. Here I mapped a couple of the top upregulated genes in both HCC1806 and MDA-MB-231 and each gene corresponded to fold change values generated by the DepMap scatterplot.

Graph generated using DepMap Visualization software.

I will return to the top-100 up-regulated genes in each cell line later, but I first wanted to briefly examine the antiviral genes evaluated in the paper. Specifically, IFIH1 (MDA5), TNF (TNF-α), IFNB1 (IFN-β), IL1A (IL-1α), and IFIT3 were used as markers of antiviral gene expression. Here, I assess the expression levels of each of these genes in MDA-MB-231 and HCC1806.

Section 3: Gene Ontology (GO) expression

Gene Ontology Biological Process (GO BP) analysis groups genes based on the broader biological functions they participate in. GO BP evaluates whether entire biological processes are up- or down-regulated between two conditions. In this analysis, GO BP terms summarize coordinated shifts in pathway-level activity between MDA-MB-231 and HCC1806, allowing us to identify which cellular programs (e.g., antiviral responses, cytokine signaling, apoptosis) are enriched across the two cell lines.

Section 3a: Gene Enrichment Analysis (GSEA) for GO

In the Functional Class Scoring Methods, I am choosing to use Gene Set Enrichment Analysis for Gene Ontology. I used clusterProfiler (documentation in references) to build GSEA for MDA-MB-231 vs HCC1806. First I built a ranked gene list sorted in descending order based on the logFC calculated in Section 1. GSEA was performed using gseGO() as a rank-based functional class scoring method, operating on the entire transcriptome ranked by log-transformed expression differences between MDA-MB-231 and HCC1806. Enrichment scores and normalized enrichment scores (NES) were computed to assess whether GO BP gene sets were preferentially enriched toward the top or bottom of the ranked list, with positive NES values for pathways enriched in MDA-MB-231 and negative NES values for pathways enriched in HCC1806.

GSEA using gseGO clusterProfiler function for GO

Section 3b: Over-Representation Analysis - Gene Ontology

Over-Representation Analysis differs from Functional Class Scoring tools like GSEA. This sections covers the use of ORA using three different methods: 1. enrichGO function from clusterProfiler 2. DAVID GO analysis from user inputted gene list (top 100 genes section 1) 3. Enrichr GO analysis from user inputted gene list (top 100 genes section 1)

ORA using enrichGO clusterProfilier function for GO

clusterProfiler provides enrichGO, performing hypergeometric testing with Benjamini–Hochberg correction.

Over-Representation Analysis (ORA) - Gene Ontology from DAVID

DAVID represents a classical, annotation-centric implementation of ORA and relies primarily on hypergeometric or Fisher’s exact tests to evaluate whether GO terms or pathways are overrepresented among an input gene list relative to a background that is typically defined as the genome or a platform-specific gene universe. DAVID emphasizes annotation clustering, grouping related GO terms to reduce redundancy and highlight broader biological themes. Multiple-testing correction is performed using the Benjamini–Hochberg procedure. While DAVID is relatively conservative and useful for identifying high-level biological patterns, it relies on older annotation snapshots and is sensitive to gene list size.

Over-Representation Analysis (ORA) - Gene Ontology from Enrichr

Enrichr performs over-representation analysis using Fisher’s exact test. Although multiple ranking metrics are available in Enrichr, in this analysis only Benjamini–Hochberg–adjusted p-values from Gene Ontology Biological Process were used to facilitate direct comparison with ORA results generated using clusterProfiler and DAVID.

Section 4: KEGG Pathway Expression

KEGG pathway expression refers to examining how groups of genes belonging to defined biological pathways (as annotated in the KEGG database) differ between conditions or cell types. Instead of looking at individual genes, KEGG organizes related genes into pathways such as metabolic processes, signaling networks, or immune responses. By mapping onto these pathways, you can identify which biological processes are more active or suppressed in one cell line compared to another.

Section 4a: GSEA for KEGG

Here I use GSEA using gseKEGG clusterProfiler function for KEGG Pathway Analysis. To run KEGG GSEA, gene symbols are first converted to ENTREZ IDs, which KEGG requires. The resulting ranked vector of logFC values (MDA – HCC) is then passed to gseKEGG, which calculates normalized enrichment scores (NES). A positive NES indicates pathways enriched in MDA-MB-231, while a negative NES indicates pathways enriched in HCC1806.

Section 4b: ORA for KEGG

KEGG pathways can be analyzed with ORA in the same way Gene Ontology terms were evaluated. ORA differs from GSEA in that it operates on a discrete gene list (e.g., genes up in each cell line), rather than on the full ranked gene set. This allows ORA to identify KEGG pathways that are statistically overrepresented among genes showing stronger fold-changes.

This sections covers the use of ORA using three different methods: 1. enrichKEGG function from clusterProfiler 2. DAVID KEGG analysis from user inputted gene list 3. Enrichr KEGG analysis from user inputted gene list

ORA using enrichKEGG clusterProfilier function for Pathway Analysis

Over-Representation Analysis (ORA) - KEGG Pathway Analysis using DAVID

Over-Representation Analysis (ORA) - KEGG Pathway Analysis from Enrichr

Section 5: MSigDB Hallmark Pathway Expression

MSigDB Hallmark pathways summarize major biological programs into 50 well-defined gene sets that reduce noise and redundancy found in larger databases. By applying GSEA to these Hallmark signatures, we can detect coordinated pathway-level differences between MDA-MB-231 and HCC1806, highlighting broad processes such as EMT, immune signaling, and proliferation that distinguish the two cell lines.

Section 5a: GSEA for MSigDB Hallmark

Section 5b: ORA for MSigDB Hallmark

Similar to ORA done before, I used multiple ORA methods for Hallmark pathway analysis, at first I used clusterProfiler::enricher + msigdbr, then after I used ORA downloaded from Enrichr.

ORA with MSigDB HALLMARK (clusterProfiler::enricher + msigdbr)

ORA MsigDB Hallmark ORA from Enrichr

Section 6: Reactome Pathway Expression

Reactome is a manually curated pathway database that organizes human biology into detailed molecular processes, including signaling cascades, immune pathways, cell-cycle events, and metabolism. Because Reactome pathways are structured hierarchically and capture mechanistic relationships between reactions, they provide more granular functional context than KEGG or Hallmark collections.

Section 6a: GSEA for Reactome

Reactome GSEA with ReactomePA + clusterProfiler functions.

Section 6b: ORA for Reactome

Using ReactomePA::enrichPathway + clusterProfiler/enrichplot viz

Enrichr ORA Reactome Pathway

References

Metrics for Ranking Genes From https://www.gsea-msigdb.org/gsea/doc/GSEAUserGuideFrame.html :

When you run the gene set enrichment analysis from the Run GSEA Page, GSEA ranks the genes in the expression dataset and then analyzes that ranked list of genes. You use the Metric for ranking genes parameter to select the metric used to score and rank the genes; the Gene list sorting mode parameter to determine whether to sort the genes using the real (default) or absolute value of the metric score; and the Gene list ordering mode parameter to determine whether to sort the genes in descending (default) or ascending order. This section describes each of the ranking metrics in the drop-down list of the Metric for ranking genes parameter. If your favorite metric is not listed here, you can rank the genes in your dataset using that metric and then use the GSEAPreranked Page to analyze your ranked list of genes. If your dataset contains only one sample, GSEA cannot rank the genes; however, you can rank the genes and then use the GSEAPreranked Page to analyze your ranked list of genes. Three settings in the Algorithms tab of the Preferences window (selected from the File menu in the application menu bar) affect the calculations shown here: ● Use median instead of mean for metrics. For categorical phenotypes, by default, GSEA calculates differential expression based on the mean expression value for each phenotype. To use the median expression value for each phenotype, set the Median for class metrics parameter on the Run GSEA page to True. To always use the median expression value for each phenotype, click on the checkbox to enable the Use median instead of mean for class metrics option in the Preferences window. The new default value takes effect when you restart GSEA. By default, this option is not selected. Note: Using median rather than mean may cause ties in the ranking. ● Fix metrics for low variance. When calculating ranking metrics, the denominator may be zero (0). A denominator of zero (0) causes an error in the analysis unless this option is selected. By default, this option is selected. ● Use biased variances. When calculating ranking metrics, GSEA uses an unbiased variance to calculate standard deviation. Select this option to have GSEA use a biased variance instead. By default, this option is not selected. For categorical phenotypes, GSEA determines a gene’s mean expression value for each phenotype and then uses one of the following metrics to calculate the gene’s differential expression with respect to the two phenotypes. To use median rather than mean expression values, set the Median for class metrics parameter to True, as described above. Diff_of_Classes uses the difference of class means to calculate fold change for log scale data:

where μ is the mean. The larger the fold change, the more distinct the gene expression is in each phenotype and the more the gene acts as a “class marker.”

Log-transformed TPM data From https://depmap.org/portal/data_page/?tab=allData

“Log-transformed TPM (Transcripts Per Million) values at gene level derived from unstranded RNA-seq for protein-coding genes in human”

clusterProfiler from https://yulab-smu.top/biomedical-knowledge-mining-book/index.html

Used enrichGO() for ORA GO and gseGO() for GSEA GO

msigdbr from https://igordot.github.io/msigdbr/articles/msigdbr-intro.html

Tutorial from https://biostatsquid.com/fgsea-tutorial-gsea/ . https://biostatsquid.com/pathway-enrichment-analysis-tutorial-clusterprofiler/

Another Tutorial from https://alexslemonade.github.io/refinebio-examples/03-rnaseq/pathway-analysis_rnaseq_01_ora.html

ReactomePA Tutorial https://bioconductor.statistik.uni-dortmund.de/packages/3.8/bioc/vignettes/ReactomePA/inst/doc/ReactomePA.html

DAVID from https://davidbioinformatics.nih.gov/

Inputted 100 top genes generated from fold change to calculate ORA using DAVID algorithim.

Enrichr from https://maayanlab.cloud/Enrichr/

Inputted 100 top genes generated from fold change to calculate ORA using Enrichr platform.

sessionInfo()

R version 4.5.1 (2025-06-13 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] lubridate_1.9.4        purrr_1.1.0            tidyr_1.3.1           
 [4] tidyverse_2.0.0        tibble_3.3.0           stringr_1.5.2         
 [7] RColorBrewer_1.1-3     readxl_1.4.5           readr_2.1.5           
[10] ReactomePA_1.52.0      pheatmap_1.0.13        org.Hs.eg.db_3.21.0   
[13] AnnotationDbi_1.70.0   IRanges_2.42.0         S4Vectors_0.46.0      
[16] Biobase_2.68.0         BiocGenerics_0.54.0    generics_0.1.4        
[19] openxlsx_4.2.8         msigdbr_25.1.1         ggupset_0.4.1         
[22] ggplot2_4.0.0          forcats_1.0.1          fgsea_1.34.2          
[25] enrichplot_1.28.4      dplyr_1.1.4            DOSE_4.2.0            
[28] clusterProfiler_4.16.0

loaded via a namespace (and not attached):
  [1] rstudioapi_0.17.1       jsonlite_2.0.0          magrittr_2.0.3         
  [4] ggtangle_0.0.9          farver_2.1.2            rmarkdown_2.30         
  [7] fs_1.6.6                vctrs_0.6.5             memoise_2.0.1          
 [10] ggtree_3.16.3           htmltools_0.5.8.1       curl_7.0.0             
 [13] cellranger_1.1.0        gridGraphics_0.5-1      sass_0.4.10            
 [16] bslib_0.9.0             plyr_1.8.9              cachem_1.1.0           
 [19] igraph_2.2.1            lifecycle_1.0.4         pkgconfig_2.0.3        
 [22] Matrix_1.7-3            R6_2.6.1                fastmap_1.2.0          
 [25] gson_0.1.0              GenomeInfoDbData_1.2.14 digest_0.6.37          
 [28] aplot_0.2.9             patchwork_1.3.2         rprojroot_2.1.1        
 [31] RSQLite_2.4.3           labeling_0.4.3          timechange_0.3.0       
 [34] httr_1.4.7              polyclip_1.10-7         compiler_4.5.1         
 [37] bit64_4.6.0-1           withr_3.0.2             graphite_1.54.0        
 [40] S7_0.2.0                BiocParallel_1.42.1     viridis_0.6.5          
 [43] DBI_1.2.3               ggforce_0.5.0           R.utils_2.13.0         
 [46] MASS_7.3-65             rappdirs_0.3.3          tools_4.5.1            
 [49] otel_0.2.0              ape_5.8-1               zip_2.3.3              
 [52] httpuv_1.6.16           R.oo_1.27.1             glue_1.8.0             
 [55] nlme_3.1-168            GOSemSim_2.34.0         promises_1.5.0         
 [58] grid_4.5.1              reshape2_1.4.4          snow_0.4-4             
 [61] gtable_0.3.6            tzdb_0.5.0              R.methodsS3_1.8.2      
 [64] hms_1.1.4               data.table_1.17.8       tidygraph_1.3.1        
 [67] XVector_0.48.0          ggrepel_0.9.6           pillar_1.11.1          
 [70] vroom_1.6.7             yulab.utils_0.2.2       babelgene_22.9         
 [73] later_1.4.4             splines_4.5.1           tweenr_2.0.3           
 [76] treeio_1.32.0           lattice_0.22-7          bit_4.6.0              
 [79] tidyselect_1.2.1        GO.db_3.21.0            Biostrings_2.76.0      
 [82] reactome.db_1.92.0      knitr_1.50              git2r_0.36.2           
 [85] gridExtra_2.3           xfun_0.53               graphlayouts_1.2.2     
 [88] stringi_1.8.7           UCSC.utils_1.4.0        workflowr_1.7.2        
 [91] lazyeval_0.2.2          ggfun_0.2.0             yaml_2.3.10            
 [94] evaluate_1.0.5          codetools_0.2-20        ggraph_2.2.2           
 [97] qvalue_2.40.0           graph_1.86.0            ggplotify_0.1.3        
[100] cli_3.6.5               jquerylib_0.1.4         Rcpp_1.1.0             
[103] GenomeInfoDb_1.44.3     png_0.1-8               parallel_4.5.1         
[106] assertthat_0.2.1        blob_1.2.4              viridisLite_0.4.2      
[109] tidytree_0.4.6          scales_1.4.0            crayon_1.5.3           
[112] rlang_1.1.6             cowplot_1.2.0           fastmatch_1.1-6        
[115] KEGGREST_1.48.1

Comment 13 Using DepMap Data BRCA1P1 Revisions

Maryam Shariff

2025-12-02