AnchorWave (Anchored Wavefront Alignment) identifies collinear regions via conserved anchors (full-length CDS and full-length exon have been implemented currently) and breaks collinear regions into shorter fragments, i.e., anchor and inter-anchor intervals. By performing sensitive sequence alignment for each shorter interval via a 2-piece affine gap cost strategy and merging them together, AnchorWave generates a whole-genome alignment for each collinear block. AnchorWave implements commands to guide collinear block identification with or without chromosomal rearrangements and provides options to use known polyploidy levels or whole-genome duplications to inform alignment.
AnchorWave takes the reference genome sequence and gene annotation as input and extracts reference full-length coding sequences (CDS) to use as anchors. We use minimap2 to lift over the position of reference full-length CDS to the query genome (step 1). AnchorWave then identifies collinear anchors using one of three user-specified algorithm options (step 2) and uses the [WFA](https://github.com/smarco/WFA) algorithm to perform alignment for each anchor and inter anchor interval (step 4). Some anchor/inter-anchor regions cannot be aligned using our standard approach due to high memory and computational time costs. For these, AnchorWave either identifies novel anchors within long inter-anchor regions (step 3), or for those that cannot be split by novel anchors, aligns using the ksw_extd2_sse function implemented in minimap2 or a reimplemented sliding window approach (step 4). AnchorWave concatenates base pair sequence alignment for each anchor and inter-anchor region and outputs the alignment in MAF format (step 5).GNU GCC >=7.0
Cmake >= 3.0
minimap2
git clone https://github.com/baoxingsong/anchorwave.git
cd anchorwave
cmake ./
make
You will get an executable file named anchorwave
. If you would like to move anchorwave
to a new place, please always put configure
and anchorwave
into the same folder.
The code has been tested under Ubuntu 20.2 and CentOS 7 with intel/AMD CPU. It should work well on other REDHAT or Debian based Linux Distributions.
In general, totally four commands are need to run through the whole pipeline.
- extract CDS
- align CDS to the reference genome
- align CDS to the query genome
- perform genome alignment
The command 4 need to be well tuned when comparing genomes with chromosomal rearrangement or whole genome duplication variation
Options:
Program anchorwave
Usage: anchorwave <command> [options]
Commands:
gff2seq get the longest full-length CDS for each gene
genoAli whole chromosome global alignment and variant calling
proali genome alignment with relocation variation, chromosome fusion or whole genome duplication
When extracting full-length CDS, if for a gene, there are multiple transcripts, only the transcript with longest full-length CDS would be used.
The gff2seq
output the concatenated full-length CDSs.
Options of the gff2seq
function:
Usage: anchorwave gff2seq -i inputGffFile -r inputGenome -o outputSequences
Options
-h produce help message
-i FILE reference genome annotation in GFF/GTF format
-r FILE reference genome sequence in fasta format
-o FILE output file of the longest CDS/exon for each gene
-x use exon records instead of CDS from the GFF file
-m INT minimum exon length to output (default: 20)
Since minimap2 could not deal with short CDS very well, and that causes error to lift over anchors to the query genome. To minimum this side effects, AnchorWave would ignore those short (-m parameter) CDS records.
Data: Arabidopsis thaliana Col-o reference genome and GFF3 annotation file from https://www.arabidopsis.org/
Arabidopsis thaliana Ler-0 accession assembly from http://www.pnas.org/content/113/28/E4052
anchorwave gff2seq -i TAIR10_GFF3_genes_no_UM.gff -r tair10.fa -o cds.fa
minimap2 -x splice -t 10 -k 12 -a -p 0.4 -N 20 tair10.fa cds.fa > ref.sam
minimap2 -x splice -t 10 -k 12 -a -p 0.4 -N 20 ler.fa cds.fa > ler.sam
This module perform base pair resolution sequence alignment for two genomes. A query chromosome sequence would be aligned against the reference chromosome with the same name.
The output would be an end-to-end sequence alignment for the whole chromosome in maf format.
Or alignment for each anchor region and inter anchor region.
A variant calling result in vcf format could be created which is derived from end-to-end alignment.
Please make sure the chromosomes from reference genome and query genomes were named in the same way. Chromosomes with the same name would be aligned.
grep ">" ler.fa
grep ">" Col.fa
Anchors lift over:
anchorwave gff2seq -i TAIR10_GFF3_genes_no_UM.gff -r tair10.fa -o cds.fa
minimap2 -x splice -t 10 -k 12 -a -p 0.4 -N 20 tair10.fa cds.fa > ref.sam
minimap2 -x splice -t 10 -k 12 -a -p 0.4 -N 20 ler.fa cds.fa > ler.sam
Perform alignment and variant calling:
anchorwave genoAli -i TAIR10_GFF3_genes.gff -as cds.fa -r tair10.fa -a ler.sam -ar ref.sam -s ler.fa -v ler.vcf -n ler.anchors -o ler.maf -f ler.f.maf > ler.log
The output file ler.anchors
could be used to visualize collinear anchors/blocks.
Please note: Under the global alignment model, if the block length is longer than a preset window size, the output alignment maybe a suboptimal alignment. While a large window size (larger values of parameters -w
and -fa3
) would cost more memory.
Options of the genoAli
function:
Usage: anchorwave genoAli -i refGffFile -r refGenome -a cds.sam -as cds.fa -ar ref.sam -s targetGenome -n outputAnchorFile -o output.maf -f output.fragmentation.maf
Options
-h produce help message
-i FILE reference GFF/GTF file
-r FILE reference genome sequence file in fasta format
-as FILE anchor sequence file. (output from the gff2seq command)
-a FILE sam file generated by mapping conserved sequence to query genome
-s FILE target genome sequence file in fasta format
Those sequences with the same name in the reference genome and query genome file would be aligned
-n FILE output anchors file
-o FILE output file in maf format
-f FILE output sequence alignment for each anchor/inter-anchor region in maf format
-t INT number of threads (default: 1)
-v FILE output variant calling in vcf format (conflict with -IV)
-m INT minimum exon length to use (default: 20, should be identical with the setting of gff2seq function)
-mi DOUBLE minimum full-length CDS anchor hit similarity to use (default:0.95)
-mi2 DOUBLE minimum novel anchor hit similarity to use (default:0.2)
-ar FILE sam file generated by mapping conserved sequence to reference genome
-w INT sequence alignment window width (default: 30000)
-fa3 INT if the inter-anchor length is shorter than this value, stop trying to find new anchors (default: 100000)
-B INT mismatching penalty (default: -2)
-O1 INT gap open penalty (default: -3)
-E1 INT gap extension penalty (default: -2)
-O2 INT gap open penalty 2 (default: -80)
-E2 INT gap extension penalty 2 (default: -1)
-IV whether to call inversions (default: false)
-IC DOUBLE penalty for having a non-linear match in inversion region (default: -1)
We use IC * alignment_similarity as the penalty in the inversion block
-I DOUBLE minimum score to keep an inversion (default: 2)
-e INT maximum expected copy number of each gene on each chromosome (default: 1)
This prevents using tandem duplicated genes to identify collinear block
-y DOUBLE minimal ratio of e+1 similarity to 1 similarity to drop an anchor (default: 0.6)
-ns do not search for new anchors (default: false)
-x use exon records instead of CDS from the GFF file
This module perform base pair resolution sequence alignment for two genomes. A query chromosome sequence would be aligned against the reference chromosome with the same name end-to-end.
The output could be end-to-end sequence alignment for the whole chromosome in maf format. Or alignment for each anchor region and inter anchor region.
data: Maize B73 v4 reference genome and GFF3 annotation file;
Maize Mo17 assembly.
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-34/gff3/zea_mays/Zea_mays.AGPv4.34.gff3.gz
gunzip Zea_mays.AGPv4.34.gff3.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-34/fasta/zea_mays/dna/Zea_mays.AGPv4.dna.toplevel.fa.gz
gunzip Zea_mays.AGPv4.dna.toplevel.fa.gz
wget https://download.maizegdb.org/Zm-Mo17-REFERENCE-CAU-1.0/Zm-Mo17-REFERENCE-CAU-1.0.fa.gz
gunzip Zm-Mo17-REFERENCE-CAU-1.0.fa.gz
Please make sure the chromosomes from reference genome and query genomes were named in the same way. Chromosomes with the same names would be aligned.
grep ">" Zm-Mo17-REFERENCE-CAU-1.0.fa
grep ">" Zea_mays.AGPv4.dna.toplevel.fa
Use this command to rename each chromosome for Zm-Mo17-REFERENCE-CAU-1.0.fa:
sed -i 's/>chr/>/g' Zm-Mo17-REFERENCE-CAU-1.0.fa
Anchors lift over:
anchorwave gff2seq -i Zea_mays.AGPv4.34.gff3 -r Zea_mays.AGPv4.dna.toplevel.fa -o cds.fa
minimap2 -x splice -t 10 -k 12 -a -p 0.4 -N 20 Zm-Mo17-REFERENCE-CAU-1.0.fa cds.fa > cds.sam
minimap2 -x splice -t 10 -k 12 -a -p 0.4 -N 20 Zea_mays.AGPv4.dna.toplevel.fa cds.fa > ref.sam
Perform alignment using the genoAli
function with parameter -IV
:
anchorwave genoAli -i Zea_mays.AGPv4.34.gff3 -as cds.fa -r Zea_mays.AGPv4.dna.toplevel.fa -a cds.sam -ar ref.sam -s Zm-Mo17-REFERENCE-CAU-1.0.fa -n anchors -o anchorwave.maf -f anchorwave.f.maf -w 38000 -fa3 200000 -B -6 -O1 -8 -E1 -2 -O2 -75 -E2 -1 -t 1 -IV
This command might cost a couple of days and 50 Gb memory. Reduce the values of parameters -w
and -fa3
could reduce memory usage but would also reduce the sequence alignment quality.
If you have larger memory and multiple CPU cores available, you could increase the value of -t
to run the command using multiple threads.
Genome alignment with relocation variation, chromosome fusion or whole genome duplication (an option of command 4)
When comparing two genomes undergone different genome duplications (goldfish genome, plant genomes and maybe hexapod genomes), this program implemented algorithm to identify collinear blocks with user specified coverage, and perform base pair resolution genome alignment for each collinear block.
The outputs are end-to-end sequence alignments for each collinear block in maf format.
Data: maize B73 genome sequence and GFF3 annotation file from http://plants.ensembl.org/Zea_mays/Info/Index;
and sorghum genome sequence from https://plants.ensembl.org/Sorghum_bicolor/Info/Index.
This function does not need the reference and query genomes having identical chromosome names.
Extract CDS of reference genome and map them to the query genome using minimap2:
anchorwave gff2seq -r Zea_mays.AGPv4.dna.toplevel.fa -i Zea_mays.AGPv4.34.gff3 -o cds.fa
minimap2 -x splice -t 10 -k 12 -a -p 0.4 -N 20 Sorghum_bicolor.Sorghum_bicolor_NCBIv3.dna.toplevel.fa cds.fa > cds.sam
minimap2 -x splice -t 10 -k 12 -a -p 0.4 -N 20 Zea_mays.AGPv4.dna.toplevel.fa cds.fa > ref.sam
Perform alignment using the proali
function:
anchorwave proali -i Zea_mays.AGPv4.34.gff3 -as cds.fa -r Zea_mays.AGPv4.dna.toplevel.fa -a cds.sam -ar ref.sam -s Sorghum_bicolor.Sorghum_bicolor_NCBIv3.dna.toplevel.fa -n anchors -R 1 -Q 2 -o alignment.maf -f alignment.f.maf -w 120000 -fa 80000 -fa2 120000 -fa3 80000
Options of the proali
function:
Usage: anchorwave proali -i refGffFile -r refGenome -a cds.sam -as cds.fa -ar ref.sam -s targetGenome -n outputAnchorFile -o output.maf -f output.fragmentation.maf
Options
-h produce help message
-i FILE reference GFF/GTF file
-r FILE reference genome sequence
-as FILE anchor sequence file. (output from the gff2seq command)
-a FILE sam file by mapping conserved sequence to query genome
-s FILE target genome sequence
-n FILE output anchors file
-o FILE output file in maf format
-f FILE output sequence alignment for each block in maf format
-t INT number of threads (default: 1)
-fa3 INT if the inter-anchor length is shorter than this value, stop trying to find new anchors (default: 10000)
-w INT sequence alignment window width (default: 30000)
-R INT reference genome maximum alignment coverage (default: 1)
-Q INT query genome maximum alignment coverage (default: 2)
-B INT mismatching penalty (default: -2)
-O1 INT open gap penalty (default: -3)
-E1 INT extend gap penalty (default: -2)
-O2 INT open gap penalty 2 (default: -80)
-E2 INT extend gap penalty 2 (default: -1)
-m INT minimum exon length to use (default: 20, should be identical with the setting of gff2seq function)
-mi DOUBLE minimum full-length CDS anchor hit similarity to use (default:0)
-mi2 DOUBLE minimum novel anchor hit similarity to use (default:0)
-e INT maximum expected copy number of each gene on each chromosome (default: 1)
This prevents using tandem duplicated genes to identify collinear block
-y DOUBLE minimal ratio of e+1 similarity to 1 similarity to drop an anchor (default: 0.6)
-ar FILE sam file by mapping conserved sequence to reference genome
this is used to improve the accurancy of anchors mapping
Following parameters are for identify collinear blocks
-d DOUBLE calculate IndelDistance (default: 3)
-OC DOUBLE chain open gap penalty (default: -0.03)
-EC DOUBLE chain extend gap penalty (default: -0.01)
-I DOUBLE minimum chain score (default: 2)
-D INT maximum gap size for chain (default: 25)
-ns do not search for new anchors (default: false)
-x use exon records instead of CDS from the GFF file
The maf-convert script could be used to reformat the MAF format file into several different formats.
For example, reformat maize B73 against Mo17 alignment into sam and bam format:
python2 maf-convert sam anchorwave.maf | sed 's/Zea_mays.AGPv4.dna.toplevel.fa.//g' | sed 's/Zm-Mo17-REFERENCE-CAU-1.0.fa.//g' | sed 's/[0-9]\+H//g' > anchorwave.sam
cat anchorwave.sam | samtools view -O BAM --reference Zea_mays.AGPv4.dna.toplevel.fa - | samtools sort - > anchorwave.bam
samtools index anchorwave.bam
If you could like to swap reference and query sequence, there is no need to re-run AnchorWave, which is computational cost.
The script anchorwave-maf-swap.py
under folder scripts
has been implemented for this aim:
cat anchorwave.maf | python2 anchorwave-maf-swap.py >anchorwave_swap.maf
maf-convert sam anchorwave_swap.maf | sed 's/Zea_mays.AGPv4.dna.toplevel.fa.//g' | sed 's/Zm-Mo17-REFERENCE-CAU-1.0.fa.//g' | sed 's/[0-9]\+H//g' > anchorwave_swap.sam
samtools view -O BAM --reference Zm-Mo17-REFERENCE-CAU-1.0.fa anchorwave_swap.sam | samtools sort - > anchorwave_swap.bam
A more detailed guideline for which algorithm to use, how to visualize the results and how to select parameters could be found at the walkthrough guidelines.
More examples and comparsion with other tools could be found at genome alignment.
The FAQ document is available at FAQ.md.
Bug report? Any question? Any suggestion? Any requirement? Want to collaborate?
Please feel free to send E-mail to [email protected].
This work is funded by NSF #1822330.
The manuscript is under preparation.