Whole genome alignment howto

The whole genome alignments are definitely the biggest mystery of the UCSC browser for me. The scripts like doBlastz and collegues automize the whole stuff for the folks working at UCSC but doesn't really make it any easier to fully understand the system as everything is buried now even more and full of parasol statements.

So I'm trying to re-create whole genome alignments to better understand how my whole-genome alignments are created. Thanks a lot to Angie for answering my questions!

Fileformats we have to know

1. lav: a compact form to store genomic pairwise alignments, using only numbers (pos of match + identities)
2. axt: a more human way to store *pairwise* alignments: positions + aligned sequences
3. maf: an extended version of axt, *multiple* genomic alignments: assemblies + positions + aligned sequences

So to convert from lav to axt/maf we need the genomic sequences from the nib/fa files.

Outline

1. MASKING: Both genomes have to be repeatmasked and masked Tandem Repeat Finder (trf) first (thanks to Hiram for pointing this out)
2. ALIGNING: The two genomes are aligned with BLASTZ (we don't use blastz's own chaining, see discussion (angie)). This generates lav-files, which have to be converted to psl (lavToPsl)
3. CHAINING: As every genomic fragment can match with several others, we keep only the best match for a given part : first do axtSort, then filter with axtBest
4. NETTING: Group blocks of chained alignments into longer stretches of synteny (netChain)
5. MAF'ING: From the synteny-files (positions), get the sequences and re-create alignments

Outline - Normal case: Many genomes, one is the reference

• BLASTZ: The reference genome is aligned with all others with BLASTZ. That creates lav-files. They are converted to psl.
  • Angie says: Our hgdownload $db/vs$OtherDb/README.txt pages give the blastz parameters used for a given run.
  • It's a good idea to do this on a compute cluster. It will take ages if you try to do this for a larger genome on one single computer. If you're Europen, you can get an account on vital-it. There are national compute cluster also in the Netherlands and the UK that I know of.
• Two matching fragments next to each other are joined into one fragment (axtChain, "chaining")
• Chains are simply better alignments: If several alignments overlap, we still don't know which the best one is. This filtering was done with axtBest before, now we use chainNet + netToAxt . (I have no clue why Kate mentioned netFilter here.)
• We feed the tree into multiz. Multiz will use many local alignments to generate a multiple local alignment.

 Example:
    A aaactg  \
    B aa--tg   \     A aaa-ctg
                ->   B aa---tg
    A aaa-ctg  /     C aattt-g
    C aattt-g /

• Multiz is not an aligner at all. It's just a "reformatter", rewriting pairwise alignments into multiple alignments. It needs a phylogenetic tree.
• TBA is more like a real aligner. But it's too slow currently for massiv genome-wide alignments.

Example - Simple case

• I want to do a whole genome alignment of C. intestinalis V2 and C.savignyi V2.
  • Downloaded ci2 from hgdownload.cse.ucsc.edu
  • Downloaded cs2 (only visible on hgwtest.cse.ucsc.edu, not on the real server nor the download server) from Ensembl.
• The Ensembl genome cioSav2 does not seem to be tandem repeat masked yet, so I do:

for i in cioSav2/*.fa; do trfBig $i trf/`basename $i`; done

Remarks

• For the history records: Axt started in mouseStuff and maf started with the ratStuff. There is an older tool called *axtBest* but Angies says: "axtBest is ancient history. It has been replaced by the chaining and netting process, which does a better job of finding the "best" alignment to cover a given region."
• Multiz contains a tool that converts lav2maf directly and UCSC includes one with lavToMaf. However, we don't care about fragments that match two times. For a whole genome, you really want "chains" of best-matching fragments, therefore we don't use these tools.

• Why the heck do mafs include the chrosome sizes?