Gene Structure Annotation and Analysis Using PASA

PASA was originally developed at The Institute for Genomic Research in 2002 as an effort to automatically improve gene structures in Arabidopsis thaliana. Since then, it has been applied to numerous Eukaryotic genome annotation projects including Rice, Aspergillus species, Plasmodium falciparum, Schistosoma mansoni, and Aedes aegypti. We have also successfully applied it to lesser extents in both Mouse and Human among others.

Functions of PASA include:

PASA is composed of a pipeline of utilities that perform the following ordered set of tasks:

PASA is only one component of a larger eukayotic annotation pipeline. Comprehensive genome annotation relies on more than transcript sequence evidence. Not all genes are expressed under assessed conditions, and some genes are expressed at low levels, which complicates their discovery and proper annotation. Other forms of evidence are required for comprehensive genome annotation, including ab initio gene predictors and homology to proteins previously discovered in other sequenced genomes. A complete annotation pipeline, as implemented at the Broad Institute, involves the following steps:

The use of PASA in both applications: first assembling transcript alignments into PASA alignment assemblies, and then later using those PASA assemblies to update EVM consensus (or other) annotations, are described below.

PASA runs on a UNIX/LINUX-based architecture (including mac-osx). PASA involves components written in Perl and C++. Utilities used by PASA, including GMAP, are wrapped by Perl code. Results are stored primarily within a MySQL database, and are available for analysis using the companion suite of Web-based tools and command-line utilities. Running PASA to generate alignment assemblies requires only two inputs: a multi-fasta file for the genome and a multi-fasta file for the transcripts; if FL-cDNAs are available, then an additional file containing listing the accessions of those FL-cDNAs can be provided as a third input. In order to compare the assemblies to existing gene structure annotations, preexisting gene structure annotations can be provided in GFF3 format, or imported by a user-customized data adapter (described below).

Sample data and a preconfigured complete PASA pipeline are available for demonstration purposes, all included in the software distribution.

Download the latest version of the PASA software straight from Sourceforge

or (preferred) pull the bleeding-edge code using CVS like so:

In addition to the PASA software obtained here, you will need the following:

Move the PASA distribution to a location on your filesystem that we can call PASAHOME, such as /usr/local/bin/PASA. From henceforth, we'll refer to this location as $PASAHOME.

The PASA distribution includes the following utilities that you should build and centrally install:

After installing each of the software tools above, all that is needed before running PASA is to configure it. The PASA configuration relies on the file: $PASAHOME/pasa_conf/conf.txt

Recursively copy (cp -r) the $PASAHOME area to the cgi-bin directory of your webserver. Change permissions on everything so that it is world executable (ie. % chmod -R 755 ./PASA ) Now, visit the URL for the status report page for the pasa database you created during the pasa run above.

http://yourServerName/cgi-bin/PASA/cgi-bin/status_report.cgi?db=$mysqldb

This will provide some summary statistics and links to additional web-based utilities for navigating the results from your pasa run.

Now that you have a URL for your base PASA url, update your original configuration file at: $PASAHOME/pasa_conf/conf.txt to set the value of BASE_PASA_URL=http://yourServerName/cgi-bin/PASA/cgi-bin/

For more info, visit the Tour of the PASA web portal

Have each of these files in the same working directory. Then, run the seqclean utility on you transcripts like so:

If you have a database of vector sequences (ie. UniVec), you can screen for vector as part of the cleaning process by running the following instead:

This will generate several output files including transcripts.fasta.cln and transcripts.fasta.clean Both of these can be used as inputs to PASA.

Sample inputs are provided in the $PASAHOME/sample_data directory. We'll use these inputs to demonstrate the breadth of the software application, including using sample DATA ADAPTERs to import existing gene annotations into the database, and tentative structural updates out.

The PASA pipeline requires separate configuration files for the alignment assembly and later annotation comparison steps, and these are configured separately for each run of the PASA pipeline, setting parameters to be used by the various tools and processes executed within the PASA pipeline. Configuration file templates are provided as $PASAHOME/pasa_conf/pasa.alignAssembly.Template.txt and $PASAHOME/pasa_conf/pasa.annotationCompare.Template.txt, and these will be further described when used below.

The next steps explain the current contents of the sample_data directory. You do NOT need to redo these operations:

The following steps, you must execute in order to demonstrate the software. (The impatient can execute the entire pipeline below by running ./run_sample_pipeline.pl. If this is your first time through, it helps to walk through the steps below instead.)

This executes the following operations, generating the corresponding output files: - aligns the all_transcripts.fasta file to genome_sample.fasta using GMAP. Files generated include: sample_mydb_pasa.gmap.tophits :the raw GMAP output sample_mydb_pasa.gmap.tophits.btab :the raw output parsed into a tab-delimited btab output format.

Incorporating PASA Assemblies into Existing Gene Predictions, Changing Exons, Adding UTRs and Alternatively Spliced Models

The PASA software can update any preexisting set of protein-coding gene annotations to incorporate the PASA alignment evidence, correcting exon boundaries, adding UTRs, and models for alternative splicing based on the PASA alignment assemblies generated above.

Comparing to and updating existing gene structure annotations requires that we import these annotations into the PASA database, and are able to extract the suggested updates. PASA utlizes annotation data adapters to achieve this. GFF3 data adapters are included in the PASA distribution, but you can write your own, and directly tie the PASA pipeline to your own informatics infrastructure (ie. other relational database). If you'd prefer to not use GFF3 and to write your own data adapters, visit the PASA data adapter cookbook.

A sample gff3-formatted annotation file is provided in our sample_data directory as orig_annotations_sample.gff3 and can be loaded like so:

Before loading your own GFF3-formatted annotation files, be sure to check them for PASA compatibility like so:

The above gff3-validator will report any entries in your gff3 file that it does not recognize, understand, or otherwise parse properly. It's not a general purpose gff3-validator since it cares only about your protein-coding genes. (note that you should only feed protein-coding genes to PASA using the loader above).

Now that the original annotations are loaded, we can perform a comparison of the PASA alignment assemblies to these preexisting gene annotations, to identify cases where updates can be automatically performed to gene structures in order to incorporate the transcript alignments.

I've copied the ../pasa_conf/pasa.annotationCompare.Template.txt file to our working directory as annotCompare.config. Then, I replaced the MYSQLDB=<MYSQLDB> line with MYSQLDB=sample_mydb_pasa as before with the alignAssembly.config file. Notice this config file contains numerous parameters that can be modified to tune the process to any genome of interest. We'll leave these values untouched for now, relying on the defaults used by PASA, and we'll revisit parameterization later. For most purposes, the defaults are well suited. Run the annotation comparison like so:

Once the annotation comparison is complete, PASA will output a new GFF3 file that contains the PASA-updated version of the genome annotation, including those gene models successfully updated by PASA, and those that remained untouched. This file will be named ${mysql_db}.gene_structures_post_PASA_updates.$pid.gff3, where $pid is the process ID for this annotation comparison computation.

You should revisit the status_report.cgi web page as described above under Setting Up the PASA Web Portal. There, you will be able to navigate the results of the comparison and examine the classifications for annotation updates assigned to each pasa alignment assembly.

Illumina RNA-Seq is quickly revolutionizing gene discovery and gene structure annotation in eukaryotes. Recent enhancements to the PASA pipeline including advancements in RNA-Seq de novo assembly now enable it to make use of these data for gene structure annotation. It is now relatively straightforward to generate strand-specific RNA-Seq data via Illumina. Given the great utility of strand-specific data in differentiating between sense and antisense transcription, plus given the great depth of transcriptome sequencing coverage and the great prevalence of antisense transcription, strand-specific RNA-Seq data is highly preferred by the PASA pipeline. PASA can still be used quite effectively in the case of non-strand-specific RNA-Seq, but the execution is quite different (see below). The dUTP strand-specific RNA-Seq method by Parkhomchuk et al., NAR, 2009 is recommended. For a comparison of strand-specific methods, see Comprehensive comparative analysis of strand-specific RNA sequencing methods. by Levin et al, Nat Methods, 2010.

The procedure for leveraging RNA-Seq in the PASA pipeline is very straightforward. First, assemble the RNA-Seq data using our new Inchworm de novo RNA-Seq assembly software. The RNA-Seq assembly process can be performed in either a genome-guided (recommended) or genome-free way. Documentation for Inchworm assembly (genome-guided or genome-free) is provided at inchworm.sourceforge.net. Instructions for assembly of strand-specific and non-strand-specific RNA-Seq are provided.

In the case of strand-specific RNA-Seq, run PASA with a set of RNA-Seq -specific parameters with the Inchworm transcript assemblies as input:

The above will cluster and assemble alignments with minimal overlap. If your gene density is high and you expect transcripts from neighboring genes to often overlap in their UTR regions, you can perform more stringent clustering of alignments like so:

Also, as an alternative, If you have existing gene structure annotations that are reasonably accurate, you can cluster Inchworm assemblies by locus (annotation-informed clustering) and further augment full-length transcript reconstruction from overlapping inchworm assemblies like so, with the alternative run command:

When PASA is invoked with the —RNA-Seq parameter, it performs filtering of overlapping alignments to exclude likely assembly artifacts based on relative expression values, and it assigns the known transcribed orientation to all alignments.

Beyond these small changes in the PASA alignment assembly invocation, no other changes are required.

In the case of non-strand-specific RNA-Seq, simply treat the Inchworm assemblies like ESTs, and use the basic pipeline command:

The results from running PASA on our sample data set can be examined via the PASA web portal. For example purposes, I've saved a few of the reports generated by the PASA web displays (note pages are generated on-the-fly, however these are provided as static only for example purposes).

The above are just a few examples. Install the PASA portal and navigate your PASA results.

If seqclean was used to clean the transcript sequences, and both the cleaned and original transcript databases were provided in the alignment assembly run of the PASA pipeline as described, then the polyadenylation sites as evidenced in the original transcript sequences and identified as part of the seqclean process were mapped to the genome. The termini of the polyadenylated transcripts are compared to the genome, and those transcripts that truly appear to be polyadenylated and not resulting from an artifact of internal priming to an A-rich region, are reported as candidate polyA sites. The genome coordinate reported as the polyA site is the nucleotide to which polyA is added, so it corresponds to the last non-polyA nucleotide of the polyadenylated transcript. An example of a candidate polyA site can be extracted from one of the log files (default pasa_run.$pid.log/polyAsite_analysis.out) like so:

An additional fasta file (default ${mysql_db}.polyAsites.fasta) summarizes all mapped polyA sites supported by the transcripts. A 100 bp segment of the genome sequence is extracted and oriented, and the last nucleotide in uppercase corresponds to the residue to which polyA is added in the processed transcript. The site corresponding to our example above is as follows:

The accession is bundled like so:

The rest of the header indicates the number of transcripts supporting this polyA site followed by the list of those transcript accessions. The examples above were extracted from our sample data set provided. A more compelling example for Arabidopsis, using spliced transcripts only, is as follows:

PASA is a tool well suited to the identification and classification of alternative splicing isoforms as evidenced by incompatible transcript alignments. Overlapping alignments found incompatible in that they have some structural difference within their overlapping region, and due to their nature of incompatibility, they are relegated to different but overlapping alignment assemblies. PASA performs and all-vs-all comparison among the clustered overlapping alignment assemblies to identify the following categories of splicing variations:

The automated alternative splicing analysis is run as part of the alignment-assembly pipeline.

The results are available in the default output files, with examples below shown from the sample_data/ pipeline run:

File ${mysql_db}.alt_splice_label_combinations.dat :a tab-delimited listing that contains all unique splicing labels for each pasa alignment assembly labeled with a variation. For example:

File ${mysql_db}.indiv_splice_labels_and_coords.dat provides the genome coordinates for each alternative splicing label applied to each corresponding pasa alignment assembly. For example:

The PASA web portal provides numerous reports, graphs, and illustrations to navigate the results of the automated alternative splicing analysis.

In our current working directory, there's a file clusters_of_valid_alignments.txt that contains all the clusters of valid alignments in a simple text format like so:

The transcribed orientation is +,-, or ?. The ? orientation should be used only for single-exon transcript alignments for which the orientation of transcription is ambiguous. By default, PASA assigns all single-exon transcripts that lack evidence of polyadenylation to the ambiguous transcribed orientation. Given this input file, we can demonstrate the pasa alignment assembler like so:

Each cluster of transcript alignments is assembled separately and the results are outputted to stdout with illustrations.

Example input

Corresponding Output

The PASA alignment assemblies can be used to automatically extract gene structures for protein coding genes to be used for training gene finding software. To do this, the longest ORF (min 100 amino acids, configurable) is found within each PASA alignment assembly, allowing for 5' and 3' partials. The top 500 longest ORFs (configurable) are extracted and, using these in addition to the full set of genome sequences, log-likelihood scores for a hexamer-based Markov model are computed. All ORFs are scored using this Markov model, and those ORFs that are most likely coding are reported as best. The scoring method for ORF identification is similar to that used in the GeneID software and nicely described in the methods section of the GeneID publication.

After running the PASA to assemble all transcript alignments as described above, you can run the following.

First, update your PERL library path setting to include the modules that form the core of PASA:

Then, run the following from your PASA working directory. The example below is what I would run in the sample_data/ directory:

This should generate a series of files, described below:

The best_candidates.pep file has headers like so:

The accession corresponds to the PASA assembly. The contig and the coordinate span of the alignment assembly is provided. The type indicator can be any of the following: complete, 5prime_partial, 3prime_partial, or internal. The 5prime_partial are missing a start codon and translate to the very 5' end, 3prime_partial are missing a stop codon and translate to their very 3' end, and internal translate from the first to the last basepair in the sequence, missing a start and a stop codon. The complete category are of greatest interest for the prospects of genefinder training. Typically, we would search these complete proteins against the nonredundant protein database at GenBank and identify those ORFs that have good database matches across most of their length. Such entries can be confidently used for training, in addition to those particularly long ORFs that do not match known proteins and are sufficiently complex in sequence.

The PASA software and its original application are described in:

The use of PASA to analyze polyadenylation signals is described in:

Enhancements to PASA that automate the identifiation and classification of alternative splicing variations are described here:

Incorporation of RNA-Seq data into gene structure annotation improvements using PASA and Inchworm: