Quantas Documentation
From Zhang Laboratory
This is a protected page.
Contents
- 1 Introduction
- 2 Versions
- 3 Important assumptions
- 4 Software installation
- 5 Mapping RNA-seq reads to the reference genome and exon junctions
- 6 Inferring transcript structure between paired end reads
- 7 Quantify alternative splicing
- 8 Quantify gene expression
- 9 Generate bedGraph files for visualization in UCSC genome browser
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 in its relatively early stage. Any comments are welcome and can be directed to Chaolin Zhang.
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. 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. 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.
- 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
http://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 (24Mb)
- human (hg19): http://zhanglab.c2b2.columbia.edu/data/Quantas/data/hg19.tgz (31M)
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. 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 15 nt with one 1-ht 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
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-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 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).
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 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. case).
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:
- 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.
- 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.