TOAST and ROAST

TABLE OF CONTENTS

• Introduction
• Installation
• Example
• Data Flow
• File and Tree Formats
• Program Descriptions

Introduction

This package contains the tools TOAST and ROAST for finding orthologous alignments of DNA sequences, plus some associated utility programs.

TOAST -- two way orthologous selection tool

ROAST -- reference dependent multiple alignment tool

Installation

Download the distribution file into an empty source directory, and unpack it using a command like:

    zcat toast-roast.tar.gz | tar xvf -

    make
    make install

Example

Here is a small example that aligns the beta-globin gene cluster from four species, with human as the reference, using sample input files included with the distribution. The input consists of sequence files in FASTA format for human, marmoset, mouse, and monodelphis, plus an annotation file chr11.genes.txt which specifies the gene locations in the human sequence. The annotation file is optional, but is useful when the genes have long introns, and may help TOAST produce better output.

To run this example for yourself, type:

    all_bz + A=2 F=human D=0 T=chr11.genes.txt "(((human marmoset) mouse) monodelphis)"
    roast + X=2 E=human "(((human marmoset) mouse) monodelphis)" *.toast2.maf example.roast.maf

Data Flow

The all_bz program automates the creation of preliminary pairwise alignments for the necessary pairs of species, using blastz, toast, and several utility tools such as single_cov2 and chain. The roast program automates the progressive creation of multi-species alignments from the pairwise ones using the multiz or multic programs, which are tools for aligning alignments.

File and Tree Formats

Input Sequences

Each species has at most one sequence file, which includes one or more sequences from the same or different chromosomes or contigs. The sequence file name is usually the species or assembly name, and it must be consistent with string1 (below).

Sequences must be supplied in FASTA format. The primary form for the headers is

    >string1:string2:int1:char:int2

• string1 is the species or assembly name;
• string2 is the chromosome or contig;
• int1 is the sequence start position (1-based, inclusive) on the chromosome/contig or its reverse complement;
• char is + or -, denoting the strand;
  (If this is -, the input sequence is in reverse complement with respect to the chromosome/contig, and int1 is specified from the end of the entire chromosome/contig back to the nearest endpoint of the supplied sequence.)
• int2 is the size of the entire chromosome/contig.

Alternatively, MSA header format is also supported. (This has been wrapped here for easier readability, but in actual use must appear all on one line.) In this format an empty field is represented by a dot ".".

    >${COMMON_NAME}|${ENCODE_REGION}|${FREEZE_DATE}|${NCBI_TAXON_ID}
    |${ASSEMBLY_PROVIDER}|${ASSEMBLY_DATE}|${ASSEMBLY_ID}|${CHROMOSOME}
    |${CHROMOSOME_START}|${CHROMOSOME_END}|${CHROM_LENGTH}|${STRAND}
    |${ACCESSION}.${VERSION}|${NUM_BASES}|${NUM_N}|${THIS_CONTIG_NUM}
    |${TOTAL_NUM_CONTIGS}|${OTHER_COMMENTS}

If there is only one sequence in the file, a header like

    >string1

In the sequence itself, TOAST and ROAST support the same characters as blastz, including lowercase letters and N to represent unsequenced positions. Other ambiguity codes such as W, R, -, etc. are not allowed.

Annotation File

Each line in this file records the location of a gene in the reference species (which must have only one sequence for this to be supported). If the FASTA header for the reference sequence specifies its location within a larger chromosome or contig, these gene positions must be relative to the entire chromosome/contig; otherwise they must be relative to the given input sequence. The intervals can be either 0-based or 1-based, and half-open or inclusive, since 1 bp either way doesn't matter much for this purpose. Note that start <= end regardless of gene orientation, and the gene intervals cannot overlap.

    start end

Species Guide Tree

A species-guide-tree is used as an argument in many of the programs to indicate the evolutionary relationships among the input species. It is a binary tree built from the species/assembly names using parentheses, and is enclosed in double quotes, e.g. "((human chimp)(rat mouse))". Note that the species/assembly names must be consistent with their sequences' file names, and also with the string1 field in the sequence header.

MAF Alignment Output

Most of these programs produce alignment output in UCSC's MAF format. Even though the raw output from blastz is in another format called LAV, all_bz converts this to MAF automatically, using the included utility tool lav2maf. If the FASTA header for an input sequence specifies its location within a larger chromosome or contig, the output coordinates for that sequence will be relative to the entire chromosome/contig; otherwise they will be relative to the given input sequence (this is true for the intermediate pairwise output as well). Except for the main programs all_bz (which creates many MAF files with names of its own choosing) and roast (which uses a final output filename that you specify), most of the other, helper programs write their output to stdout so it can be piped, redirected, etc.

Program Descriptions

For all programs below, square brackets [] denote optional parameters, while question marks ? and items in angle brackets <> represent meta-syntactic variables that should be replaced with your numbers or names. Don't type any of the brackets or the question marks. Some of the command lines have been wrapped here for easier readability, but they should always be typed all on one line; the order of the arguments is not important. In the parameter descriptions, default values are given in parentheses (). Not all of the programs are listed here, but in general, running any of them with no arguments will print a help message with command-line syntax, similar to the following.

all_bz -- runs blastz and post-processing commands for all necessary pairs of the specified sequences.

    all_bz [+-] [b=?] [A=?] [F=<reference-species>] [T=<annotation-file>]
      [h=?] [q=?] [D=?] [c=?] [f=?] <species-guide-tree> [<blastz-specfile>]

toast -- Two-way Orthologous Assignment Selection Tool; produces orthologous alignments from blastz alignments.

    toast [q=?] [h=?] [s=?] [f=?] [A=<annotation-file>] <seq1-seq2-maf>
      <seq1-seq1-maf> <seq2-seq2-maf>

roast -- reference dependent multiple aligner; uses pairwise alignments produced by all_bz.

    roast [+-] [R=?] [M=?] [P=?] [T=<temp-dir>] [X=?] [C=?] E=<reference-species>
      <species-guide-tree> <maf-source> <destination>

blastz_clean -- screens out close overlapped alignments: whenever two alignments overlap in one sequence and are nearby in the other sequence, the overlapping part of the shorter one is removed.

    blastz_clean [c=?] <pairwise-maf-file>

chain -- groups alignment blocks into chains (see publication) and projects the chains onto a reference species.

    chain [R=<reference-species>] [A=<annotation-file>] <pairwise-maf-file>

maf_order -- orders the rows in MAF blocks according to a given list of species/assemblies (regardless of chromosome/contig); rows not specified are removed.

    maf_order <maf-file> <species1> <species2> ... [nohead] [all]

maf_project -- extracts MAF blocks that include a given reference sequence or species.

    maf_project <maf-file> <reference> [<from> <to>] [<filename-for-removed>]
      [<species-guide-tree>] [nohead]

maf2lav -- converts a MAF file to the format used by PipMaker and the Laj alignment viewer. (This is not needed for the Gmaj alignment viewer, which reads MAF files directly.)

    maf2lav <pairwise-maf-file> <seq1-file> <seq2-file>

single_cov2 -- screens out all overlapped alignments: whenever two alignments overlap in the specified species, the overlapping part of the shorter one is removed.

    single_cov2 <pairwise-maf-file> [R=<species>] [F=<filename-for-removed>]

transformAnnotation -- transforms a gene annotation file in a format like the knownGene or refGene tables at the UCSC Table Browser to the simpler format used by argument A of toast. It removes overlaps in a heuristic fashion while attempting to have each resulting interval still correspond to a gene. (Note that this program extracts the transcription start and end fields without doing any coordinate conversion or chromosome filtering; if you are starting with the full UCSC table files, you will need to grep for your chromosome of interest before running this.)

    transformAnnotation <known-gene-file>