ClustalOmega alignment

Question

I am performing multiple sequence alignment using Clustal Omega for 600 RNA covid-sequences, i.e. cDNA on Genbank, with input characters around 30574 characters for each sequence.

I am running it on windows cmd. I defined parameters in input of maxseqle = 37000, however, the output each time gives duplicated length around 75000 charcters per each sequence.

clustalo.exe -i Allseq.fasta --is-profile --use-kimura --seqtype DNA --maxseqlen 37000 --threads 8 -o myclustalv3.fasta

How can I solve this problem, what will I define in input parameters. This problem occurs with large no of sequences 600 sequences. In test of 40 seq, it gave me normal output length around 32000 characters but when align all 600 sequences together, it gave this duplicated length results.

Answer 1

Have you had a look at your alignment? I personally use Jalview but any alignment viewer should work. Colour > Percentage Identity (or Nucleotide) should help with visualization.

Does the alignment look reasonable (and gaps are just inserted at the end for some reason), or is it pretty messy? If it looks messy, is it just a few outlier sequences that cause the problem, or does it look like you have two main groups of sequences that align between themselves but not with the other? In the latter case, could it be possible that some of your initial sequences were inputed in the reverse orientation? I'm not sure about the command line implementation of Clustal Omega, but I know that at least the versions I use don't automatically reverse-complement sequences to be aligned.

Answer 2

Firstly its good you've done some proofing of your output wherein you spotted,

75000

The alignment output is approximately double the genome size, just over double the input.

An off-question point, I assume the computer you are working on has at least 8-cores because thats what you are requesting.

Clustal omega is a very solid algorhithm and has been around for a long time in its current omega version. The error you are encountering is an input sequence bug. @Laura's suggest is correct and is formal practice: you must manually verify bioinformatics output and in this case an alignment viewer is the correct proceedure.

Thereafter I suspect you would find the answer, such as omitting carriage returns in the fasta format, seeing taxa names caught up post 30K.

If you run,

wc ./myclustalv3.fasta # output file
wc ./Allseq.fasta  # input file
wc ./testgenome.fasta

You will probably find that the number lines per file does not agree between the input and output format OR the number of lines of both input and output is <<600 or <<1200 (600 x2 ), i.e. the expected number of sequences. This points to parsing errors within the input data. Fasta has a line for the sequence identifier and a line for the sequence. It is using > for the sequence identifier, but will use a carriage return for sequence.

Again manually assessing an alignment is baseline proceedure regardless of output.

An example (guess) about what might have happened is that every other input has omitted ">", omega will then see the first ">" then read everything until it sees another ">" as sequence data. Thus its doubled the sequence length and dragged a sequence identifier exactly in the centre of the doubled up genomes.

---

An issue with fasta is there is no internal proofing of the format. In other formats the number of sequences and their length is specified within the file, so it is easy to internally verfiy the data prior a calculation. Code-wise it is known as 'setting' and you can't do that with fasta.