Whole genome alignment howto

The whole genome alignments are definitely the biggest mystery of the UCSC browser for me. The scripts like doBlastzChainNet.pl 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. In addition, it's using the new hgAutomate libraries. That said, doBlastzChainNet.pl was very helpful when writing this just as was Conservation Track and this page which is very similar.

Thanks a lot to Angie and Hiram for answering all my questions!

Before you start

To understand the following you need to:

1. Know about shell scripting. Have a look on the Bash Scripting Guide
2. Have your own browser installed and running, see Learn_about_the_Browser and AngieTest
3. Have created some tracks yourself, loaded them with hgLoadBed and configured them in trackDb.ra, see Learn_about_the_Browser and Adding New Tracks to a browser installation

Fileformats we have to know

I already knew the basic file formats like bed/gff/etc, but hadn't crossed some of these yet:

1. lav: a compact form to store genomic pairwise alignments created by Blastz, using only numbers (pos of match + identities)
2. psl: similar to lav, but the most compact format at UCSC to store alignments, easier to parse than lav (positions of match + identities), e.g. blat is creating this
3. axt: a much more readable way to store *pairwise* alignments: positions + aligned sequences (It seems that you can forget about this format, it's rarely used, but the tools still contain axt... in their name although most also accept psl. )
4. maf: an extended version of axt, *multiple* genomic alignments: assemblies + positions + aligned sequences

You see: To convert from lav to axt or 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: Two matching fragments next to each other are joined into one fragment (axtChain). 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 (more info on the mailing list)
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
6. PhastCons: Calculate the strength of conservation for every base, similar to a Vista- or protein Conservation plot, but applicable to multiple alignments

Outline - More Details

• 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 took ages on my own computer. I'm Europen, so I got an account vital-it, anyone can apply who is working in the EU. There are national compute cluster also in the Netherlands and the UK that I know of.
• 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. (Sometimes UCSC also uses netfilter to filter out nets with a gapsize <= 12)
• We feed the tree into multiz. Multiz will use many local alignments to generate a multiple local alignment.

 Example: A-B and A-C have been aligned with blastz...
    A aaactg  \
    B aa--tg   \     A aaa-ctg
                ->   B aa---tg
    A aaa-ctg  /     C aattt-g
    C aattt-g /

• Multiz is not a global aligner, as it uses Blatz
• 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. (according to Angie) and is only used for the encode regions currently.
• The only other aligner on whole genomes that I know of is SLAGAN. In theory, it's output can be converted to lav but I doubt that anyone has ever tried this. SLAGAN is not related to Blastz. I have alignments for Ci/Cs made by both algorithms and they differ quite a lot, however, I've never explored the subject in detail in the end and will probably stick with Blastz

Example, step 1: Alignments with Blastz

• I want to do a whole genome alignment of C. intestinalis V2 and C.savignyi V2. This is rather small genome (150 MB).
  • Downloaded ci2 from hgdownload.cse.ucsc.edu into the directory ci2
  • Downloaded cs2 (only visible on hgwtest.cse.ucsc.edu, not on the real server nor the download server) from Ensembl into the directory cs2
• Sometimes genomes come in one big fasta file, we want to split it into smaller ones:

  faSplit byName cs2/cs2.fa cs2/
  rm -f cs2/cs2.fa

• The Ensembl genome cs2 does not seem to be tandem repeat masked yet, so I do:

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

• Blastz seems to ignore lower case characters (took me one afternoon to figure this out), we simply convert everything to nib-format:

  for i in in ci2/* cs2/*; do faToNib $i `echo $i | sed -e s/.fa/.nib/`; done

• Then, to blastz everything onto everything, in theory we could issue this command:

  for i in cs2/*.nib; do echo 'for j in ci2/*.nib; do blastz '$i' $j  H=2000 Y=3400 L=6000 K=2200 Q=HoxD55.q > lav/`basename '$i' .nib`-`base name $j .nib`.lav; done'; done

• However, I better use the cluster at vital-it.ch: I rsync the nib files onto it and add a cluster instruction:

  for i in cs2/*.nib; do echo 'for j in ci2/*.nib; do bsub blastz '$i' $j  H=2000 Y=3400 L=6000 K=2200 Q=HoxD55.q > lav/`basename '$i' .nib`-`base name $j .nib`.lav; done'; done

• • In my case I have 4000 scaffolds on one genome and 400 reftigs on the other. That was still doable with the cluster.
  • However, if both genomes are spread across many more scaffolds you have to play around with the sequences, otherwise it will be millions of blastz-runs. At UCSC they first join the smaller scaffolds into one big chromsome and then run blastz on the "ChrUn"-virtual-chromosome. (this will create some false nets, going from one scaffold to another, but we don't have a choice here (hiram))

Example, step 2: Chaining

• Convert the lav files into the more compact psl format:

  for i in *.lav; do echo $i; lavToPsl $i `basename $i .lav`.psl; done; 
  rm -f *.lav
  less chr01p.psl

• Oups. Now, once again, I've messed up target and query. My psl-files actually look like "... chrom01 ... reftig1..." (the format is query, target). However, since I want to annotate Ci2, which is assembled into chromosomes, chrom01 should really be the target. I have to swap these fields in the current psl files and re-create new psl-files, split by the name of the target:

 cat *.psl > ../all.psl
 pslSwap ../all.psl ../all-swap.psl
 pslSplitOnTarget all-swap.psl psl/ -lump

• Now the chaining (this is quite fast for my little ciona-genome):

 for i in psl/*.psl; do echo $i; axtChain $i ci2 cs2 chain/`basename $i .psl`.chain -linearGap=loose -psl; done

• For the filtering of the chains, we need the size of each chromosome:

 faSize ../genome/ci2.fa -detailed > ci2.sizes
 faSize ../../cs/cs2/cioSav2.fa -detailed > cs2.sizes

• Sort and filter the chains (more info in this post of the mailing list):

  chainMergeSort chain/*.chain > all.chain
  chainPreNet all.chain ci2.sizes cs2.sizes all.pre.chain

• We load these chains into the browser (hm...wouldn't it be more interesting to have the complete chains in the browser?):

 hgLoadChain ci2 chainCioSav2 all.pre.chain

• Adding a section similar to the following to trackDb.ra, running make alpha DBS=<yourdbname> ZOO= there, the track should then appear on the browser:

 track chainCioSav2
 shortLabel C. savignyi chain
 longLabel C. savignyi chain
 group compGeno
 priority 125
 visibility hide
 color 100,50,0
 altColor 255,240,200
 spectrum on
 type chain cioSav2
 otherDb cioSav2

Example, step 3: Netting

• Then the netting itself, combine the chains into nets and add synteny information (not sure where this is displayed, probably needed for double/single-lines display of net-tracks) to them:

  chainNet all.pre.chain -minSpace=1 ci2.sizes cs2.sizes stdout /dev/null | netSyntenic stdin noClass.net

• We need to add additional information to this nets (shown on the details page):

  netClass -noAr noClass.net ci2 cioSav2 cioSav2.net

NB: I got an error message here because of the gap.bin field in cioSav2 which does not seem to exist in lib/agpGap.c, function agpGapLoad, but is retrieved by "select * from gap order by chrom". I don't get it, I just updated my source tree. I hacked around the problem but there shouldn't really be any error message.

• Finally, we can load these nets into our own browser:

  hgLoadNet ci2 netCioSav2 cioSav2.net

adding something like this to trackDb.ra:

 track netCioSav2
 shortLabel C. savignyi Net
 longLabel $o_Organism ($o_date/$o_db) Alignment Net
 group compGeno
 priority 134
 visibility dense
 spectrum on
 type netAlign cioSav2 chainCioSav2
 otherDb cioSav2

Example, step 4: Maffing

• I needs mafs for the browser as I want to display the alignments:

 netToAxt cioSav2.net all.pre.chain ci2/ cs2/ ci2-cioSav2.axt
 axtToMaf ci2-cioSav2.axt ci2.sizes cs2.sizes ci2-cioSav2.maf -tPrefix=ci2. -qPrefix=cioSav2.

• Load the data into the brower. It needs an mysql-index to find the right position. The alignment files go into /gbdb:

 hgLoadMafSummary ci2 mafCioSav2Summary ci2-cioSav2.maf
 sudo mkdir /gbdb/ci2/mafCioSav2
 sudo cp ci2-cioSav2.maf /gbdb/ci2/mafCioSav2/mafCioSav2.maf
 hgLoadMaf ci2 mafCioSav2

Example, step 5: Phastcons

• Phastcons doesn't seem to like big maf files that include many chromsomes, so I split the maf files by chromosome:

  mkdir maf
  mafSplit -byTarget dummy.bed maf/ ci2-cioSav2.maf

• I take the biggest chromsome and calculate the background model from it:

  PHAST=~/bio/phastCons
  ${PHAST}/bin/phyloFit -i MAF maf/002.maf

• Run phastCons on every maf file with this background model and some obscure parameters that phastCons needs. I'm parsing the chrosomome name here from the maf as it was lost during the mafSplitByTarget (there should be an easier way but it works):

 mkdir wig
 mkdir mostCons
 for i in maf/*.maf; do \
         ~/bio/phastCons/bin/phastCons --target-coverage 0.25 --expected-length 12 \
          --rho 0.4 --msa-format MAF $i phyloFit.mod \
          --seqname `cat $i | head -n 3  | tail -n 1 | tr -s ' ' | cut -f 2 -d ' ' | cut -d. -f2`
          --most-conserved mostCons/`basename $i .maf`.bed > wig/`basename $i .maf`.wig; \
 done

This took 7 minutes for a maf file of 113 MB

Hiram: How can I check if these parameters are "good"? Use CDS overlap?

• Load the wiggle files into the browser:

 cat wig/*.wig > ../ci2-cioSav2.pp
 wigEncode ci2-cioSav2.pp ci2-cioSav2.wig ci2-cioSav2.wib
 sudo mkdir /gbdb/ci2/wib/
 sudo cp ci2-cioSav2.wib /gbdb/ci2/wib/
 hgLoadWiggle  -pathPrefix=/gbdb/ci2/wib ci2 phastCons ci2-cioSav2.wig

• Load the "most Conserved" regions into the browser:

 cat mostCons/*.bed > mostCons.bed
 hgLoadBed ci2 mostConserved mostCons.bed

• ... add something like this entry in trackDb.ra:

 track mafCioSav2
 shortLabel Conservation
 longLabel Blastz/Phastcons conservation C. savignyi
 group compGeno
 visibility full
 type wigmaf 0.0 1.0
 speciesOrder ci2 cioSav2
 wiggle phastCons
 viewLimits 0:1
 autoScale off
 maxHeightPixels 50:50:11

• Refreshing your trackDb (make ...) should then display the track

Remarks

• The original authors already wrote similar tools: 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.
• Later on, 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."
• Why do mafs include the chrosome sizes?
• GERP is an alternative for phastCons
• I still wonder if filtering for synteny isn't a bit too stringent sometimes.. The "pure" chains, however seem to include a lot of noise. Isn't there something in between? Chains filtering with their score might be an option.