read in gene annotations

library('GenomicFeatures')
library('rtracklayer')

library(GenomicFeatures)
txdbfile <- '/n/projects/CompGenomics/Data/M.musculus.Ens_91.mm10.db'
txdb <- loadDb(txdbfile)

# get tss
tss.gr <- promoters(txdb,0,1,columns=c("tx_id", "tx_name","gene_id"))

# keep only the stuff on cannonical chromosomes
#tss.gr <- keepStandardChromosomes(tss.gr, pruning.mode="coarse")

Map oct4 summit to nearest TSS

# find nearest TSS
ntss.idx <- nearest(oct4.summits,tss.gr)

# keep summits that have a distance (not NA)
tss_anchored_summits.gr <- oct4.summits[!is.na(ntss.idx)]

# get a copy of the near TSS (those that map to a summit)
# in the order in which they map (see help on na.omit)
ntss.gr <- tss.gr[na.omit(ntss.idx)]

# transfer the strand of the nearest gene to the summit
strand(tss_anchored_summits.gr) <- strand(ntss.gr)

# make sure order of seqlevels is the same in both of these things
all(seqlevels(ntss.gr) == seqlevels(tss_anchored_summits.gr))

## [1] FALSE

# seqlevels have to be in same order for the operation below
# so put them in order
seqlevels(ntss.gr) <- sort(seqlevels(ntss.gr))
seqlevels(tss_anchored_summits.gr) <- sort(seqlevels(tss_anchored_summits.gr))

# find the upstream or downstream distance
tss_anchored_summits.gr$tss_nearest_distance <-
  distance(tss_anchored_summits.gr,ntss.gr)  *
  ifelse(tss_anchored_summits.gr < ntss.gr,-1,1) *
  ifelse('-' == strand(tss_anchored_summits.gr),-1,1)

# write data to the meta data setction of our summits
mcols(tss_anchored_summits.gr) <- cbind(mcols(tss_anchored_summits.gr),mcols(ntss.gr))

What are the distances between oct4 binding sites and nearest TSSs?

summary(tss_anchored_summits.gr$tss_nearest_distance)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -925038   -8496      -7    1042   10127  887272

ggplot(data=data.frame(dist=tss_anchored_summits.gr$tss_nearest_distance),
  aes(dist)) + geom_histogram() + xlim(-10000,10000)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 12160 rows containing non-finite values (stat_bin).

What is the distribution of peak scores?

## Examine histogram of oct4 peak scores
hist(oct4.summits$score, breaks=50)
abline(v=-10 * log10(0.01), col="red")

Let’s filter our data to reduce the number of peaks

## Filtering only summits within 10kb
is.good <- tss_anchored_summits.gr$tss_nearest_distance < abs(10^3)
oct4.topSet <- tss_anchored_summits.gr[is.good]

# how many were returned?
sum(is.good)

## [1] 14401

There are still quite a few. Let’s partition the data into 4 quartiles based on score.

# first let's examine a summary
summary(oct4.topSet$score)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.009   7.370  12.182  25.576  26.371 540.763

# then let's grab the value at the 90th percentile
quantile(oct4.topSet$score, 0.9)

##      90% 
## 59.16802

quantile(oct4.topSet$score, c(0.25,0.5,0.75))

##      25%      50%      75% 
##  7.37034 12.18200 26.37071

scoresByQ <- cut(oct4.topSet$score,
                     breaks = quantile (oct4.topSet$score, c (0, .25, .5, .75, 1)),
                     include.lowest = TRUE)

Using Factors in R, an Example

What is a factor? It’s simply a description of group memberships for our data.

If our data is categorized, or labelled in groups, R is good at utilizing this. Here’s a simple example of getting the average height of jockeys or ball players from a set of height measurements in inches.

# heights in inches of Jockeys or Basketball players
heights <- c(66,64,82,62,78,80,67)
playerType <- c("jock","jock","bask","jock","bask","bask","jock")

We can create a factor for the playerType variable, and use it to automatically group the data using an R function. In this case we can get the average height of each group using tapply()

# create the factor
playerType.f <- factor(playerType)
# apply the mean function
tapply(heights, playerType.f, mean)

##  bask  jock 
## 80.00 64.75

Splitting summits into quartiles by score

Back to our data set, once we have a factor representing a grouping label for each of our data points, we can split the data into those groups.

oct4ByQ <- split(oct4.topSet, scoresByQ)

The result of this is a list of GRanges objects (a GrangesList), and this is complicated. So rather than split our data this way, let’s simply split the index positions of our original object - that way we can use those to grab out pieces of our object at will.

# split the index of our object into groups
oct4ByQ.iv <- split(1:length(oct4.topSet), scoresByQ)

Take a look at the result, and see if this makes any more sense. The result is still a list, but the values in each section represent the elements of the oct4.topSet GRanges in each category. Let’s examine it.

# what are the names of the list elements?
names(oct4ByQ.iv)

## [1] "[2.01,7.37]" "(7.37,12.2]" "(12.2,26.4]" "(26.4,541]"

We see the score values that make up the boundaries of each category.

How many elements are in each category?
Examine a few elements of each.

# because it's a list, we can use lapply on it.
# How long is each group?
lapply(oct4ByQ.iv, length)

## $`[2.01,7.37]`
## [1] 3613
## 
## $`(7.37,12.2]`
## [1] 3594
## 
## $`(12.2,26.4]`
## [1] 3594
## 
## $`(26.4,541]`
## [1] 3600

# Examine the first few elements of each group
lapply(oct4ByQ.iv, head)

## $`[2.01,7.37]`
## [1]  5  9 15 18 22 23
## 
## $`(7.37,12.2]`
## [1]  2  6  7  8 11 12
## 
## $`(12.2,26.4]`
## [1] 19 27 28 29 30 32
## 
## $`(26.4,541]`
## [1]  1  3  4 10 13 14

Let’s extract 500 randomly chosen elements from each quantile. We reference the 1st list element, which contain 3613 items and choose 500 at random.

Summarize the scores for the first quantile:

# take a random sample of 500 ids
q1.ids <- sample(oct4ByQ.iv[[1]], 500)

summary(oct4.topSet$score[q1.ids])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.311   4.775   5.694   5.635   6.616   7.370

Pick random sets for each quantile.

q1.idx <- sample(oct4ByQ.iv[[1]], 500)
q2.idx <- sample(oct4ByQ.iv[[2]], 500)
q3.idx <- sample(oct4ByQ.iv[[3]], 500)
q4.idx <- sample(oct4ByQ.iv[[4]], 500)

boxplot(list(q1=oct4.topSet$score[q1.idx],
             q2=oct4.topSet$score[q2.idx],
             q3=oct4.topSet$score[q3.idx],
             q4=oct4.topSet$score[q4.idx]), main="Peak Scores by Quartile")

Do peak widths differ by quantile?

uh, they’re summits - they’re all the same width. We have to load in the actual peaks to test this. Skip this chunk.

boxplot(list(q1=width(oct4.topSet[q1.idx]),
             q2=width(oct4.topSet[q2.idx]),
             q3=width(oct4.topSet[q3.idx]),
             q4=width(oct4.topSet[q4.idx])), main="Peak Width by Quartile")

Get the genes associated with each quantile

# get the gene ids
goi.q1 <- unique(as.character(oct4.topSet$gene_id[q1.idx]))
goi.q2 <- unique(as.character(oct4.topSet$gene_id[q2.idx]))
goi.q3 <- unique(as.character(oct4.topSet$gene_id[q3.idx]))
goi.q4 <- unique(as.character(oct4.topSet$gene_id[q4.idx]))

# grab the coordinates for each gene TSS
q1.tss.gr <- tss.gr[match(goi.q1, as.character(tss.gr$gene_id))]
q2.tss.gr <- tss.gr[match(goi.q2, as.character(tss.gr$gene_id))]
q3.tss.gr <- tss.gr[match(goi.q3, as.character(tss.gr$gene_id))]
q4.tss.gr <- tss.gr[match(goi.q4, as.character(tss.gr$gene_id))]

# resize them 
q1.tss.gr <- resize(q1.tss.gr, width=1000, fix="center")
q2.tss.gr <- resize(q2.tss.gr, width=1000, fix="center")
q3.tss.gr <- resize(q3.tss.gr, width=1000, fix="center")
q4.tss.gr <- resize(q4.tss.gr, width=1000, fix="center")

# make sure they didn't extend past boundaries
q1.tss.gr <- trim(q1.tss.gr)
q2.tss.gr <- trim(q2.tss.gr)
q3.tss.gr <- trim(q3.tss.gr)
q4.tss.gr <- trim(q4.tss.gr) 

Now we can extrcat out the RNA Polymerase II values and make metaplots for each quartile.

# set the location of our RNA polII bigwWig file
bwFile <- "/n/projects/CompGenomics/Data/GSE115436/R1_WT_ES_RPB1.bw"

polII.q1 <- grangesMatrix(q1.tss.gr, bwFile)
polII.q2 <- grangesMatrix(q2.tss.gr, bwFile)
polII.q3 <- grangesMatrix(q3.tss.gr, bwFile)
polII.q4 <- grangesMatrix(q4.tss.gr, bwFile)

Let’s make a plot and see what it looks like!

plot(colMeans(polII.q4), type="l", main="RNA Polymerase II by OCT4 quartile", ylab="mean RPM", xaxt="n")
lines(colMeans(polII.q3), col="blue")
lines(colMeans(polII.q2), col="green")
lines(colMeans(polII.q1), col="red")
# add and axis
axis(1, at=c(0,500,1000), labels=c("-500","0","+500"))
# add a vertical line at TSS
abline(v=500, col="red", lty="dashed", lwd=2)

legend("topright", c("q1", "q2","q3","q4"), col=c("red","green","blue","black"), lwd=2)