The RNA Centric Analysis System Report

Bora Uyar, Dilmurat Yusuf, Ricardo Wurmus, Altuna Akalin

2018-03-23

library(RCAS)

Warning: In this vignette, due to space limitations, we demonstrate the functions of RCAS using static images. In order to see how an interactive report from RCAS looks see RCAS::runReport().

For the most up-to-date functionality, usage and installation instructions, and example outputs, see our github repository here.

Introduction

RCAS is an automated system that provides dynamic genome annotations for custom input files that contain transcriptomic regions. Such transcriptomic regions could be, for instance, peak regions detected by CLIP-Seq analysis that detect protein-RNA interactions, RNA modifications (alias the epitranscriptome), CAGE-tag locations, or any other collection of target regions at the level of the transcriptome.

RCAS is designed as a reporting tool for the functional analysis of RNA-binding sites detected by high-throughput experiments. It takes as input a BED format file containing the genomic coordinates of the RNA binding sites and a GTF file that contains the genomic annotation features usually provided by publicly available databases such as Ensembl and UCSC. RCAS performs overlap operations between the genomic coordinates of the RNA binding sites and the genomic annotation features and produces in-depth annotation summaries such as the distribution of binding sites with respect to transcript features (exons, introns, 5’/3’ UTR regions, exon-intron boundaries, promoter regions, and whole transcripts). Moreover, by detecting the collection of targeted transcripts, RCAS can carry out functional annotation tables for enriched gene sets (annotated by the Molecular Signatures Database) and GO terms. As one of the most important questions that arise during protein-RNA interaction analysis; RCAS has a module for detecting sequence motifs enriched in the targeted regions of the transcriptome. The final report of RCAS consists of high-quality dynamic figures and tables, which are readily applicable for publications or other academic usage.

Data input

RCAS minimally requires as input a BED file and a GTF file. The BED file should contain coordinates/intervals of transcriptomic regions which are located via transcriptomics methods such as Clip-Seq. The GTF file should provide reference annotation. The recommended source of GTF files is the ENSEMBLE database.

For this vignette, in order to demonstrate RCAS functionality, we use sample BED and GTF data that are built-in the RCAS library, which can be imported using a common R function: data(). To import custom BED and GTF files, the user should execute two RCAS functions called importBed() and importGtf().

Importing sample data

Importing custom data

To use importBed() and importGtf(), the user should provide file paths to the respective BED file and GTF file. To reduce memory usage and time consumption, we advise the user to set sampleN=10000 to avoid huge input of intervals.

Summarizing Overlaps of Query Regions with Genomic Annotation Features

Querying the annotation file

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.

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.

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.)

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.

Motif Analysis using motifRG

Calculating enriched motifs

With the RCAS package, a motif analysis is also possible. RCAS uses motifRG library to find enriched motifs among the query regions.

## AAGGAG 9.456508e-07 
## Skip pattern  ATTTTT 
##  Refine  AAGGAG 12.26311 : 11.85714 11.31425 11.95031 10.97143 12.26579 TRUE 482 147 467 139 
## New motif:  AAGGAG 
## match range  624 
## [1] "Rescore"
## [1] "Finished Rescore"
## TGGAGA 1.020551e-06 
## Skip pattern  TTTCTT TTTTCT TTTTTA TTTTAA 
##  Refine  TGGAGA 12.23307 : 12.35601 13.16996 11.58157 12.41869 12.47023 TRUE 555 197 534 190 
## New motif:  TGGAGA

motif analysis: getting motif summary statistics

A summary table from the motif analysis results can be obtained

patterns scores fgHits bgHits fgSeq bgSeq ratio fgFrac bgFrac
AAGGAG 12.3 482 147 467 139 3.3 0.0467 0.0139
TGGAGA 12.3 558 197 536 190 2.8 0.0536 0.0190

GO term analysis

Biological processes enriched among targeted genes

RCAS can perform GO term enrichment analysis to find out enriched functions in genes that overlap the query regions. Below is demonstrated how to get biological processes terms (‘BP’) enriched in the genes that overlap query regions and the top 10 GO terms with most fold change increase relative to the background are provided.

Term Significant Expected bonferroni foldEnrichment
GO:0043487 regulation of RNA stability 16 4.92 0.0006496 3.25
GO:0006402 mRNA catabolic process 26 8.28 0.0000005 3.14
GO:0006401 RNA catabolic process 28 9.17 0.0000003 3.05
GO:0006403 RNA localization 15 4.92 0.0054880 3.05
GO:0051168 nuclear export 15 4.92 0.0054880 3.05
GO:0061013 regulation of mRNA catabolic process 14 4.70 0.0179200 2.98
GO:0015931 nucleobase-containing compound transport 15 5.14 0.0123200 2.92
GO:0006611 protein export from nucleus 13 4.47 0.0537600 2.91
GO:1903311 regulation of mRNA metabolic process 20 7.16 0.0010752 2.79
GO:0034655 nucleobase-containing compound catabolic… 28 10.96 0.0001008 2.55

Gene Set Enrichment Analysis

MSIGDB gene sets enriched among targeted genes

RCAS can use gene sets from Molecular Signatures Database and calculate gene set enrichment analysis (GSEA) to find out which gene sets are enriched among the genes targeted by the query regions.

Below we demonstrate a GSEA case using randomly generated gene sets (in order not to breach MSIGDB licence agreement) that are provided as built-in data in RCAS. The actual MSIGDB gene set annotations must be downloaded by the user from the MSIGDB website. RCAS provides functions to parse the annotations (RCAS::parseMsigdb) and map them to other species via orthology (RCAS::createOrthologousGeneSetList) to enable GSEA on other species such as mouse and fly.

treatment treatmentSize expectedInTreatment fisherPVal BH bonferroni foldEnrichment
randomGeneSet52 9 415 3.5 0.0215522 0.6129037 1 2.57
randomGeneSet16 8 415 3.2 0.0318513 0.6129037 1 2.50
randomGeneSet99 7 415 3.0 0.0577985 0.6129037 1 2.33
randomGeneSet87 10 415 4.4 0.0269477 0.6129037 1 2.27
randomGeneSet95 10 415 4.4 0.0269477 0.6129037 1 2.27
randomGeneSet42 9 415 4.1 0.0391768 0.6129037 1 2.20
randomGeneSet8 7 415 3.2 0.0696890 0.6129037 1 2.19
randomGeneSet53 7 415 3.2 0.0696890 0.6129037 1 2.19
randomGeneSet11 10 415 4.6 0.0323445 0.6129037 1 2.17
randomGeneSet13 4 415 2.0 0.1903736 0.6129037 1 2.00

RCAS also provides functions to map the MSIGDB annotations from human to fly and mouse.

#parse human annotations
refGeneSets <- parseMsigdb(filePath = <path to MSIGDB annotation file>)

#Map the gene sets to other species using orthologous relationships of genes between
#the reference genome (human) and the target genome (e.g. mouse)
orthGeneSets <- createOrthologousGeneSetList(referenceGeneSetList = refGeneSets, 
                                                refGenomeVersion = 'hg19', 
                                                targetGenomeVersion = 'mm9')
#the mapped gene sets can be used for GSEA analysis using the runGSEA command.

Generating a full report

The users can use the runReport() function to generate full custom reports including all the analysis modules described above. There are four main parts of the analysis report.

By default, runReport() function aims to run all four modules, while the user can turn off these individual modules.

Below are example commands to generate reports using these functionalities.

A test run for human

runReport()

A custom run for human

runReport( queryFilePath = 'input.BED',
            gffFilePath = 'annotation.gtf',
            msigdbFilePath = 'human_msigdb.gmt')

To turn off certain modules of the report

runReport( queryFilePath = 'input.BED',
            gffFilePath = 'annotation.gtf',
            msigdbFilePath = 'human_msigdb.gmt',
            motifAnalysis = FALSE,
            goAnalysis = FALSE )

To run the pipeline for species other than human

If the msigdb module is needed, the msigdbFilePath must be set to the MSIGDB annotations for ‘human’. MSIGDB datasets for other species will be calculated in the background using the createOrthologousMsigdbDataset function

runReport( queryFilePath = 'input.mm9.BED',
            gffFilePath = 'annotation.mm9.gtf',
            msigdbFilePath = 'human_msigdb.gmt',
            genomeVersion = 'mm9' )

To turn off verbose output and progress bars

runReport(quiet = TRUE)

Printing raw data generated by the runReport function

One may be interested in printing the raw data used to make the plots and tables in the HTML report output of runReport function. Such tables could be used for meta-analysis of multiple analysis results. In order to activate this function, printProcessedTables argument must be set to TRUE.

runReport(printProcessedTables = TRUE)

Acknowledgements

RCAS is developed in the group of Altuna Akalin (head of the Scientific Bioinformatics Platform) by Bora Uyar (Bioinformatics Scientist), Dilmurat Yusuf (Bioinformatics Scientist) and Ricardo Wurmus (System Administrator) at the Berlin Institute of Medical Systems Biology (BIMSB) at the Max-Delbrueck-Center for Molecular Medicine (MDC) in Berlin.

RCAS is developed as a bioinformatics service as part of the RNA Bioinformatics Center, which is one of the eight centers of the German Network for Bioinformatics Infrastructure (de.NBI).