target.R
In target: Predict Combined Function of Transcription Factors

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE
)

## ----load_libraries-----------------------------------------------------------
# load reguired libraries
library(target)
library(GenomicRanges)

## ----load_data----------------------------------------------------------------
# load peaks and transcripts data
data("real_peaks")
data("real_transcripts")

## ----associated_peaks---------------------------------------------------------
# get associated peaks
ap <- associated_peaks(real_peaks, real_transcripts, 'ID')
ap

## ----direct_tragets-----------------------------------------------------------
# get direct targets
dt <- direct_targets(real_peaks, real_transcripts, 'ID', 't')
dt

## ----show_relation,fig.align='center',fig.width=7,fig.height=3----------------
par(mfrow = c(1, 3))

# show peak distance vs score
plot(ap$distance, ap$peak_score,
     pch = 19, cex = .5, 
     xlab = 'Peak Distance', ylab = 'Peak Score')
abline(v = 0, lty = 2, col = 'gray')

# show gene stat vs score
plot(dt$stat, dt$score,
     pch = 19, cex = .5, 
     xlim = c(-35, 35),
     xlab = 'Gene t-stats', ylab = 'Gene Score')
abline(v = 0, lty = 2, col = 'gray')

# show gene regulatory potential ecdf
groups <- c('Down', 'None', 'Up')
colors <- c('darkgreen', 'gray', 'darkred')

fold_change <- cut(dt$logFC,
                   breaks = c(min(dt$logFC), -.5, .5, max(dt$logFC)),
                   labels = groups)

plot_predictions(dt$score_rank,
                 fold_change,
                 colors,
                 groups,
                 xlab = 'Gene Regulatory Potential Rank',
                 ylab = 'ECDF')

## ----test_groups--------------------------------------------------------------
# test up-regulated transcripts are not random
test_predictions(dt$score_rank,
                 group = fold_change,
                 compare = c('Up', 'None'),
                 alternative = 'greater')

# test down-regulated transcripts are not random
test_predictions(dt$score_rank,
                 group = fold_change,
                 compare = c('Down', 'None'),
                 alternative = 'greater')

## ----top_gene_transcript------------------------------------------------------
# show the top regulated transcript, gene name and its associated peaks
top_trans <- unique(dt$ID[dt$rank == min(dt$rank)])
top_trans

unique(dt$name2[dt$ID == top_trans])
unique(ap$peak_name[ap$assigned_region == top_trans])

## ----session_info-------------------------------------------------------------
sessionInfo()