Using factR
Fursham Hamid
2024-05-01
1 Introduction
Many eukaryotic genes give rise to multiple RNA isoforms, increasing the protein-coding capacity of the genome and extending the range of post-transcriptional regulation possibilities. High-throughput sequencing is often used to deduce repertoires of transcripts expressed in specific biological samples by aligning the data to genomic sequences and assembling the alignments into transcript architectures. This typically outputs Gene Transfer Format (GTF) files describing newly identified transcripts as sets of exonic coordinates, but lacking the information about their coding sequences (CDSs; also known as ORFs) and possible biological functions.
To this end, we developed a package for functional annotation of custom-assembled transcriptomes in R (factR). factR predicts CDSs for novel RNA isoforms using a reference-guided process and then determines domain organisation of the protein products and possible susceptibility of transcripts to nonsense-mediated decay (NMD; a pathway destabilizing mRNAs with premature translation termination codons). factR also provides supporting tools for matching new transcripts to “official” gene IDs, visualizing transcript architectures and annotating alternatively spliced segments.
2 Getting started
2.1 Installing factR
To install factR, enter the following commands in your R environment:
2.2 Materials needed
2.2.1 Custom-assembled transcriptome
factR requires a custom transcriptome file in the GTF format as an input, and we provide the following three sample custom transcriptome files that can be used to test factR tools:
- bulk_merged_sample.gtf.gz assembled using the HISAT2-StringtTie2 pipeline (Pertea et al. 2016) from bulk RNA-seq data for mouse embryonic stem cells treated with the NMD inhibitor cycloheximide (CHX) or left untreated as a control.
- sc_merged_sample.gtf.gz assembled using the HISAT2-StringtTie2 pipeline from single-cell RNA-seq data for glutamatergic neurons, GABAergic neurons, astrocytes or endothelial cells from the mouse visual cortex (Tasic et al. 2018). A simplified version of this GTF containing the first 500 genes on chr15 is available as part of factR.
- lr_merged_sample.gtf.gz assembled using minimap2 (Li 2018) and StringtTie2 (Kovaka et al. 2019) from brain-specific long-read Oxford Nanopore RNA sequencing data (Sessegolo et al. 2019).
Method to download the above GTF files is described in the “Importing and inspecting GTF data” section below.
Users may alternatively prepare their own GTF files making sure that these contain both gene_id and transcript_id attributes.
2.2.2 Reference transcriptome
Several factR functions require reference transcriptome files (GTF or GFF3) as a guide. Such files can be accessed from within the R environment using e.g. the R package AnnotationHub or downloaded from an external database such as GENCODE or Ensembl. Both possibilities are described in the “Updating gene info” section below.
2.2.3 Genome
factR needs genomic DNA sequence to predict CDSs. Users may obtain genome files using e.g. BSgenome or AnnotationHub or download them from an external database such as GENCODE or Ensembl. We describe this in more detail in the “Constructing CDS information” section below.
3 Using factR
Load factR to the R environment as follows:
3.1 Importing and inspecting GTF data
factR handles transcriptome information in the form of GenomicRanges objects containing genomic interval data and relevant metadata. To create such an object from a GTF file, we use the importGTF function:
gtf <- system.file("extdata", "sc_merged_sample.gtf.gz", package = "factR")
custom.gtf <- importGTF(gtf) The imported GTF file is stored as a GenomicRanges object.
Contents of the object can be examined using head:
head(custom.gtf)
#> GRanges object with 6 ranges and 9 metadata columns:
#> seqnames ranges strand | source type score
#> <Rle> <IRanges> <Rle> | <factor> <factor> <numeric>
#> [1] chr15 3180731-3180944 * | StringTie transcript 1000
#> [2] chr15 3180731-3180944 * | StringTie exon 1000
#> [3] chr15 3217391-3219698 - | StringTie transcript 1000
#> [4] chr15 3217391-3219698 - | StringTie exon 1000
#> [5] chr15 3268547-3277274 + | StringTie transcript 1000
#> [6] chr15 3268547-3268768 + | StringTie exon 1000
#> phase gene_id transcript_id exon_number gene_name
#> <integer> <character> <character> <character> <character>
#> [1] <NA> MSTRG.14523 MSTRG.14523.1 <NA> <NA>
#> [2] <NA> MSTRG.14523 MSTRG.14523.1 1 <NA>
#> [3] <NA> MSTRG.14524 ENSMUST00000227053.1 <NA> Gm7962
#> [4] <NA> MSTRG.14524 ENSMUST00000227053.1 1 Gm7962
#> [5] <NA> MSTRG.14525 ENSMUST00000160787.8 <NA> Selenop
#> [6] <NA> MSTRG.14525 ENSMUST00000160787.8 1 Selenop
#> ref_gene_id
#> <character>
#> [1] <NA>
#> [2] <NA>
#> [3] ENSMUSG00000114999.1
#> [4] ENSMUSG00000114999.1
#> [5] ENSMUSG00000064373.12
#> [6] ENSMUSG00000064373.12
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengthsIn addition to genomic coordinates (seqnames and ranges), a typical GTF file contains metadata describing the feature type (e.g. transcript or exon), transcript IDs and some information on their parental genes.
Use the following command to calculate the total number of transcripts in the input transcriptome:
Note that none of the transcripts in the custom.gtf object contain CDS information:
3.2 Plotting transcript structures
Users may visualize specific sets of transcripts using viewTranscripts. For example, the following will plot transcripts from the Zfr gene encoding a conserved zinc finger-containing RNA-binding protein with known neuronal functions :
StringTie2, a popular transcript assembler used to generate our custom GTF, typically assigns arbitrary names to newly identified transcripts (e.g. “MSTRG.x.y”) and uses the same prefix for their gene IDs (e.g. “MSTRG.x”). Incidentally, this is why the output above contains 10 previously known Zfr transcripts but lacks any novel entries.
3.3 Updating gene info
factR can update gene metadata in custom transcriptome objects using a reference annotation as guide. Below, we describe two alternative ways to obtain mouse reference transcriptome data from GENCODE.
- Using AnnotationHub package
# query database for mouse gencode basic annotation
library(AnnotationHub)
ah <- AnnotationHub()
query(ah, c('Mus musculus', 'gencode', 'gff'))
#> AnnotationHub with 9 records
#> # snapshotDate(): 2024-04-29
#> # $dataprovider: Gencode
#> # $species: Mus musculus
#> # $rdataclass: GRanges
#> # additional mcols(): taxonomyid, genome, description,
#> # coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
#> # rdatapath, sourceurl, sourcetype
#> # retrieve records with, e.g., 'object[["AH49545"]]'
#>
#> title
#> AH49545 | gencode.vM6.2wayconspseudos.gff3.gz
#> AH49546 | gencode.vM6.annotation.gff3.gz
#> AH49547 | gencode.vM6.basic.annotation.gff3.gz
#> AH49548 | gencode.vM6.chr_patch_hapl_scaff.annotation.gff3.gz
#> AH49549 | gencode.vM6.chr_patch_hapl_scaff.basic.annotation.gff3.gz
#> AH49550 | gencode.vM6.long_noncoding_RNAs.gff3.gz
#> AH49551 | gencode.vM6.polyAs.gff3.gz
#> AH49552 | gencode.vM6.primary_assembly.annotation.gff3.gz
#> AH49553 | gencode.vM6.tRNAs.gff3.gz
# Download full annotation
ref.gtf <- ah[['AH49546']]
#> loading from cache- Downloading from a GENCODE database (ref.gtf generated from this method will be used for the remaining of the workflow)
tmp <- tempfile()
download.file("ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M25/gencode.vM25.annotation.gtf.gz",
destfile = tmp)
ref.gtf <- importGTF(tmp) When choosing a reference, users should consider using one with the same chromosome naming style (e.g. “chr1” or “1”). Alternatively, the styles can be matched using the matchChromosomes function (see help(matchChromosomes) for more detail).
Once the ref.gtf object is ready, novel transcripts in custom.gtf can be assigned to “official” gene IDs, whenever possible, using matchGeneInfo. By default, this function matches the query (custom.gtf) to the reference (ref.gtf) by finding overlapping coordinates:
# matching gene metadata
custom_matched_1.gtf <- matchGeneInfo(custom.gtf, ref.gtf)
#> Number of mismatched gene_ids found: 500
#> ---> Attempting to match gene_ids by finding overlapping coordinates...
#> ---> 386 gene_id matched
#> Total gene_ids corrected: 386
#> Remaining number of mismatched gene_ids: 114To tune the performance of matchGeneInfo, we provide additional arguments specifying the name of columns containing “primary” and potentially “secondary” gene IDs from the query (custom.gtf). For further information, see help(matchGeneInfo).
# matching gene metadata
custom_matched_2.gtf <- matchGeneInfo(custom.gtf, ref.gtf,
primary_gene_id = "gene_id",
secondary_gene_id = "ref_gene_id")
#> Number of mismatched gene_ids found: 500
#> -> Attempting to correct gene ids by replacing gene_id with ref_gene_id...
#> -> 212 gene_ids matched
#> --> Attempting to match ensembl gene_ids...
#> --> All ensembl gene ids have been matched
#> ---> Attempting to match gene_ids by finding overlapping coordinates...
#> ---> 174 gene_id matched
#> Total gene_ids corrected: 386
#> Remaining number of mismatched gene_ids: 114Note that custom.gtf updated by matchGeneInfo now contains 10 known and 7 newly Zfr assembled transcripts:
3.4 Shortlisting novel transcripts
As seen in the above example, custom transcriptomes typically combine both new and previously annotated transcripts. To select only newly predicted transcripts, run the following:
custom_new.gtf <- subsetNewTranscripts(custom_matched_2.gtf, ref.gtf)
#> Removing transcripts with exact exon coordinates
viewTranscripts(custom_new.gtf,"Zfr") 
This will subset custom.gtf transcripts with distinct exonic coordinates compared to ref.gtf and will store these transcripts in the custom_new.gtf object. Some custom-built transcripts may differ from their reference counterparts by having different start or/and end coordinates, with otherwise similar exon-intron structure. To shortlist novel transcripts with distinct intronic coordinates only, simply set the “refine.by” argument to “intron”:
custom_new.gtf <- subsetNewTranscripts(custom_matched_2.gtf, ref.gtf, refine.by = "intron")
#> Removing transcripts with exact exon coordinates
#> Removing transcripts with exact intron coordinates
viewTranscripts(custom_new.gtf, "Zfr") 
We will use the custom_new.gtf object in the rest of the workflow.
3.5 Constructing CDS information
Functional annotation of newly assembled transcripts in factR begins by deducing their protein-coding sequences (CDSs). To search for putative CDSs, factR requires a genome sequence file, which can be obtained using R packages or downloaded from online databases (e.g. UCSC, GENCODE or Ensembl). Three alternative ways to retrieve mouse genomic sequence are described below.
- Using BSgenome (Mmusculus generated from this method will be used for the remaining of the workflow) This package supports most sequenced genomes. Mouse mm10 sequence can be downloaded as follows:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("BSgenome.Mmusculus.UCSC.mm10") and loaded into R environment:
- Using AnnotationHub
library(AnnotationHub)
ah <- AnnotationHub()
query(ah, c("mm10","2bit"))
#> AnnotationHub with 1 record
#> # snapshotDate(): 2024-04-29
#> # names(): AH14005
#> # $dataprovider: UCSC
#> # $species: Mus musculus
#> # $rdataclass: TwoBitFile
#> # $rdatadateadded: 2014-12-15
#> # $title: mm10.2bit
#> # $description: UCSC 2 bit file for mm10
#> # $taxonomyid: 10090
#> # $genome: mm10
#> # $sourcetype: TwoBit
#> # $sourceurl: http://hgdownload.cse.ucsc.edu/goldenpath/mm10/bigZips/mm10.2bit
#> # $sourcesize: NA
#> # $tags: c("2bit", "UCSC", "genome")
#> # retrieve record with 'object[["AH14005"]]'- Downloading from a database (e.g. GENCODE)
tmp <- tempfile()
download.file("ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M25/GRCm38.primary_assembly.genome.fa.gz",
tmp)
Mmusculus <- importFASTA(paste0(tmp, "/GRCm38.primary_assembly.genome.fa.gz")) Once the genome sequence object is ready, factR can predict CDSs using its buildCDS() function and reference transcriptome data as a guide. buildCDS() first generates a database of previously annotated ATGs and uses this information to search for a potential translation start sites in query transcripts. buildCDS() then deduces the CDS and appends its coordinates to the custom transcriptome object. Let’s run this function for our novel transcripts:
custom_new_CDS.gtf <- buildCDS(custom_new.gtf, ref.gtf, Mmusculus)
#> Searching for reference mRNAs in query
#> No reference mRNAs found
#> Building database of annotated ATG codons
#> Selecting best ATG start codon for remaining transcripts and determining open-reading frame
#> 114 new CDSs constructed
#>
#> Summary: Out of 380 transcripts in `custom_new.gtf`,
#> 114 transcript CDSs were builtNote that the novel Zfr transcripts have been updated with information about likely CDSs (dark blue) and untranslated regions (light blue) :