Revision as of 10:06, 2 November 2013

This is a protected page.

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.

Versions

• v1.0.2 (05-31-2013), current
  • 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, although the pipeline can be expanded to examine novel transcripts as well.

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 available at http://zhanglab.c2b2.columbia.edu/index.php/OLego. You can download the source codes or precompiled binaries and the documentation.

The program takes i) read file(s) in either FASTA or FASTQ format, together with ii) the indexed genome and iii) (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. 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 include 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.

• Install the perl library files czplib

czplib is a perl library with various functions for genomic/bioinformatic analysis. (download from SourceForge.net)

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

https://sourceforge.net/projects/ngs-quantas/

$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
• human (hg19): http://zhanglab.c2b2.columbia.edu/data/Quantas/data/hg19.tgz

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.

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

olego -v -t 16 -w 14 -M 4 --max-multi 10 -r olego_dir/models/mm9.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 are used, the memory footprint is in general < 4Gb. That said, the exact amount of memory also depends on the seed size and the repetitive nature of reads. In terms of speed, our estimate is that it can align 200 M 100x2 PE reads in ~29 hrs with 8 cores (2.0G Hz).

This will map reads to the exons, known junctions as well as novel junctions. By default OLego searches with perfect matches in seeds of 14 nt, allowing for a maximum of 4 substations or indels along the whole read (which should be good 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

Important note: this should be 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 package.

perl sam2bed.pl --uniq -v sample.r1.sam sample.r1.uniq.bed
perl sam2bed.pl --uniq -v sample.r2.sam sample.r2.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-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 ~/czsrc/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.

-uniq means only uniquely mapped reads (single hits) were 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).

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 was 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

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 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, and the -weight option is specified because each fragment from gapless now has a probability score.

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 meaningless, 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 not to run all types of AS events.

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 ~/czsrc/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 change.

Two notes:

1. A caveat of the Fisher exact test described above is that it does not take into account the within-group variation. This is in general a not big issue for lab mice if RNA-seq libraries are prepared properly because the heterogeneity between lab mice is minimal. So far I have not seen tools that deal with this issue effectively, but this could be due to the fact that I did not track RNA-seq literature very closely.
2. Also, a recent paper described a method MATS based on MCMC to estimate p-values and claims that the performance is somewhat better than Fisher's exact test in simulation. I have not tested this, but if this is true, the test can be swathed very easily.

Quantify gene expression

This is also done with the countit package, but first one needs to download gene definition files from http://geller.rockefeller.edu/Quantas/data/mm9.gene.tar.gz. There are two files in the tarball that will be used in the next step.

perl summarize_expression_wrapper.pl -big --cache ../cache -exon mm9.exon.uniq.core.bed -e2g mm9.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.

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 was done for all samples. Then an R package edgeR was used to do statistical test of differential gene expression. An example script can be downloaded from: http://geller.rockefeller.edu/Quantas/expr.edgeR.tagwise.R, which will report log2FC, p-value, and benjamini FDR

RPKM for each group needs to be recalculated after combining all samples in the same group.

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 tag2profile.pl -big -weight -exact -of bedgraph -n "Sample_bed_graph" -v 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 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.