ATAC-seq and RNA-seq Analysis of Triple-Negative Breast Cancer

Introduction

The goal of this project is to use raw dataset to perform upstream analysis and bulk ATAC analysis, then analyze the peaks to identify genes that may play important roles in triple-negative breast cancer (TNBC). RNA sequencing was further used to validate the findings from ATAC sequencing. To start this project, the raw data was obtained from NCBI project PRJNA1069018, and a custom Nextflow pipeline was used to perform the upstream bulk ATAC analysis, including preprocessing of raw reads, alignment, peak calling, and generation of peak files for downstream analysis. After that, differential accessibility analysis and peak annotation were performed to identify regions with altered chromatin accessibility and connect these regions to nearby genes. Functional enrichment analysis and gene set enrichment analysis were then used to examine the biological functions related to these genes. Lastly, RNA-seq differential expression analysis was performed and compared with the ATAC-seq results as a validation step.

Loading libraries

library(DiffBind)
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)
library(GenomeInfoDb)
library(GenomicRanges)
library(AnnotationDbi)
library(clusterProfiler)
library(DESeq2)
library(enrichplot)
library(ggplot2)
library(fgsea)
library(BiocParallel)

Loading files

atac_sample_id <- setdiff(seq(27726627, 27726645), 27726635)
peak_path <- "C:/Users/fchen/OneDrive/Desktop/TNBC_data/peak_files/SRR"
bam_path <- "C:/Users/fchen/OneDrive/Desktop/TNBC_data/bam_files/SRR"
peak_suffix <- "_peaks.narrowPeak"
bam_suffix <- "_final_sorted.bam"

peak_files <- paste0(peak_path, atac_sample_id, peak_suffix)

bam_files <- paste0(bam_path, atac_sample_id, bam_suffix)

Loading meta file

meta_file <- read.csv("C:/Users/fchen/OneDrive/Desktop/TNBC_data/TNBC_meta.csv")

dba_sheet <- data.frame(
  SampleID = paste(meta_file$Patient, meta_file$Treatment, sep = "_"),
  Condition = meta_file$Treatment, ## Condition are normal or tumor
  Tissue = meta_file$Patient, ## Each patient will have two sample, one normal and one tumor
  Factor = "ATAC",
  Replicate = meta_file$Replicate,
  bamReads = bam_files,
  Peaks = peak_files,
  PeakCaller = "macs"
)

Create dba object and peak counts

atac_dba <- dba(sampleSheet = dba_sheet)
dba.show(atac_dba)

atac_count <- dba.count(
                        atac_dba,
                        summits = 250,
                        minOverlap = 2
                        )
dba.show(atac_count)

Heatmap and PCA plot

dba.plotHeatmap(atac_count)

dba.plotPCA(atac_count, DBA_CONDITION, label = DBA_ID) ## PCA plot based on the condition

The heatmap shows the correlation between samples based on their chromatin accessibility signals. Most samples show relatively high correlation with each other, which indicates that the overall chromatin accessibility patterns are generally similar across samples. In addition, the clustering result shows that most tumor and normal samples tend to group into two clusters based on their condition, suggesting that chromatin accessibility profiles are different between tumor and normal tissues.

The PCA plot also shows two conditions. Most of the normal samples form one cluster, while tumor samples tend to spread out more across the plot. This suggests that tumor samples have more diverse chromatin accessibility patterns compared with normal samples. This matches the characteristic of triple-negative breast cancer (TNBC), which is known to have higher heterogeneity.

Create dba contrast for future comparison analysis

atac_contrast <- dba.contrast(
                              atac_count,
                              categories = DBA_CONDITION,
                              block = DBA_TISSUE
                              )
dba.show(atac_contrast, bContrasts = TRUE) ## Total of 18 samples, 9 tumor and 9 normal

Perform comparison analysis between tumor and normal cells

atac_analyze <- dba.analyze(atac_contrast, method = DBA_DESEQ2)
dba.plotMA(atac_analyze)

dba.plotVolcano(atac_analyze)

Both the MA plot and volcano plot indicate that tumor samples contain more regions with increased chromatin accessibility compared with regions showing decreased accessibility. In total, 1435 regions show significant differential accessibility (FDR < 0.05) between tumor and normal samples. Only 5 out of 1435 show decreased accessibility in tumor samples. This is consistent with many previous studies suggesting that cancer cells often exhibit more open chromatin regions.

Split genes with increased accessibility and genes with decrease accessibility

atac_df <- as.data.frame(dba.report(atac_analyze, th = 1))
closed_region <- atac_df[atac_df$Fold < 0 & atac_df$FDR <= 0.05, ] ## significant closed region
open_region <- atac_df[atac_df$Fold > 0 & atac_df$FDR <= 0.05, ] ## significant open region


all_gr <- makeGRangesFromDataFrame(atac_df, keep.extra.columns = TRUE)
closed_gr <- makeGRangesFromDataFrame(closed_region, keep.extra.columns = TRUE)
open_gr <- makeGRangesFromDataFrame(open_region, keep.extra.columns = TRUE)

Prepare datasets for functional enrichment analysis

## set all granges to same style as the reference data
seqlevelsStyle(all_gr) <- seqlevelsStyle(TxDb.Hsapiens.UCSC.hg38.knownGene)
seqlevelsStyle(closed_gr) <- seqlevelsStyle(TxDb.Hsapiens.UCSC.hg38.knownGene)
seqlevelsStyle(open_gr) <- seqlevelsStyle(TxDb.Hsapiens.UCSC.hg38.knownGene)

## run peak annotation for all genes
anno_all <- annotatePeak(all_gr,
             TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene, 
             annoDb = "org.Hs.eg.db",
             tssRegion = c(-3000, 3000), 
             verbose = FALSE)

## >> Using Genome: hg38 ...

## >> Using Genome: hg38 ...

## >> Using Genome: hg38 ...

anno_closed <- annotatePeak(closed_gr,
             TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene, 
             annoDb = "org.Hs.eg.db",
             tssRegion = c(-3000, 3000), 
             verbose = FALSE)

## >> Using Genome: hg38 ...

## >> Using Genome: hg38 ...

## >> Using Genome: hg38 ...

anno_open <- annotatePeak(open_gr,
                          TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene,
                          annoDb = "org.Hs.eg.db",
                          tssRegion = c(-3000, 3000),
                          verbose = FALSE)

## >> Using Genome: hg38 ...

## >> Using Genome: hg38 ...

## >> Using Genome: hg38 ...

## convert to dataframe for later enrichment analysis 
df_all <- as.data.frame(anno_all)
df_closed <- as.data.frame(anno_closed)
df_open <- as.data.frame(anno_open)
plotAnnoPie(anno_closed)

plotAnnoPie(anno_open)

plotDistToTSS(anno_closed)

plotDistToTSS(anno_open)

The annotation pie charts show the genomic composition of the peaks. For peaks with increased accessibility, around 20–25% of the peaks are located in promoter regions, while a larger portion of peaks are located in intronic regions (around 40%), with the remaining peaks distributed across exons and distal intergenic regions. For regions with decreased accessibility, the annotation pie chart is less detailed. This is expected because only a small number of peaks show decreased accessibility, and the genomic distribution of closed regions appears generally similar to that of open regions.

The TSS distribution plot for peaks with increased accessibility shows that many peaks are located several kilobases away from the TSS, while around one-third are located close to promoter regions. For regions with decreased accessibility, because only a small number of peaks are observed, the plot is simplified and shows a more general distribution relative to the TSS.

Overrepresentation (ORA) functional enrichment analysis

## filter out the gene symbols locate within the promoter region
symbol_closed <- unique(df_closed$SYMBOL[!is.na(df_closed$SYMBOL) & abs(df_closed$distanceToTSS) <= 3000])
symbol_open <- unique(df_open$SYMBOL[!is.na(df_open$SYMBOL) & abs(df_open$distanceToTSS) <= 3000])
symbol_background <- unique(df_all$SYMBOL[!is.na(df_all$SYMBOL) & abs(df_all$distanceToTSS) <= 3000])

## load hallmark gmt
hall_gmt <- read.gmt("C:\\Users\\fchen\\Downloads\\h.all.v2026.1.Hs.symbols.gmt")

## create two hall mark enrichment 
hallmark_closed <- enricher(symbol_closed, 
                            TERM2GENE = hall_gmt, 
                            universe = symbol_background)
hallmark_open <- enricher(gene = symbol_open, 
                          TERM2GENE = hall_gmt, 
                          universe = symbol_background,
                          pvalueCutoff = 1,
                          qvalueCutoff = 1)

df_hallmark_open <- as.data.frame(hallmark_open)

## dot plot for hallmark function enrichment for open regions
dotplot(hallmark_open, showCategory = 10) +
  ggtitle("Hallmark enrichment (Open promoters)")

## count how many rows are significant for hallmark analysis
nrow(df_hallmark_open[df_hallmark_open$p.adjust < 0.05,])

## [1] 0

##  filter out the ensembl id locate within the promoter region
open_ensembl <- unique(na.omit(df_open$ENSEMBL[!is.na(df_open$geneId) & abs(df_open$distanceToTSS) <= 3000]))
background_ensembl <- unique(na.omit(df_all$ENSEMBL[!is.na(df_all$geneId) & abs(df_all$distanceToTSS) <= 3000]))

## convert the ensembl id into entrez id for KEGG and GO analysis
open_entrezid <- mapIds(org.Hs.eg.db,
                 keys = open_ensembl,
                 column = "ENTREZID",
                 keytype = "ENSEMBL")

background_entrezid <-  mapIds(org.Hs.eg.db,
                               keys = background_ensembl,
                               column = "ENTREZID",
                               keytype = "ENSEMBL")


## KEGG functional enrichment analysis
kegg_open <- enrichKEGG(gene = open_entrezid,
                    universe = background_entrezid,
                    organism = "hsa",
                    pAdjustMethod = "BH",
                    pvalueCutoff = 1,
                    qvalueCutoff = 1)

dotplot(kegg_open, showCategory = 10) +
  ggtitle("KEGG enrichment (Open promoters)")

df_kegg_open <- as.data.frame(kegg_open)
nrow(df_kegg_open[df_kegg_open$p.adjust < 0.05,])

## [1] 1

## GO functional enrichment analysis
go_open <- enrichGO(gene = open_entrezid,
                    universe = background_entrezid,
                    OrgDb = org.Hs.eg.db,
                    ont = "BP",              # BP/CC/MF
                    pAdjustMethod = "BH",
                    pvalueCutoff = 1,
                    qvalueCutoff = 1)

dotplot(go_open, showCategory = 10) +
  ggtitle("GO enrichment (Open promoters)")

df_go_open <- as.data.frame(go_open@result)
nrow(df_go_open[df_go_open$p.adjust < 0.05,])

## [1] 0

Since there are only a very small number of genes related to regions showing decreased accessibility, performing functional enrichment analysis on those genes will not perform well. Therefore, I skipped the enrichment analysis for the closed regions and focused only on regions with increased accessibility. For the open regions, I performed three different enrichment analyses: Hallmark pathway enrichment, Gene Ontology (GO) enrichment, and KEGG pathway enrichment. Only genes located in promoter regions (within 3 kb of the TSS) were used for the analysis.

From the results, only one pathway appears to be statistically significant, while most pathways do not pass the adjusted p-value threshold (0.05). This suggests that overrepresentation analysis (ORA) may not be the best method for this dataset.

Although most of the results are not statistically significant, we can still see some cancer-related pathways appearing in the results. For example, the estrogen response pathway from Hallmark, the tight junction pathway from KEGG, and some stress-related pathways from GO appear among the top results. Therefore, to further verify the results, I also performed GSEA analysis using ranked genes.

Gene-Set Enrichment Analysis (GSEA)

## replace all 0 p.value with 1e-100 (extreme small number) so that -log10(p.value) can still return numeric value
df_all$p.value <- ifelse(
  df_all$p.value == 0, 1e-100, df_all$p.value
)

## define ranking score for each peak by combining fold direction and significance
df_all$ranking <- sign(df_all$Fold) * (-log10(df_all$p.value))

## create data frame grouped by gene symbol, and keep the highest ranking score for each gene, This avoids one gene being linked to multiple peaks
df_gene_rank <- data.frame(
  SYMBOL = names(tapply(df_all$ranking, df_all$SYMBOL, max, na.rm = TRUE)),
  rank_score = tapply(df_all$ranking, df_all$SYMBOL, max, na.rm = TRUE)
)

## sort genes by ranking score from low to high
df_gene_rank <- df_gene_rank[order(df_gene_rank$rank_score), ]

## plot overall distribution of ranking score
plot(df_gene_rank$rank_score)

## dot plot for first 50 genes with lowest ranking score
ggplot(df_gene_rank[1:50, ], aes(x = SYMBOL, y = rank_score)) +
  geom_point() +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

## split hallmark gmt file into list format, each pathway contains its own genes
hallmark_gene_sets <- split(hall_gmt$gene, hall_gmt$term)

## create named numeric vector for GSEA input and sort from high to low
ranks <- df_gene_rank$rank_score
names(ranks) <- df_gene_rank$SYMBOL
ranks <- sort(ranks, decreasing = TRUE)

## run hallmark GSEA using ranked gene list
gsea_hallmark_df <- fgsea(
  pathways = hallmark_gene_sets,
  stats = ranks,
  scoreType = "std",
  minSize = 10,
  nproc = 1
)

##   |                                                                              |                                                                      |   0%  |                                                                              |======                                                                |   9%  |                                                                              |=============                                                         |  18%  |                                                                              |===================                                                   |  27%  |                                                                              |=========================                                             |  36%  |                                                                              |================================                                      |  45%  |                                                                              |======================================                                |  55%  |                                                                              |=============================================                         |  64%  |                                                                              |===================================================                   |  73%  |                                                                              |=========================================================             |  82%  |                                                                              |================================================================      |  91%  |                                                                              |======================================================================| 100%

## count how many hallmark pathways are significant
sum(gsea_hallmark_df[, padj < 0.05])

## [1] 7

## top 5 positively enriched pathways
top5 <- gsea_hallmark_df[ES > 0][head(order(padj), n = 5), ]
pathway_open_hallmark <- top5[, pathway]

## check how many pathways are selected for open and closed regions
length(pathway_open_hallmark)

## [1] 5

## show names of top positive hallmark pathways
pathway_open_hallmark

## [1] "HALLMARK_ESTROGEN_RESPONSE_LATE"  "HALLMARK_ESTROGEN_RESPONSE_EARLY"
## [3] "HALLMARK_ANDROGEN_RESPONSE"       "HALLMARK_IL2_STAT5_SIGNALING"    
## [5] "HALLMARK_KRAS_SIGNALING_UP"

## create GSEA summary table for top open hallmark pathways
plotGseaTable(
  hallmark_gene_sets[pathway_open_hallmark],
  stats = ranks,
  fgseaRes = gsea_hallmark_df,
  gseaParam = 0.5
)

## plot enrichment curve for the most significant open hallmark pathway
plotEnrichment(hallmark_gene_sets[[pathway_open_hallmark[1]]], ranks) + 
  labs(title = pathway_open_hallmark[1])

## plot cell growth related pathway
plotEnrichment(hallmark_gene_sets[["HALLMARK_KRAS_SIGNALING_UP"]], ranks) +
  labs(title = "HALLMARK_KRAS_SIGNALING_UP")

## plot inflammatory related pathway
plotEnrichment(hallmark_gene_sets[["HALLMARK_INFLAMMATORY_RESPONSE"]], ranks) +
  labs(title = "HALLMARK_INFLAMMATORY_RESPONSE")

## visualize hallmark GSEA result
ggplot(top5, aes(x = NES, y = pathway)) +
  geom_point(aes(size = size, color = padj)) +
  scale_color_gradient(low = "red", high = "blue") +
  labs (title = "Hallmark GSEA")

## convert SYMBOL into ENTREZ ID for GO and KEGG
symbol_gsea <- mapIds(
  org.Hs.eg.db,
  keys = df_gene_rank$SYMBOL,
  column = "ENTREZID",
  keytype = "SYMBOL",
  multiVals = "first"
)

## add ENTREZ ID into ranked gene table
df_gene_rank_entrezid <- df_gene_rank
df_gene_rank_entrezid$ENTREZID <- unname(symbol_gsea[df_gene_rank_entrezid$SYMBOL])

## remove genes without ENTREZ ID
df_gene_rank_entrezid <- df_gene_rank_entrezid[
  !is.na(df_gene_rank_entrezid$ENTREZID),
]

## create ranked numeric vector using ENTREZ IDs
## if multiple symbols map to same ENTREZ ID, keep the highest rank score
ranks_entrez <- tapply(
  df_gene_rank_entrezid$rank_score,
  df_gene_rank_entrezid$ENTREZID,
  max
)

ranks_entrez <- setNames(as.numeric(ranks_entrez), names(ranks_entrez))
ranks_entrez <- sort(ranks_entrez, decreasing = TRUE)

## register serial parameter to avoid parallel issue
register(SerialParam())

## run GO GSEA using biological process ontology
gsea_go <- gseGO(
  geneList = ranks_entrez,
  OrgDb = org.Hs.eg.db,
  keyType = "ENTREZID",
  ont = "BP",
  minGSSize = 10,
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH"
)

## check the minimum enrichment score
## if minimum ES is still positive, then all enriched GO pathways are on open chromatin side
min(gsea_go$enrichmentScore)

## [1] 0.3244336

## convert GO GSEA result into data frame
gsea_go_df <- as.data.frame(gsea_go)

## keep only positive GO pathways
open_go <- gsea_go_df[gsea_go_df$enrichmentScore > 0, ]

## select top 5 open GO pathways by adjusted p value
pathway_open_go <- head(open_go[order(open_go$p.adjust), "Description"], n = 5)

## visualize GO GSEA result
dotplot(gsea_go, showCategory = pathway_open_go, x = "NES") +
  ggtitle("GO GSEA")

## run KEGG GSEA using ranked ENTREZ gene list
gsea_kegg <- gseKEGG(
  geneList = ranks_entrez,
  organism = "hsa",
  pAdjustMethod = "BH",
  minGSSize = 10,
  pvalueCutoff = 0.05
)

## convert KEGG GSEA result into data frame
gsea_kegg_df <- as.data.frame(gsea_kegg)

## keep only positive KEGG pathways
open_kegg <- gsea_kegg_df[gsea_kegg_df$enrichmentScore > 0, ]

## select top 5 open KEGG pathways by adjusted p value
pathway_open_kegg <- head(open_kegg[order(open_kegg$p.adjust), "Description"], n = 5)

## visualize KEGG GSEA result
dotplot(gsea_kegg, showCategory = pathway_open_kegg, x = "NES") +
  ggtitle("KEGG GSEA")

To perform GSEA analysis, I first defined a ranking score for each gene based on both the fold change and the p-value. The equation used is sign(Fold change) × −log10(p-value). This allows genes associated with increased chromatin accessibility to receive higher ranking scores.

Similar to what was done in the ORA analysis, genes associated with decreased accessibility were not included in this analysis. This is because the number of genes related to regions with decreased accessibility is very small, and performing enrichment analysis on such a small gene set would not produce reliable results. Therefore, the analysis focused only on genes associated with regions showing increased accessibility.

Compared to the ORA analysis, GSEA identified more statistically significant pathways. The results suggest that regions with increased chromatin accessibility in tumor samples are associated with pathways related to cell signaling, hormone response, immune response, and cell adhesion or migration, which are consistent with known characteristics of triple-negative breast cancer from previous studies.

RNA sequencing

## load RNA count matrix
df_featurecount <- read.csv("RNA_Counts.csv")

## replace "." in sample names with "-" so that sample names match metadata
colnames(df_featurecount) <- gsub("\\.", "-", colnames(df_featurecount))

## move gene id column into row names
rownames(df_featurecount) <- df_featurecount$X
df_featurecount$X <- NULL

## load RNA metadata
metadata <- read.csv("meta_data.csv")

## set sample id as row names so that metadata matches count matrix columns
rownames(metadata) <- metadata$sample_id

## keep only condition column for DESeq2 analysis
metadata <- subset(metadata, select = "condition")

## set Normal as reference group and Tumor as comparison group
metadata$condition <- factor(metadata$condition, levels = c("Normal", "Tumor"))

## create DESeq2 object using raw count matrix and sample condition
des <- DESeqDataSetFromMatrix(
  countData = df_featurecount,
  colData = metadata,
  design = ~condition
)

## remove very low count genes to reduce noise
des <- des[rowSums(counts(des)) > 10, ]

## run DESeq2 differential expression analysis
des <- DESeq(des)

## get differential expression result for Tumor versus Normal
result <- results(
  des,
  contrast = c("condition", "Tumor", "Normal"),
  alpha = 1e-5
)

## remove ensembl version number so that RNA gene ids match ATAC gene ids
rownames(result) <- sub("\\..*$", "", rownames(result))

## define significantly upregulated and downregulated RNA genes
rna_up <- rownames(result)[which(result$padj < 0.05 & result$log2FoldChange > 0)]
rna_down <- rownames(result)[which(result$padj < 0.05 & result$log2FoldChange < 0)]

## apply variance stabilizing transformation for PCA plot
variance <- vst(des, blind = FALSE)

## check overlap between ATAC open genes and RNA downregulated genes
intersect(open_ensembl, rna_down)

##  [1] "ENSG00000203721" "ENSG00000127914" "ENSG00000167291" "ENSG00000070731"
##  [5] "ENSG00000133056" "ENSG00000260710" "ENSG00000139631" "ENSG00000134313"
##  [9] "ENSG00000203985" "ENSG00000140443" "ENSG00000166780" "ENSG00000143507"
## [13] "ENSG00000009830" "ENSG00000235750" "ENSG00000065809" "ENSG00000132554"
## [17] "ENSG00000027075" "ENSG00000165084" "ENSG00000167861" "ENSG00000130779"
## [21] "ENSG00000167202" "ENSG00000150764" "ENSG00000050405" "ENSG00000205726"
## [25] "ENSG00000167549" "ENSG00000116977" "ENSG00000058091" "ENSG00000115648"
## [29] "ENSG00000245750" "ENSG00000170190" "ENSG00000140577" "ENSG00000132746"
## [33] "ENSG00000263586" "ENSG00000182704" "ENSG00000124749" "ENSG00000233559"
## [37] "ENSG00000137486" "ENSG00000197496" "ENSG00000132561" "ENSG00000108846"
## [41] "ENSG00000198369" "ENSG00000100968" "ENSG00000126368" "ENSG00000166548"
## [45] "ENSG00000047346" "ENSG00000141905" "ENSG00000204913" "ENSG00000157985"
## [49] "ENSG00000169018" "ENSG00000173334" "ENSG00000108582" "ENSG00000144476"
## [53] "ENSG00000196535" "ENSG00000141298" "ENSG00000106100" "ENSG00000185630"
## [57] "ENSG00000137266" "ENSG00000143412" "ENSG00000137815" "ENSG00000139629"
## [61] "ENSG00000228013" "ENSG00000196639" "ENSG00000139865" "ENSG00000233806"
## [65] "ENSG00000109339" "ENSG00000123472" "ENSG00000177182" "ENSG00000285043"
## [69] "ENSG00000168792" "ENSG00000144228" "ENSG00000160551" "ENSG00000133789"
## [73] "ENSG00000203799" "ENSG00000114248" "ENSG00000065989" "ENSG00000006282"
## [77] "ENSG00000183077" "ENSG00000234690" "ENSG00000048707" "ENSG00000010704"
## [81] "ENSG00000203666" "ENSG00000126653" "ENSG00000107263"

## check overlap between ATAC open genes and RNA upregulated genes
intersect(open_ensembl, rna_up)

##   [1] "ENSG00000141837" "ENSG00000249859" "ENSG00000140548" "ENSG00000159259"
##   [5] "ENSG00000227496" "ENSG00000034677" "ENSG00000104442" "ENSG00000196684"
##   [9] "ENSG00000117597" "ENSG00000159915" "ENSG00000101391" "ENSG00000064313"
##  [13] "ENSG00000171302" "ENSG00000285677" "ENSG00000169896" "ENSG00000147687"
##  [17] "ENSG00000198125" "ENSG00000163545" "ENSG00000230526" "ENSG00000128039"
##  [21] "ENSG00000166278" "ENSG00000157637" "ENSG00000176920" "ENSG00000141564"
##  [25] "ENSG00000006453" "ENSG00000167881" "ENSG00000152457" "ENSG00000123992"
##  [29] "ENSG00000116260" "ENSG00000173227" "ENSG00000149346" "ENSG00000127948"
##  [33] "ENSG00000022567" "ENSG00000196187" "ENSG00000253802" "ENSG00000172965"
##  [37] "ENSG00000167543" "ENSG00000143156" "ENSG00000225656" "ENSG00000206053"
##  [41] "ENSG00000196890" "ENSG00000135127" "ENSG00000008283" "ENSG00000173715"
##  [45] "ENSG00000181085" "ENSG00000227036" "ENSG00000185467" "ENSG00000167522"
##  [49] "ENSG00000204387" "ENSG00000223882" "ENSG00000253302" "ENSG00000179583"
##  [53] "ENSG00000109062" "ENSG00000156639" "ENSG00000182446" "ENSG00000154380"
##  [57] "ENSG00000108669" "ENSG00000115970" "ENSG00000143512" "ENSG00000147509"
##  [61] "ENSG00000130402" "ENSG00000198520" "ENSG00000221955" "ENSG00000175115"
##  [65] "ENSG00000177398" "ENSG00000278677" "ENSG00000197153" "ENSG00000134824"
##  [69] "ENSG00000197472" "ENSG00000167011" "ENSG00000100401" "ENSG00000261716"
##  [73] "ENSG00000265817" "ENSG00000215148" "ENSG00000117450" "ENSG00000167703"
##  [77] "ENSG00000254337" "ENSG00000087157" "ENSG00000117009" "ENSG00000132698"
##  [81] "ENSG00000125864" "ENSG00000147697" "ENSG00000152969" "ENSG00000156515"
##  [85] "ENSG00000284690" "ENSG00000143375" "ENSG00000005001" "ENSG00000164932"
##  [89] "ENSG00000103522" "ENSG00000141741" "ENSG00000198182" "ENSG00000263961"
##  [93] "ENSG00000171813" "ENSG00000286522" "ENSG00000160439" "ENSG00000238160"
##  [97] "ENSG00000158769" "ENSG00000162078" "ENSG00000173638" "ENSG00000273983"
## [101] "ENSG00000066827" "ENSG00000175866" "ENSG00000133816" "ENSG00000261150"
## [105] "ENSG00000170873" "ENSG00000124215" "ENSG00000173156" "ENSG00000105289"
## [109] "ENSG00000141738" "ENSG00000249816" "ENSG00000142002" "ENSG00000168014"
## [113] "ENSG00000233791" "ENSG00000161929" "ENSG00000143367" "ENSG00000253644"
## [117] "ENSG00000104824" "ENSG00000122643" "ENSG00000163462" "ENSG00000163159"
## [121] "ENSG00000180999" "ENSG00000213937" "ENSG00000207606" "ENSG00000161405"
## [125] "ENSG00000105676" "ENSG00000104529" "ENSG00000184357" "ENSG00000104450"
## [129] "ENSG00000117222" "ENSG00000251002" "ENSG00000101109" "ENSG00000140506"
## [133] "ENSG00000048052" "ENSG00000167889" "ENSG00000005486" "ENSG00000164880"
## [137] "ENSG00000285758" "ENSG00000282849" "ENSG00000173120" "ENSG00000230479"
## [141] "ENSG00000116580" "ENSG00000067057" "ENSG00000225783" "ENSG00000163162"
## [145] "ENSG00000095380" "ENSG00000026508" "ENSG00000246526" "ENSG00000255004"
## [149] "ENSG00000125450" "ENSG00000167880" "ENSG00000204380" "ENSG00000128340"
## [153] "ENSG00000130182" "ENSG00000118194" "ENSG00000136450" "ENSG00000187605"
## [157] "ENSG00000167642" "ENSG00000236790" "ENSG00000141526" "ENSG00000167281"

## plot PCA and MA plot for RNA differential expression result
dev.new()
plotPCA(variance)
plotMA(result)

## calculate the probability of overlap between ATAC open genes and RNA upregulated genes happening by chance
phyper(length(intersect(open_ensembl, rna_up)) - 1, length(unique(open_ensembl)), nrow(result) - length(unique(open_ensembl)), length(unique(rna_up)), lower.tail = FALSE)

## [1] 3.738562e-05

## calculate the probability of overlap between ATAC open genes and RNA downregulated genes happening by chance
phyper(length(intersect(open_ensembl, rna_down)) - 1, length(unique(open_ensembl)), nrow(result) - length(unique(open_ensembl)), length(unique(rna_down)), lower.tail = FALSE)

## [1] 0.001046141

atac_rna <- data.frame(
  Pathway = character(),
  Length = integer(),
  Genes = character(),
  stringsAsFactors = FALSE
)

## convert between gene ENSEMBL, SYMBOL and ENTREZID
map_function <- function(x, column, keytype) {
  return(unname(mapIds(
    org.Hs.eg.db,
    keys = x,
    column = column,
    keytype = keytype
  )))
}

## extract core enrichment genes from GO or KEGG GSEA result, ENTREZ ids separated by "/"
cord_enrichment <- function(x, y) {
  cord <- x$core_enrichment[x$Description == y]
  cord <- unlist(strsplit(cord, "/"))
  return(cord)
}

## convert overlap result into one row of data frame
to_data_frame <- function(x, i, intersection) {
  x <- rbind(x, data.frame(
    "Pathway" = i,
    "Length" = length(intersection),
    "Genes" = paste(intersection, collapse = ", ")
  ))
  return(x)
}

## calculate overlap between pathway genes and RNA upregulated genes
atac_rna_overlap <- function(x, dataset, column, keytype, rna = rna_up, cord = FALSE) {
  for (i in x) {
    if (cord) {
      id <- cord_enrichment(dataset, i)
    } else {
      id <- dataset[[i]]
    }
    new_id <- map_function(id, column, keytype)
    intersection <- intersect(new_id, rna)
    atac_rna <- to_data_frame(atac_rna, i, intersection)
  }
  return(atac_rna)
}

atac_rna <- atac_rna_overlap(pathway_open_hallmark, hallmark_gene_sets, "ENSEMBL", "SYMBOL")
atac_rna <- atac_rna_overlap(pathway_open_go, gsea_go_df, "ENSEMBL", "ENTREZID", cord = TRUE)
atac_rna <- atac_rna_overlap(pathway_open_kegg, gsea_kegg_df, "ENSEMBL", "ENTREZID", cord = TRUE)

Just to further validate the results from the ATAC analysis, I performed RNA-seq differential expression analysis using the dataset from the TCGA-BRCA project. Before running the differential expression analysis, I removed all genes with counts less than 10 to reduce noise. Then DESeq2 was used to identify genes that are upregulated or downregulated in tumor samples compared to normal samples.

The RNA PCA plot shows two very distinct groups, which means the gene expression patterns are very different between the two conditions, which is what we expect to see. The tumor samples also appear to spread out a little more than the normal samples, which is similar to what we observed in the ATAC PCA plot. For the MA plot, there are a large number of genes significantly changed in both directions, meaning many genes are expressed differently between tumor samples and normal samples.

To compare the RNA results with the ATAC results, I looked at the intersection between ATAC open genes and both RNA upregulated genes and RNA downregulated genes. A hypergeometric test was then used to calculate the probability that these overlaps happened by chance. The probability is extremely low (3.74e-5 for ATAC open and RNA upregulated genes, and 0.00105 for ATAC open and RNA downregulated genes). Genes associated with open chromatin regions show more overlap with RNA upregulated genes than with RNA downregulated genes. This suggests that genes with increased accessibility in tumor samples are more likely to also have increased expression.

I also checked how many genes from the top five significantly enriched pathways for each test overlap with the RNA upregulated genes. The overlap results show that several pathways identified from the ATAC-seq GSEA analysis also contain a substantial number of RNA upregulated genes, which further supports the results.

References

NCBI Sequence Read Archive. BioProject: PRJNA1069018. Available at: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1069018

The Cancer Genome Atlas Research Network. Comprehensive molecular portraits of human breast tumours. Nature. 2012;490:61–70.