ORFik 1.27.0
Welcome to the ORFik package.
ORFik is an R package for analysis of transcript and translation features through manipulation of sequence data and NGS data.
This vignette will walk you through how to how to download annotation and align data with ORFik.
Here we will show a full example of aligning RNA-seq from yeast using the SacCer3 genome.
Working with NGS data and genome annotations, results ++, can be hassle to organize and to remember where everything is stored, ORFik tries to simplify this job for you by creating a fixed set of how files are organized. Using ORFiks organizing system is optional (you can specify all the paths manually if you want to), but using it will most likely make you code faster.
Usually fastq files, processed data (bam files) and annotation (genome references) are stored in separate locations, for better clarity and re-usability of files.
Given a single parrent folder (for example “~/Bio_data/”), where all your Biological analysis is stored, ORFik uses a folder structure of separating these 3 types:
This is defined in the config() function, and only is needed to run once (config is stored for later experiments), so all other experiments will reuse the parent file structure. ORFik will try to automatically create suitable locations if the following function is not run;
library(ORFik) # This package
# Here the default values are shown:
where_to_save_config <- config_file() # Bioc Cache directory default
parent_folder <- "~/bio_data/" # Change if you want it somewhere else
fastq.dir <- file.path(parent_folder, "raw_data")
bam.dir <- file.path(parent_folder, "processed_data")
reference.dir <- file.path(parent_folder, "references")
exp.dir <- file.path(base.dir, "ORFik_experiments/")
config.save(where_to_save_config,
fastq.dir, bam.dir, reference.dir, reference.dir)Either way if you made a custom config or let ORFik do it for you in next step, you can check what it is by now running config()
config()We have created the general folder structure, we must now create the folder structure for a specific experiment, like our yeast RNA-seq experiment from Chalmers et al.
conf <- config.exper(experiment = "CHALMERS_Yeast", # short experiment info: here I used author + species
assembly = "Yeast_SacCer3", # Reference folder
type = c("RNA-seq")) # fastq and bam type
confconf will show where the specific experiment will be stored.
We will now download the NGS data, if you have in-lab data, you don’t need this step, since you already have access to the fastq files.
On the other hand, if you want to use published data, you need to download it. I here show what would work for the paired end RNA-seq experiment SRP012047.
ORFik comes with a read archive run downloader (Archives supported are: SRA, ERA, DRA or GSE). We will now show how to get metadata and data.
info <- download.SRA.metadata("SRP012047", outdir = conf["fastq RNA-seq"])
# Let's take 2 first runs in this experiment:
info <- info[1:2,]Now let’s download the subset, using the metadata table (ORFik will automatically extract SRR numbers and download):
# Subset
# 18 MB, ~ 40 sec download time ->
download.SRA(info, conf["fastq RNA-seq"], subset = 50000)
# Or you could do the full libraries
# 1.6 GB, ~ 100 sec download time
# (Full size files have faster MB/s download speed than subsets) ->
# download.SRA(info, conf["fastq RNA-seq"])We now have the RNA-seq run, separated into 2 files stored in the raw data fastq folder, there are 2 files since this is paired end data. We could for ease also just have specified the SRR number directly in download.SRA, but then we get no meta-data csv file, which is handy for auto-detection of paired end data, the organism name etc. This is shown below:
organism <- info$ScientificName[1]
is_paired_end <- all(info$LibraryLayout == "PAIRED")There are two ways to use annotation:
To download the genome annotation we use the getGenomeAndAnnotation function. We need to decide 3 things:
The uncompressed toplevel file of human genome is > 70 GB, so for large genomes (like for human) you should usually set to primary_assembly (3 GB for human, no chromosome haplotypes). For small organisms like yeast, they don’t have a primary assembly, so use “toplevel”.
annotation <- getGenomeAndAnnotation(
organism = organism,
output.dir = conf["ref"],
assembly_type = "toplevel" # Usually leave this to "auto"
)The function will also create a txdb object to speed up loading of gtf annotation, and index your genome to a .fai file (fasta index), for faster access to genomic regions.
If you rerun this function after is has completed, it will not re-download, only output the correct object paths, this makes it easy to rerun the script (for other experiments of same species), when you have some steps already completed before.
If you you want to remove contaminants: phix, non coding RNAs, ribosomal RNAs, or tRNAs, also specify these in the function. By default it will download phix from refseq and the other contaminants are usually defined within the genome annotation of the species, so they are extracted from the .gtf file.
annotation <- getGenomeAndAnnotation(
organism = organism,
genome = TRUE, GTF = TRUE,
phix = TRUE, ncRNA = TRUE, tRNA = TRUE, rRNA = TRUE,
output.dir = conf["ref"],
assembly_type = "toplevel"
)Note that some species does not have well annotated rRNAs, tRNAs etc in their gtf files. For this ORFik contains some clever tricks:
If you are not downloading annotation through ORFik, we need to create the Txdb object (Bioconductor optimized annotation format)
This will save you a lot of the common warnings and problems downstream: - correct seqlevels - correct seqlengths - correct organism assigned - faster loading
gtf <- "/path/to/local/annotation.gtf"
genome <- "/path/to/local/genome.fasta"
makeTxdbFromGenome(gtf, genome, organism = "Saccharomyces cerevisiae")The txdb is saved in same directory as the gtf, with an extension “.db”
It should be noted that Bioconductor BSGenome packages are not currently supported as annotation for ORFik alignments, since they lack some of the functionality needed for ORFik. They might be supported in the future and can already be used for all other parts of ORFik.
If you do not have defined UTRs in your annotation, this will give you considerable less quality results (no results using leaders can be created). So optimial is always to find an annotation with UTRs. If that is not possible ORFik contains two ways to create 5’ UTR annotation, 3’ UTR annotation will be included soon.
If you have CAGE libraries, ORFik can in a simple way create 5’ UTRs for your annotation:
txdb <- "/path/to/local/annotation.gtf.db"
cage <- "/path/to/CAGE.ofst" # Can be bed, wig etc (ofst is fastest)
reassigned_txdb <- assignTSSByCage(txdb, cage)
AnnotationDbi::saveDb("/path/to/local/annotation_cage.db")If you do not have CAGE, you can also just specify some fixed 5’ UTR lengths, for all CDS transcripts that does not have a 5’ UTR.
txdb <- "/path/to/local/annotation.gtf.db"
reassigned_txdb <- assignTSSByCage(txdb, cage = NULL, pseudoLength = 100)
AnnotationDbi::saveDb(reassigned_txdb,
"/path/to/local/annotation_pseudo_leaders.db")
leaders <- loadRegion(reassigned_txdb, "leaders")If you have a non reference-genome gff, your gff might be malformed. Your gtf/gff needs to follow the gtf/gff standardizations: specifications.
Most importantly it need the 8 first coordinate columns correctly defined, and the 9th attribute column needs at least gene_id and transcript_id! Else it will be useless in most of R analysis in ORFik. Even better if it also has gene_biotype column defined, this will make you able to load for example all long-non-coding RNA coordinates directly.
ORFik uses the STAR aligner, which is splice aware and fast and shown to usually map the highest amount of reads. The sad thing is that it only works on unix systems (Linux or Mac).
For windows users: installing and running ORFik on Windows subsystem for Linux (WSL) will work, if you want to align data with ORFik.
To align the data we need two steps, the indexing of the genome step and the alignment to the genome step.
To index the genome just give it the annotation output from previous step. This will also make an index for each of the depletion steps like phix, if you specified them in the earlier step.
index <- STAR.index(annotation)If you run this function again after index exists in current file location, it will not re-index, but just output the correct object paths. Do remake = TRUE if you want to re-index.
ORFik uses the fastp for trimming reads, this also only works on unix (Linux or Mac OS). If you are on windows, or you want to trim the reads yourself, just run the trimming and give the folder with the trimmed reads as input in next step.
Most Illumina sequencing contains 3’ adapters. This is because let’s say you sequence 100 bases, but the actual RNA read is only 28 nt long, you still need 72 bases as fillers, to not screw up that all should be 100 bases long. The great thing with fastp is that it has auto detection and removal of adapters, if you check out the resulting files you will see fastp has auto removed the Illumina adapters. In rare cases some libraries are have a double set of adapters, or low quality adapters, in those cases run fastqc (or another 3’ adapter analyzer) to see what it missed.
The final rule is that for Ribo-seq, your trimmed reads should be no longer than 55 nt, else STAR (using the default local alignment) will not be able to find the 28nt of the read that align (STAR must by default align ~ 50% of the read). For RNA-seq, adapter trimming is not really that crucial, as local alignment will find the hit anyway.
Now let’s see what we need as inputs for the trimming / alignment pipeline: We need usually 9 arguments (more are possible if you need them):
alignment <-
STAR.align.folder(conf["fastq RNA-seq"], conf["bam RNA-seq"], index,
paired.end = is_paired_end,
steps = "tr-ge", # (trim needed: adapters found, then genome)
adapter.sequence = "auto",
max.cpus = 30, trim.front = 3, min.length = 20)If you used the fastp (tr step), you will get a pre-alignment QC report. Just like FastQC in html format. You will also get a MultiQC report from STAR runs made by ORFik for you.
Some of the newer NGS libraries contains Unique molecular identifiers. If you need to store the UMI in the read header, fastp supports this. A direct wrapper for this in STAR.align.folder is currently not implemented yet, but will be soon! For now you can call fastp directly through a system call:
dir.create("~/UMIandTrimmed")
system(paste(install.fastp(),
"-i", "~/sample1.fastq.gz",
"-o", "~/UMIandTrimmed/sample1.UMIandTrimmed.fastq.gz"
"--adapter_sequence=AGATCGGAAGAGC",
"--umi",
"--umi_loc", "read1",
"--umi_len", 12))
# Read fastp documentation for info about umi argumentsA speed up, if UMIs are not needed is also to just trim the UMIs away. Using the “trim.front” argument in STAR.align.folder. This can be done when sequencing is not done on very low input samples and sequencing <= 80 million reads. A good intro can be found here: dnatech.genomecenter.ucdavis.edu
Collapsing of duplicated reads is supported (all identical reads merged, with a replicate counter in the header), to do this, first only trim the reads to remove adapters. If your reads are < 50 bases long on average, it is a good idea to collapse duplicated reads (much faster alignment and processing). Note in this case the reads are paired end, so then collapsing is very dangerous. But for single end reads, you could do this:
alignment <-
STAR.align.folder(conf["fastq RNA-seq"], conf["bam RNA-seq"], index,
paired.end = is_paired_end,
steps = "tr", # (Only trim)
adapter.sequence = "auto",
max.cpus = 30, trim.front = 3, min.length = 20)
trimmed_files <- file.path(conf["bam RNA-seq"], "trim")
# Check that trimmed reads are average < 50 bases, else collapsing makes little sense for speedup.
trim_table <- ORFik:::trimming.table(trimmed_files)
stopifnot(all(trim_table$trim_mean_length < 50))
trimmed_fastq <- list.files(path = trimmed_files, pattern = "fastq", full.names = TRUE)
ORFik::collapse.fastq(files = trimmed_fastq, outdir = file.path(conf["bam Ribo-seq"], "collapsed") # Collapse all files in trim folder
# Then align using collapsed reads
alignment <-
STAR.align.folder(trimmed_files, conf["bam RNA-seq"], index,
paired.end = is_paired_end,
steps = "ge", # (Only genome alignment from collapsed)
max.cpus = 30, min.length = 20)Usually, RNA-seq is > 50 bases, so only do this for short reads like Ribo-seq!
Running STAR is a heavy job for most machines, so you may experience an error if your machine is not good enough or you have a non-default OS config (Aligning human genome on 32GB memory computer will make it useless for any other task until alignment is done, so if you have access to a multi-core server with more memory, you should use it).
We will now go through some possible problems and how to fix them.
STAR is very memory hungry, therefor if you want to index and align large genomes like human (3 GB usually), you might need to adjust 2 parameters. If you have less than 40 GB free memory (32 GB might work), adjust this during indexing:
index <- STAR.index(annotation, max.ram = 20, SAsparse = 2)STAR can use a lot of threads, each threads makes many several files open. Some systems have a restriction on how many files you can have open. This is found by doing “ulimit -Hn” in the terminal. If STAR crashes from this error, you need either increase amount of open files allowed (requires root access) or decrease the amount of cores used by STAR (does not require root access):
STAR.align.folder(conf["fastq RNA-seq"], conf["bam RNA-seq"], index,
max.cpus = 4) # Reduce cores to 4 usually works for most systemsSTAR can sometimes crash, maybe your disc becomes full or you specified a wrong parameter etc. ORFik has some clever tricks to continue alignment from where it crashed. Let’s say the trimming is done and it fails during genome alignment. Then of course if you rerun the function as it was, it would re-trim the raw fastq-reads. To continue using the completed trim files instead and only do genome alignment, do:
STAR.align.folder(conf["fastq RNA-seq"], conf["bam RNA-seq"], index,
max.cpus = 4, steps = "tr-ge", resume = "ge") # Resume ge using completed trThis would save you the time to re-trim, or update the path of input files to the fastq files of the trimmed folder.
The next step is to make an ORFik experiment (check out the ORFik experiment vignette if you are unfamiliar with this class).
If you renamed the files during the fastq step properly, this will make the ORFik experiment able to guess correctly what the data is. If there are replicates etc. If you did not, you can still rename the bam files now.
We can now easily make an ORFik experiment from the data we have:
txdb_file <- paste0(annotation["gtf"], ".db") # Get txdb file, not raw gtf
fa <- annotation["genome"]
create.experiment(exper = "yeast_exp_RNA",
dir = paste0(conf["bam RNA-seq"], "/aligned/"),
txdb = txdb_file,
fa = fa,
organism = organism,
viewTemplate = FALSE,
pairedEndBam = is_paired_end # True/False per bam file
)The files is now saved to default directory, which is: saveDir = “~/Bio_data/ORFik_experiments/”
df <- read.experiment("yeast_exp_RNA")If you are not happy with the libtype, stage, replicates and so on for the file, you can edit the ORFik experiment in R (recreate experiment, wanted slots) or edit in Libre Office, Excel or another spreadsheet viewer.
See ?QCreport for details of what you will get as output
QCreport(df)Now you have an experiment, but bam files are big and slow to load. Let’s convert to some faster formats.
If you want optimzed format identical to bam file (contains cigar information), use .ofst. (Fastest, not readable in IGV) (ofst files are made when running ORFikQC)
remove.experiments(df) # Remove loaded libraries
convertLibs(df, type = "ofst")If you want peaks only, use wig files (Fast, readable in IGV)
remove.experiments(df)
convertLibs(df, type = "wig")As an example of how to load the data to R in the optimized format .ofst.
This will output the libraries to the environment specified, default .GlobalEnv (the default R environment). The files are named from the experiment table RNA_1_WT, RNA_1_treated etc.
remove.experiments(df)
outputLibs(df, type = "ofst")After you have run QCreport you will have count tables of peaks over the mrna’s, 5’ UTRs, CDS and 3’ UTRs.
Let’s do an example to find the ratio between fpkm of between the CDS and mRNAs transcript regions.
mrna <- countTable(df, region = "mrna", type = "fpkm")
cds <- countTable(df, region = "cds", type = "fpkm")
ratio <- cds / mrnaWe now have a ratio of fpkm values between CDS and mrna.
You can now continue to the Ribo-seq pipeline to see a more advanced example.