Latest revision as of 16:29, 7 October 2022

Introduction

This document describes the analysis of RNA-Seq data using OLego (for alignment), gapless (for inference of transcript structure), and countit (for quantification of gene expression and alternative splicing). gapless and countit are collectively named Quantas. The pipeline has been tested quite extensively, but the document is still at its relatively early stage. Any comments are welcome and can be directed to Chaolin Zhang.

Versions

• v1.0.9 (12-01-2015), current
  • minor change of output format of expression summary files of individual samples
  • minor bug fixes
• v1.0.8 (03-10-2015)
  • improved efficiency and memory usage in junction tag counting
  • minor bug fixes
• v1.0.7 (08-14-2014)
  • enhanced support for different stranded library types
  • included additional filtering metric in differential expression analysis
• v1.0.6 (04-22-2014)
  • improved memory usage/speed in splicing quantification
  • included reports of read counts of individual junctions for tandem cassette (taca)
  • fixed a rare bug fix in FDR calculation for diff splicing
  • included a glm for differential splicing test
• v1.0.4 (11-10-2013)
  • added several scripts into the package
• v1.0.3 (06-09-2013)
  • bug fixes
• v1.0.2 (05-31-2013)
  • add all scripts required for polyA-seq analysis; will be documented later
• v1.0.1 ( 5-23-2013 )
  • expanded methods to test differential splicing
  • included script to test differential expression

• v1.0.0 ( 04-8-2013 )
  • The initial Public release

Important assumptions

• This document is for paired-end (PE) mRNA-seq libraries that have NO strand information, and the two mates are on the opposite strands (the standard Illumina protocol).
• It is IMPORTANT to note that read1 and read2 in the same pair must have the same and unique read name.
• The focus of this document is to quantify splicing and gene expression change of annotated transcripts. We are in the process of extending it to analysis of novel transcripts and alternative splicing events.

Single-end (SE) mRNA-seq libraries or strand-specific SE or PE RNA-seq libraries can also be analyzed by the same packages, but the parameters are slightly different.

Software installation

OLego installation and preparation

The program we use for read mapping is OLego, a program we developed recently and made available at http://zhanglab.c2b2.columbia.edu/index.php/OLego. There you can download the source codes or precompiled binaries and access the more detailed documentation.

The program takes

• read file(s) in either FASTA or FASTQ format
• the indexed genome
• (optionally) a database of known exon junctions

For the reference genome, it is recommended that one should remove random chromosomes as well as the mitochondria genome. The program required to build the index is included in the package. Assuming one has concatenated all chromosomes in a FASTA file mm10.fa (note that the name of each chromosome is chrN instead of N to be consistent with the other library files). The following command will index the genome.

olegoindex -a bwtsw mm10.fa

This will produce several files with a prefix mm10.fa. Move all the output files to a directory such as genome_index/olego/mm10

A comprehensive database of exon junctions (mm10) can be downloaded from http://zhanglab.c2b2.columbia.edu/data/OLego/mm10.intron.hmr.bed.gz. Download and decompress this file.

See the OLego website for more details.

Quantas installation and preparation

Quantas includes a package to infer transcript structure from paired-end RNA-seq data (gapless), and a package to count RNA-seq reads for each alternative splicing isoform (countit). It also depends on a perl library czplib, which contains various functions for genomic/bioinformatic analysis.

Through anaconda

Below are the installation instructions for Quantas, CZPLIB (dependency) and OLego (dependency) through Anaconda.

conda config --env --append channels conda-forge
conda config --env --append channels bioconda

conda install --yes -c chaolinzhanglab quantas

• Install the perl library files czplib

Download from github

Decompress it and move it to a place you like

$tar zxvf czplib.v1.0.x.tgz
$mv czplib /usr/local/lib

Add the library path to the environment variable, so perl can find it.

PERL5LIB=/usr/local/lib/czplib

One might need to install Math:CDF package, if not already. It is available at http://search.cpan.org/~callahan/Math-CDF-0.1/CDF.pm

• Download gapless & countit [Quantas]

Download from github

$tar zxvf Quantas.v1.x.x.tgz

One also needs to download annotation files used by Quantas at

• mouse (mm10): http://zhanglab.c2b2.columbia.edu/data/Quantas/data/mm10.tgz (29Mb)
• human (hg19): http://zhanglab.c2b2.columbia.edu/data/Quantas/data/hg19.tgz (37M)

(A previous version of the two packages missed some files for gene expression analysis, which is corrected now - 02/06/2014)

Mapping RNA-seq reads to the reference genome and exon junctions

Mapping

Assuming one has a paired-end sequencing lane of 100 nt x2 reads in sample.r1.fa and sample.r2.fa for the first and second read mates, respectively, derived from mouse.

olego -v -t 16 -r olego_dir/models/mm.cfg -j mm10.intron.hmr.bed -o sample.r1.sam genome_index/olego/mm10.fa sample.r1.fa
# do the mapping with 15nt seed with 1nt overlap, allowing a total of 4 mismatches or indels (default), with 16 cpu cores, output to sample.r1.sam

olego -v -t 16 -r olego_dir/models/mm.cfg -j mm10.intron.hmr.bed -o sample.r2.sam genome_index/olego/mm10.fa sample.r2.fa
# do the same thing for the other file

No matter how many threads (-t 16) are used, the memory footprint is in general < 4Gb. That said, the exact amount of memory also depends on the seed size, the repetitive nature of reads and the exon junction database.

This will map reads to the exons, known junctions as well as novel junctions. By default OLego searches with perfect matches in seeds of 15 nt with one 1-nt overlap. The maximal mismatches or indels along the whole read allowed in a match is determined by the read length(4 nt for 100 nt illumina reads). For known junctions, OLego requires 5 or more nt overlap on either side of the junction. For novel junctions, OLego requires 8 nt overlap, and considers only GT/AG splice sites, to improve the accuracy.

Convert SAM to BED file (single-end data ONLY)

Important note: this should be ONLY done for single end reads. For paired end reads, gapless will take sam files and output BED files.

After the alignment step, the SAM files are converted into BED files using the script sam2bed.pl included in the OLego package.

perl sam2bed.pl --uniq -v sample.sam sample.uniq.bed

This will keep only reads that are unambiguously mapped to the reference genome (single hits), which are identified by the tag “XT:A:U” in the SAM file.

Inferring transcript structure between paired end reads (paired-end data ONLY)

Paired-ends reads could be located in different exons and since the reads are relatively short, the transcript structure between the sequenced ends could be ambiguous when alternative splicing occurs. The missing information, leveraging on the size constraints of each cDNA fragment and prior isoform abundance estimated from directly mapped junction reads, is inferred using a simple Bayes model.

perl gapless/gapless_huge_file.pl -v -sam -uniq --split-size 10000000 -isoform mm10.exon.trio.hmr.nr.bed -E 400 -big --print-singleton -o gapless_out sample.r1.sam sample.r2.sam

This script takes the sam files of read1 and read2 as input (-sam). Here, again it is IMPORTANT to note that read1 and read2 in the same pair must have the same read name and order.

-uniq means only uniquely mapped reads (single hits) are kept. They are identified by the tag “XT:A:U” in the SAM file. -E 400 means large exons with 400 nt or more are used to estimate the distribution of insert cDNA fragment size in the library (adaptor sequences are not included in calculation here).

Important: from v1.0.7, an additional option --library-type is introduced to indicate different library types. This is to replace the original option --ss.

unstranded:  unstranded libraries assuming the two reads are on the opposite strand
stranded-as: read1 is antisense and read2 is sense
stranded-sa: read1 is sense and read2 is antisense.

This step will try to combine read1 and read2 in the same pair if possible (i.e., they are possibly sampled from the same isoform in our database), and output two files in the gapless_out dir, a bed file pair.gapless.bed with combined reads or fragments, and a text file size_dist.txt with the estimated distribution of insert cDNA fragment size in the library. When alternative splicing is observed between a pair, each possible isoform, or fragment, is assigned a probability score, so that one fragment can potentially have multiple lines. This score is recorded in the 5th column of the BED file.

Quantify alternative splicing

Run countit

First, one needs to prepare a configuration file cz_annotation.loc that specifies the location of AS database files. These files can be downloaded from http://zhanglab.c2b2.columbia.edu/data/Quantas/data/mm10.tgz, if not already.

mm10    as  annotation/Mm.seq.all.cass.chrom.can.bed    cass
mm10    as  annotation/Mm.seq.all.taca.chrom.can.bed    taca
mm10    as  annotation/Mm.seq.all.alt5.chrom.can.bed    alt5
mm10    as  annotation/Mm.seq.all.alt3.chrom.can.bed    alt3
mm10    as  annotation/Mm.seq.all.mutx.chrom.can.bed    mutx
mm10    as  annotation/Mm.seq.all.iret.chrom.can.bed    iret

A template of this file (mm10.conf) is included in the archive above. The third column of this file has to be modified to point to the location of the other AS annotation files in the archive.

perl countit/summarize_splicing_wrapper.pl -c ./cache -v -big -weight -conf mm10.conf -dbkey mm10 -cass -taca -alt5 -alt3 -mutx -iret gapless_out/pair.gapless.bed countit_out

Note that the -c option specifies temporary output dir, which should have sufficient space. The -weight option is specified because each fragment from gapless now has a probability score, and this option should NOT be used for single end RNA-Seq data.

In the output directory, there is a file for each AS event (e.g., cass.count.txt). In these files, each line is an alternative splicing event in comparison of two isoforms. The first 10 columns in output file for each type of AS are the same:

chrom chromStart chromEnd name score strand type isoformIDs isoform1Tags isoform2Tags

The score column above is reserved, and the first 6 columns provide coordinate information so that one can quickly go to the region in UCSC genome browser (see visualization below).

The last three columns are the two isoforms compared, the read count for isoform1 and the read count for isoform2, which are the most critical for quantification of splicing level.

Each output file has extra columns, which are different depending on the type of AS events, and also important for filtering in downstream analysis, as summarized below:

cass: exonTags inclusionJuncction1Tags inclusionJunction2Tags skippingJunctionTags
taca: exonTags inclusionJuncctionTags skippingJunctionTags
alt5: altSSDistance exonTags proximalJunctionTags distalJunctionTags
alt3: altSSDistance exonTags proximalJunctionTags distalJunctionTags
mutx: 5'ExonTags 5'ExonJunction1Tags 5'ExonJunction2Tags 3'ExonTags 3'ExonJunction1Tags 3'ExonJunction2Tags
iret: retainedIntronTags junctionTags

In the command line above, you can also choose to run only certain types of AS events (e.g. -cass for cassette exons).

Summarize alternative splicing results and perform statistical tests

After one runs OLego, gapless and countit for each sample, the final step of AS analysis is to test differential splicing between two groups of samples. For this purpose one needs to prepare a configuration file to specify which sample belongs to each group.

#dataset.group.conf
sample1_countit_out<tab>group1
sample2_countit_out<tab>group1
sample3_countit_out<tab>group1
sample4_countit_out<tab>group2
sample5_countit_out<tab>group2
sample6_countit_out<tab>group2

The first column of the file above is the path to the output directory of countit. One can then run the following line for cassette exons (and similarly for other types of AS events):

perl countit/test_splicing_diff.pl -type cass -v --min-cov 20 --id2gene2symbol annotation/Mm.seq.all.AS.chrom.can.id2gene2symbol dataset.group.conf dataset.diff.txt

This script will pool the number of reads for each isoform from each of the individual samples in the same group. It will then perform Fisher's exact test using the total number of reads supporting each isoform, calculate exon inclusion level and the change between the two groups, etc.

Finally, it will filter AS events and consider those events with sufficient read coverage for multiple test correction (using the Benjamini method).

The option (--min-cov 20) means the following requirement:

junction_in + junction_skip≥20 AND junction_group1+junction_group2≥20

The output of the file has the following columns:

column 1: gene, Entrez gene id and gene symbol

columns 2-9: chrom   chromStart  chromEnd    name    score   strand  type    isoformIDs: the same as described in the output of countit  
column 10: coverage    I_g1(group1) I_g2(group2) dI_g1_vs_g2 pvalue  FDR

In general, we require FDR<0.05 (or 0.01), |dI|>0.1 (or 0.2) to call significant splicing changes.

A caveat of the Fisher exact test described above is that it does not take into account the within-group variation. Particular caution is warranted for samples with a substantial amount of heterogeneity (e.g. see visualization below). We are working on this to address the issue.

Quantify gene expression

This is also done with the countit package.

perl countit/summarize_expression_wrapper.pl -big --cache ../cache -exon mm10.exon.uniq.core.bed -e2g mm10.exon.uniq.core.id2gene2symbol -weight -v sample_out/pair.gapless.bed sample.expr.txt

Again, the -weight option is specified because each fragment from gapless now has a probability score. This option should not be used for single end RNA-seq data.

The output is like this

#gene_id    gene_symbol tag_num exon_len    RPKM

Here exon_len is the total exonic length of the gene, and tag_num is the total number of exonic reads (or more precisely fragments) in the gene. RPKM is always calculated using total number of exonic tags summed over all genes in the output file as a denominator, to avoid the complication of rRNA reads, etc.

This should be done for all samples.

Differential gene expression analysis can be performed with different tools, such as edgeR (http://www.bioconductor.org/packages/2.12/bioc/html/edgeR.html). The package includes a script for this purpose. One needs to prepare a configuration file to specify which sample belongs to each group.

#expr.diff.conf
sample1.expr.txt<tab>group1
sample2.expr.txt<tab>group1
sample3.expr.txt<tab>group1
sample4.expr.txt<tab>group2
sample5.expr.txt<tab>group2
sample6.expr.txt<tab>group2

#this script will invoke R and requires edgeR to be installed
perl countit/test_expr_diff.pl -MAplot -v expr.diff.conf group1_vs_group2.expr.diff.txt

We typically call significant differential expression by requiring FDR<0.01 (or 0.001), |log2FC|>1 (or 0.6 or 2), and max(group1RPKM, group2RPKM) in the top half.

Generate bedGraph files for visualization in UCSC genome browser

In most cases, bedGraph format provides excellent visualization, especially for alternative splicing events. The BED file after gapless analysis can be converted into a bedGraph file with a script in the countit package.

perl countit/tag2profile.pl -v -big -weight -c ./cache -exact -of bedgraph -n "Sample_bed_graph" sample_out/pair.gapless.bed sample.bedGraph

In some cases, one might also want to show the summarized number of junction reads for each exon junction. This can be done by

perl countit/tag2junction.pl -v -big -weight -c ./cache sample.combined.junction.bed sample.junction.count.bed

The 5th column now contain the read count for each junction, which is also shown as part of the name column. This file can also be loaded to the UCSC genome browser for visualization.

For single end RNA-Seq, the -weight option should not be used.