Lecture 11

Set up

library(data.table)
library(tidyverse)
library(lattice)
library(gridExtra)
library(rtracklayer)
library(idr2d)
library(patchwork)
library(reactable)
library(ggplot2)
library(DESeq2)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
library(plotly)
library(ChIPseeker)
library(annotate)
library(gprofiler2)
library(clusterProfiler)
options(ChIPseeker.ignore_1st_exon = TRUE)
options(ChIPseeker.ignore_1st_intron = TRUE)
options(ChIPseeker.ignore_downstream = TRUE)
options(ChIPseeker.ignore_promoter_subcategory = TRUE)
# set your working directory
setwd("~/Documents/HGEN_663/extra/lec11")

Import

# Read in and process the feature counts from /home/hgen_share/lec11 as previously shown
mat <- read.table() %>%
  column_to_rownames()

DESeq2

mdat <- data.frame()
mdat$cond <- c("mut","mut","rescue","rescue") 
mdat <- mdat %>%
  column_to_rownames("kind")
dds <- DESeqDataSetFromMatrix(countData = ?,
                              colData = ?,
                              design = ?)
dds <- DESeq(?)
# get results dataframe
res <- results(dds,contrast = c("cond","mut","rescue")) %>%
  data.frame() %>%
    na.omit() %>%
  rownames_to_column("Geneid")

PCA

The code to run a PCA is given. However, you will first need to run variance stabilizing transformation as previously shown.

vst <- ?
# subset data - get the top 500 most variable peaks
rv <- rowVars(assay(vst))
select <- order(rv, decreasing=T)[seq_len(min(500,length(rv)))]
# get pc scores
pc <- prcomp(t(assay(vst)[select,]))
condition <- mdat$cond
scores <- data.frame(pc$x, condition)
percentage <- round(pc$sdev / sum(pc$sdev) * 100, 2)
percentage <- paste( colnames(scores), "(", paste( as.character(percentage), "%", ")", sep="") )
# plot
p <- plot_ly(scores,x=scores$PC1,y=scores$PC2,text=rownames(scores),mode="markers",color=factor(condition),marker=list(size=11),colors = "Dark2")
p <- layout(p,title="",   xaxis = list(title = percentage[1]),
            yaxis = list(title = percentage[2]))
p

Volcano plot

volc <- function(r, x, y, ttl) {
  d <- as.data.frame(r) %>%
    dplyr::rename(x = !!x, y = !!y) %>%
    mutate(kind = case_when(abs(x) > 1 & y < .05 ~ '|log2FC|>1,FDR<0.05',
                            abs(x) < 1 & y < .05 ~ '|log2FC|<1,FDR<0.05',
                            y > .05 ~ 'FDR > 0.05'),
           y = -log10(y))
  ct <- na.omit(d) %>% 
    dplyr::filter(kind == '|log2FC|>1,FDR<0.05') %>%
    mutate(up = x > 0) %>%
    dplyr::count(up) %>%
    mutate(x = ifelse(up, Inf, -Inf),
           y = Inf,
           h = as.numeric(up))
  ggplot(d, aes(x, y, color = kind)) +
    geom_vline(xintercept = c(-1, 1), linetype = 'dashed') +
    geom_hline(yintercept = -log10(0.05), linetype = 'dashed') +
    geom_point(alpha = .5) +
    geom_label(aes(x = x, y = y, label = n, hjust = h),
              vjust = 1, data = ct, inherit.aes = F) +
    scale_y_continuous() +
    scale_color_manual(values = c('royalblue', 'red2', 'grey30')) +
    labs(x = x, y = "-log10(FDR)", title = ttl) +
    theme_bw()
}
volc(res,'log2FoldChange','padj','mutant (H3K36me3 -) versus rescue (H36K36me3 +)')

Peaks

Processing

# filter 'res' object for open and closed ATAC-seq peaks as well as separate the genomic intervals for each peak into 'chr', 'start' and 'end'
open_peaks <- res[res$log2FoldChange > 1 & res$padj < 0.05,] %>%
  tidyr::separate(Geneid,c("chr","start","end"))
closed_peaks <- res[res$log2FoldChange < -1 & res$padj < 0.05,] %>%
  tidyr::separate(Geneid,c("chr","start","end"))
# get rid of NAs and non-numeric values
value <-  as.numeric( as.character(open_peaks$start) ) # get the numbers
noninteger <- value %% 1 != 0   # see if there's a fractional part
noninteger <- noninteger | is.na(noninteger)  # get rid of the NA's 
open_peaks <- open_peaks[!noninteger,]
# make granges object
open_peaks_gr <- makeGRangesFromDataFrame(open_peaks)

# Do the same for closed peaks
# get rid of NAs and non-numeric values
value <-  as.numeric( as.character(closed_peaks$start) ) # get the numbers
noninteger <- value %% 1 != 0   # see if there's a fractional part
noninteger <- noninteger | is.na(noninteger)  # get rid of the NA's 
closed_peaks <- closed_peaks[!noninteger,]
# make granges object
closed_peaks_gr <- makeGRangesFromDataFrame(closed_peaks)

# Lastly,we can do the same for the original pooled peaks
pooled <- res %>%
  tidyr::separate(Geneid,c("chr","start","end"))
# get rid of NAs and non-numeric values
value <-  as.numeric( as.character(pooled$start) ) # get the numbers
noninteger <- value %% 1 != 0   # see if there's a fractional part
noninteger <- noninteger | is.na(noninteger)  # get rid of the NA's 
pooled <- pooled[!noninteger,]
# make granges object
pooled_peaks_gr <- makeGRangesFromDataFrame(pooled)

Distribution

We can use ChIPSeeker to annotate our peaks. Below, I’ve shown how to annotate the open peaks for you. Annotate the closed and pooled peaks as well as the open peaks. Combine all three types of peaks and compare their distribution.

# load known genes for hg19
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

# annotate the peaks
ann <- annotatePeak(open_peaks_gr,TxDb = txdb,tssRegion = c(-3000, 3000))

## >> preparing features information...      2023-03-20 01:07:06 
## >> identifying nearest features...        2023-03-20 01:07:07 
## >> calculating distance from peak to TSS...   2023-03-20 01:07:07 
## >> assigning genomic annotation...        2023-03-20 01:07:07 
## >> assigning chromosome lengths           2023-03-20 01:07:28 
## >> done...                    2023-03-20 01:07:28

# plot the peak distribution
plotAnnoBar(ann)

# plot the percentage per each annotation, here it's shown for open peaks
plotAnnoPie(ann)

# plot distribution of TF-binding loci relative to TSS. Since we observed almost most of the open peaks to be in distal intergenic regions, it’s not surprising that a large percentage of sites are far away from the TSS
plotDistToTSS(ann)

Pathway enrichment analysis

We can investigate whether genes nearest each peak are enriched for a pathway. Considering this mutation is associated with kidney cancer, some of these enriched pathways makes sense, such as ‘renin secretion’ and ‘calcium signaling pathway’. Perform KEGG enrichment analysis for the closed peaks.

# you will need a background list, which will be all the genes nearest the peaks from your pooled peak list
pooled <- ann$pooled_peaks %>% as.data.frame()
bg <- pooled$geneId %>%
  as.character() %>%
  unique()
# genes for open peaks
open_peak_df <- ann$open_peaks %>%
  as.data.frame() 
entrezids <- open_peak_df$geneId %>%
  as.character() %>%
  unique()
#
ekegg <- enrichKEGG(gene = entrezids,
                 organism = 'hsa',
                 pvalueCutoff = 0.05,
                 pAdjustMethod = "fdr",
                 universe = bg)
dotplot(ekegg)

Extra

Split ATAC-seq peaks according to their annotation: promoter, intron or intergenic

promoter <- ann$open_peaks %>% as.data.frame() %>%
  dplyr::filter(str_detect(annotation, "Promoter")) %>%
  dplyr::select(1:3)
intergenic <- ann$open_peaks %>% as.data.frame() %>%
  dplyr::filter(str_detect(annotation, "Distal Intergenic")) %>%
  dplyr::select(1:3)
intron <- ann$open_peaks %>% as.data.frame() %>%
  dplyr::filter(str_detect(annotation, "Intron")) %>%
  dplyr::select(1:3)
export.bed(object = promoter,con="open_peaks.promoter.bed")
export.bed(object = intergenic,con="open_peaks.intergenic.bed")
export.bed(object = intron,con="open_peaks.intron.bed")

Convert peaks into sequence strings. Essentially, going from .bed files to .fasta files

# load required libraries
library(BSgenome.Hsapiens.UCSC.hg19)
library(seqinr)

# load genome
Hsapiens <- BSgenome.Hsapiens.UCSC.hg19

# extract sequence from 'granges' object
seq <- getSeq(Hsapiens,open_peaks_gr)

# you need to create an 'id' for each sequence, here we are naming them peak1, peak2, etc.
names <- paste0("peak",1:8228)

# export fasta file
write.fasta(seq,names = names,file.out = "open_peaks.fasta")