Summarizing Overlaps of Query Regions with Genomic Annotation Features
Querying the annotation file
Finding targeted gene types
To find out the distribution of the query regions across gene types:
biotype_col <- grep('gene_biotype', colnames(overlaps), value = T)
df <- overlaps[,length(unique(queryIndex)), by = biotype_col]
colnames(df) <- c("feature", "count")
df$percent <- round(df$count / length(queryRegions) * 100, 1)
df <- df[order(count, decreasing = TRUE)]
ggplot2::ggplot(df, aes(x = reorder(feature, -percent), y = percent)) +
geom_bar(stat = 'identity', aes(fill = feature)) +
geom_label(aes(y = percent + 0.5), label = df$count) +
labs(x = 'transcript feature', y = paste0('percent overlap (n = ', length(queryRegions), ')')) +
theme_bw(base_size = 14) +
theme(axis.text.x = element_text(angle = 90))
Extending the annotation feature space
GTF files contain some annotation features (e.g. exons, transcripts) that are usually explicitly defined, however, some transcript features such as introns, exon-intron boundaries, promoter regions are only implicitly defined. Such implicit features can be extracted from a GTF file using makeTxDb family of functions from the GenomicFeatures library.
First we create a list of GRanges objects, where each list element contains all the available coordinates of transcript features such as transcripts, exons, introns, 5’/3’ UTRs, exon-intron boundaries, and promoter regions.
Plotting overlap counts between query regions and transcript features
To have a global overview of the distribution of query regions across gene features, we can use the summarizeQueryRegions function. If a given query region does not overlap with any of the given coordinates of the transcript features, it is categorized under NoFeatures.
summary <- summarizeQueryRegions(queryRegions = queryRegions,
txdbFeatures = txdbFeatures)
df <- data.frame(summary)
df$percent <- round((df$count / length(queryRegions)), 3) * 100
df$feature <- rownames(df)
ggplot2::ggplot(df, aes(x = reorder(feature, -percent), y = percent)) +
geom_bar(stat = 'identity', aes(fill = feature)) +
geom_label(aes(y = percent + 3), label = df$count) +
labs(x = 'transcript feature', y = paste0('percent overlap (n = ', length(queryRegions), ')')) +
theme_bw(base_size = 14) +
theme(axis.text.x = element_text(angle = 90))
Obtaining a table of overlap counts between query regions and genes
To find out which genes overlap with how many queries and categorise overlaps by transcript features; we use getTargetedGenesTable function, which returns a data.frame object.
dt <- getTargetedGenesTable(queryRegions = queryRegions,
txdbFeatures = txdbFeatures)
dt <- dt[order(transcripts, decreasing = TRUE)]
knitr::kable(dt[1:10,])| tx_name | transcripts | exons | promoters | fiveUTRs | introns | cds | threeUTRs |
|---|---|---|---|---|---|---|---|
| ENST00000317713 | 33 | 28 | 0 | 0 | 5 | 24 | 4 |
| ENST00000361689 | 33 | 28 | 0 | 0 | 5 | 24 | 4 |
| ENST00000372915 | 33 | 28 | 0 | 0 | 5 | 24 | 4 |
| ENST00000539005 | 33 | 28 | 0 | 0 | 5 | 24 | 4 |
| ENST00000545844 | 33 | 28 | 0 | 0 | 5 | 24 | 4 |
| ENST00000564288 | 33 | 28 | 0 | 0 | 5 | 24 | 4 |
| ENST00000567887 | 33 | 28 | 0 | 0 | 5 | 24 | 4 |
| ENST00000372925 | 28 | 23 | 0 | 0 | 5 | 19 | 4 |
| ENST00000289893 | 27 | 22 | 0 | 0 | 5 | 18 | 4 |
| ENST00000367142 | 14 | 14 | 0 | 0 | 0 | 3 | 12 |
Profiling the coverage of query regions across transcript features
Coverage profile of query regions at feature boundaries
It may be useful to look at the distribution of query regions at the boundaries of transcript features. For instance, it may be important to see the relative signal at transcript ends (transcription start sites versus transcription end sites). Or, it may be important to see how the signal is distributed at exon boundaries, which may give an idea about the regulation of the transcript. Here we demonstrate how to get such signal distributions at transcription start/end sites. The same approach can be done for any other collection of transcript features (exons, introns, promoters, UTRs etc.)
cvgF <- getFeatureBoundaryCoverage(queryRegions = queryRegions,
featureCoords = txdbFeatures$transcripts,
flankSize = 1000,
boundaryType = 'fiveprime',
sampleN = 10000)
cvgT <- getFeatureBoundaryCoverage(queryRegions = queryRegions,
featureCoords = txdbFeatures$transcripts,
flankSize = 1000,
boundaryType = 'threeprime',
sampleN = 10000)
cvgF$boundary <- 'fiveprime'
cvgT$boundary <- 'threeprime'
df <- rbind(cvgF, cvgT)
ggplot2::ggplot(df, aes(x = bases, y = meanCoverage)) +
geom_ribbon(fill = 'lightgreen',
aes(ymin = meanCoverage - standardError * 1.96,
ymax = meanCoverage + standardError * 1.96)) +
geom_line(color = 'black') +
facet_grid( ~ boundary) + theme_bw(base_size = 14) 
Coverage profile of query regions for all transcript features
Coverage profiles can be obtained for a single type of transcript feature or a list of transcript features. Here we demonstrate how to get coverage profile of query regions across all available transcript features. It might be a good idea to use sampleN parameter to randomly downsample the target regions to speed up the calculations.
cvgList <- calculateCoverageProfileList(queryRegions = queryRegions,
targetRegionsList = txdbFeatures,
sampleN = 10000)
ggplot2::ggplot(cvgList, aes(x = bins, y = meanCoverage)) +
geom_ribbon(fill = 'lightgreen',
aes(ymin = meanCoverage - standardError * 1.96,
ymax = meanCoverage + standardError * 1.96)) +
geom_line(color = 'black') + theme_bw(base_size = 14) +
facet_wrap( ~ feature, ncol = 3) 
