6
$\begingroup$

I want to perform a multi-sequence alignment on 12 samples that clustered based on cgMLST. Ultimate goal is to find out whether they differ by the presence of a few genes.

I have performed multi-sequence alignment in the past using MAFFT, the main difference being that, back then, I was aligning protein sequences, and now I want to align multiple denovo assembled genomes, with multiple contigs each

I first had a look at MAFFT, the program that I used in the past for protein multi-sequence alignment, though I noticed that it considers each contig in a denovo assembly as independent sequence, which is not what I want

Then I had a look at Mauve, but was discouraged by installation trouble-shooting. Also I read on their webpage:

When an individual file contains several sequence entries, they will be concatenated and the whole concatenated sequence will be aligned to the sequences in the other files. This behavior allows multi-chromosomal genomes to be aligned by including all chromosomes in a single sequence file. Similarly, incomplete genomes that consist of several sequence contigs can be aligned, though beware that incorrectly ordered contigs will appear as genome rearrangements in the Mauve viewer.

So I guess concatenating all contigs in a single fasta file is not an option, as artefacts (genome rearrangements) might appear

It turns out that MAFFT has also an experimental feature to multi-align fastq reads, which sounds attractive to me. Therefore I would like to ask if anyone can recommend this new/ experimental feature. Also if this option is available on the command-line tool, as the mafft help is quite limited (it does not specify the format of input and output files):

  MAFFT v7.505 (2022/Apr/10)
  https://mafft.cbrc.jp/alignment/software/
  MBE 30:772-780 (2013), NAR 30:3059-3066 (2002)
------------------------------------------------------------------------------
High speed:
  % mafft in > out
  % mafft --retree 1 in > out (fast)

High accuracy (for <~200 sequences x <~2,000 aa/nt):
  % mafft --maxiterate 1000 --localpair  in > out (% linsi in > out is also ok)
  % mafft --maxiterate 1000 --genafpair  in > out (% einsi in > out)
  % mafft --maxiterate 1000 --globalpair in > out (% ginsi in > out)

If unsure which option to use:
  % mafft --auto in > out

--op # :         Gap opening penalty, default: 1.53
--ep # :         Offset (works like gap extension penalty), default: 0.0
--maxiterate # : Maximum number of iterative refinement, default: 0
--clustalout :   Output: clustal format, default: fasta
--reorder :      Outorder: aligned, default: input order
--quiet :        Do not report progress
--thread # :     Number of threads (if unsure, --thread -1)
--dash :         Add structural information (Rozewicki et al, submitted)

Any other recommendation on another strategy to align multiple denovo assembled genomes with multiple contigs each are welcome

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ Thanks its a very good question, for alignments MAFFT in very recent history was the optimal aligner on 2 criteria. $\endgroup$
    – M__ ♦
    Commented Jun 29, 2022 at 13:37
  • $\begingroup$ Were you able to resolve your issue? And which tool did you use? I have a similar problem, though the number of sequences I have to align are in the thousands. $\endgroup$
    – pseudomonas
    Commented Nov 16, 2022 at 22:24
  • $\begingroup$ @Raghuram in the end I just loaded the .bam files into artemis and inspected the regions of interest. Looking at the depth, I could inspect the presence/ absence of the genes I wanted to detect. $\endgroup$
    – BCArg
    Commented Nov 17, 2022 at 9:39

4 Answers 4

Reset to default
2
$\begingroup$

MUMmer4 is a versatile alignment tool for DNA and protein sequences. It supports one reference genome and up to 32 query genomes. MUMmer4 will align every contig in each genome to the reference genome.

Improve this answer
$\endgroup$
4
  • $\begingroup$ do you concatenate all query genomes in one fasta file? If so, how does MUMmer identify which contig belongs to which genome? Based on a header identifier, for example? $\endgroup$
    – BCArg
    Commented Jun 29, 2022 at 14:45
  • $\begingroup$ @BCArg I did not concatenate all query genomes in one fasta file. $\endgroup$
    – Forrest Vigor
    Commented Jun 30, 2022 at 7:21
  • $\begingroup$ @ForrestVigor does it support anything like nucmer -p <prefix> ref.fa qry1.fa qry2.fa ... qryn.fa ? on the help page it reads <query> specifies the multi-FastA sequence file that contains the query sequences, to be aligned with the references $\endgroup$
    – BCArg
    Commented Jun 30, 2022 at 9:54
  • 2
    $\begingroup$ @BCArg Yes. You can use mummer ref.fa qry1.fa qry2.fa ... $\endgroup$
    – Forrest Vigor
    Commented Jun 30, 2022 at 13:16
2
$\begingroup$

If you have a reference genome (or are willing to designate one of your de novo assemblies a reference), you may find QUAST helpful.

QUAST will perform an alignment of all genomes against the reference genome and use that as the basis for comparisons. I believe that under the hood it is running minimap2 for this purpose.

It includes a large-scale alignment viewer ("Icarus") to investigate possible differences within contigs, computes high-level statistics on genome completeness/etc, and has other handy features.

You can see the most recent QUAST paper or the manual for more information.

Improve this answer
$\endgroup$
3
  • $\begingroup$ thanks for the comment, though I am looking for a program to perform the multi-sequence alignment, not to visualise the alignment, for which purpose QUAST is indeed excellent $\endgroup$
    – BCArg
    Commented Jul 4, 2022 at 6:49
  • 1
    $\begingroup$ @BCArg QUAST will generate the alignment. Updating answer accordingly. $\endgroup$
    – Maximilian Press
    Commented Jul 4, 2022 at 20:03
  • $\begingroup$ @Maximillian Press you are right, and obviously it makes sense that Quast needs to perform the alignment to evaluate/ compare two assemblies. Will explore it further $\endgroup$
    – BCArg
    Commented Jul 5, 2022 at 11:30
2
$\begingroup$

in the end, I just loaded the .bam files into artemis and, by inspecting the depth (heat-map), I could check which samples had the genes I was looking for:

enter image description here

Improve this answer
$\endgroup$
1
$\begingroup$

It won't work on more than 3 genomes at a time, but SynMap3D from the Comparative Genomics webserver (https://genomevolution.org/coge/) should do the trick. I'm mentioning it here because the advantage, for me (in contrast to, say, Mauve), was user-friendliness.

Small note: Be mindful of what you name your uploaded samples. You can set your sequence not to be publicly available, but if you give it a new name, as opposed to a species name already in the database, that will be visible.

Improve this answer
$\endgroup$
1
  • $\begingroup$ thanks, though I am looking for a MS aligner that supports > 3 sequences $\endgroup$
    – BCArg
    Commented Jun 29, 2022 at 14:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.