Contaminants

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:

• rRNA: You can set rRNA = “silva” to download the Silva database (~ 2GB file)
• ncRNA: If the gtf does not have Non-coding RNAs, they can be extracted by setting ncRNA = “auto”, it will then check if the organism exists in the NONCODE database and automatically download them for you if they exists.
• tRNA: Genomes without defined tRNAs are rare, and no web API to access tRNAs for all common experimental organisms exists. So in that rare case, you must manually download and add the sequences from tRNAs using tRNA scan or similar databases (then add this file to the getGenomeAndAnnotation call using argument tRNA = “path/to/downloaded/tRNAs.fasta”)

Create 5’ UTR annotation from CAGE

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

Create Pseudo 5’ UTR annotation

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

Fixing malformed gtf/gff

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.

Trimming data

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.

Aligning the data

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

• input.dir.rna: directory with raw fastq files (or user pre-trimmed files)
• output.dir.rna: output directory for bam files
• index: the STAR index from previous step
• paired.end: “yes” in this case, or “no” if single end.
• steps: steps wanted from trimming, depletion and alignment wanted: (a string: which steps to do? (default: “tr-ge”, write “all” to get all: “tr-co-ge”) tr: trimming (only for unix / WSL), co: deplete contaminants included, ph: phix depletion, rR: rrna depletion, nc: ncrna depletion, tR: trna depletion, ge: genome alignment) Write your wanted steps, seperated by “-”. Order does not matter. To just do trim and alignment to genome write “tr-ge”
• adapter.sequence “auto”, or if you know add it, usually more secure with manual. Presets are “illumina”, “small_RNA” and “nextera”.
• max.cpus How many cpus maximum to use
• trim.front How many bases to trim front. Only if you believe there are low quality reads in front.
• min.length minimum length of reads that pass to the bam file.

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.

UMIs

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 arguments

A 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 duplicated reads

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!

RAM usage warnings

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:

• 20 GB max ram usage during genome generation
• 2 SA sparse (suffix array should be sparse, this will give slower mapping)

index <- STAR.index(annotation, max.ram = 20, SAsparse = 2)

Systems that restricts max open files

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 systems

Restart STAR from crashed step

STAR 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 tr

This would save you the time to re-trim, or update the path of input files to the fastq files of the trimmed folder.