Closed Augustin-Zidek closed 4 months ago
Those two alignments don't appear to contain the same sequence data, so I wouldn't expect them to give the same profile. They aren't just the same alignment in two different formats.
Where you say "Input a3m created by removing all columns where the query had a gap and also removing deletions (lower-case letters) in non-query sequences so that everything has exactly the same length", you mean you're deleting residues from the (Stockholm version of the) alignment? That's changing the sequences; the two alignments are different after those steps are taken.
I've only glanced at this quickly, so if I've misunderstood something, let me know.
Thanks for the blazingly fast reply!
Where you say "Input a3m created by removing all columns where the query had a gap and also removing deletions (lower-case letters) in non-query sequences so that everything has exactly the same length", you mean you're deleting residues from the (Stockholm version of the) alignment? That's changing the sequences; the two alignments are different after those steps are taken.
Yes, you are right that the inputs are not exactly the same. However, if I use the original (uncleaned) a3m:
>query
FHQECSLQSCTQHQPYVVDDPCPIHFYSKWYIRVGARKSAPLIELCVDEAGSKSPIQYIDIGNYTVSCLPFTINCQEPKLGSLVVRCSFYEDFLEYHDVRVVLDFI
>UniRef90_A0A2R3SUZ9/16-121 [subseq from] Uncharacterized protein n=1 Tax=Bat SARS-like coronavirus TaxID=1508227 RepID=A0A2R3SUZ9_CVHSA
FHQECSLQSCAQHQPYVVDDPCPIHFYSRWYIRVGARKSAPLIELCVDEVGSKSPIQYIDIGNYTVSCSPFTINCQEPKLGSLVVRCSYYEDFLEYHDIRVVLDFI
>UniRef90_A0A023PSW4/17-121 [subseq from] Protein 10b n=1 Tax=Rhinolophus affinis coronavirus TaxID=1487703 RepID=A0A023PSW4_CVHSA
-HKECTIQECCENQPYILEDPCPIHYYSDWYLKIGPRKSARLLQLCAGEYGKRLPVQYEKLGNYTINCEPFEINCQTPPVGSLIVRCSYDYDFIEYHDVRVVLDFI
>UniRef90_U5WLD0/17-121 [subseq from] Uncharacterized protein n=11 Tax=Betacoronavirus TaxID=694002 RepID=U5WLD0_CVHSA
-HKECSIQECCENQPFQLEDPCPIHYYSDWFVKIGPRKSARLVQLCAGEYGHRVPIHYEMFGNYTISCEPLEINCQNPPVGSLIVRCSYDVDFMEYHDVRVVLDFI
>UniRef90_Q0Q469/17-121 [subseq from] Non-structural protein 8 n=17 Tax=Betacoronavirus TaxID=694002 RepID=NS8_BC279
-HKECSIQECCENQPYQIEDPCPIHYYSDWFIKIGSRKSARLVQLCEGDYGKRIPIHYEMFGNYTISCEPLEINCQAPPVGSLIVRCSYDYDFVEHHDVRVVLDFI
>UniRef90_A0A0U1UYY4/17-121 [subseq from] Uncharacterized protein n=1 Tax=BtRs-BetaCoV/HuB2013 TaxID=1503302 RepID=A0A0U1UYY4_CVHSA
-HKECSIQECCEIQPSQIEDPSPIHYYSDWFIQIGYRKSARLVQLCEGDYGKRIPIHYEMFGNYSMYCEPLEINCQAPPVGRLIVRWLHDYESAEHHDVRVVLDFI
>UniRef90_D5HJV2/17-51 [subseq from] Uncharacterized protein n=1 Tax=Bat SARS coronavirus HKU3-8 TaxID=742005 RepID=D5HJV2_BCHK3
-HKECSIQECCENQPYQIEDPCPIHYYSDWFIKIGS----------------------------------------------------------------------
>UniRef90_Q3ZTE5/20-117 [subseq from] Orf10 n=10 Tax=Human SARS coronavirus TaxID=694009 RepID=Q3ZTE5_CVHSA
------VQRCASNKPHVLEDPCPTGYQPEWNIRYKTRGNTySTAWLCA--LGKVLPFHRW--HTMVQTCTPnVTINCQDPAGGALIARCWYLHEghqTAAFRDVLVVL---
>UniRef90_A0A0K1Z007/20-103 [subseq from] Uncharacterized protein n=6 Tax=Human SARS coronavirus TaxID=694009 RepID=A0A0K1Z007_CVHSA
------VQRCVSNTPYVLENPCPTGYRPEWNIRYNTRGNTyNTARLCA--LGK--VLSFHRWHTMVQACTPnITINCQDPVGGALVARCWYFHK--------------
I get the following error:
Alignment input parse error:
sequence UniRef90_Q3ZTE5/20-117 has alen 111; expected 106
while reading aligned FASTA file orig_input.a3m
at or near line 17
Our format autodetector is detecting the file as aligned FASTA. If you give it --informat a2m
or make the file have a suffix of .a2m
instead of .a3m
, it should work fine.
more detail: It's not possible for our format detector to accurately distinguish aligned FASTA from A2M (or A3M) formats, because we have a requirement that we must not misread either the alignment or any associated annotation. There are ambiguous cases where the alignment is valid both as AFA or A2M/A3M, but A2M/A3M also implies annotation of consensus columns and AFA does not.
So we only autodetect AFA format, and are unable to autodetect A2M/A3M. A2M/A3M requires some sort of other affirmative hint from the user; either the file suffix is .a2m
or --informat a2m
is provided as a commandline option.
A3M format is the same as "dotless" A2M format, and unfortunately we're a little behind the times (our parsers predate the naming of A3M, I think) and we currently call both formats A2M. I've made a note to update the code so we can see the .a3m
suffix too, which would also fix the problem you're seeing.
Oh, I see, thank you very much! --informat a2m
indeed addressed the error I was seeing.
However, I think we will have to stick with using Stockholm nevertheless as a2m/a3m doesn't allow us to define consensus columns as far as I know.
I.e. we need the consensus columns to be all columns that don't have a gap in the query (the first sequence in the input MSA).
> query
AB--C
> seq 1
ABMMC
> seq 2
ABMMCdd
> seq 3
A--MkC
We can easily achieve that by setting GC RF=xx..x
in Stockholm.
Is there a way to do this using a2m/afa?
Yes, a2m/a3m use uppercase and '-' to denote consensus positions, whereas lowercase and '.' are insertions or gaps aligned to insertions. Each aligned sequence must have the same number of uppercase and '-' characters; that number is the number of consensus positions in the profile.
Our parser picks up a reference annotation line from a2m/a3m format based on this.
Thanks again, closing this issue.
Hmmbuild gives different results depending on the input format. For Stockholm, it gives the output of the same length as the aligned sequences. For A3M, it sometimes does not.
In order to get the correctly aligned result, we therefore have to convert a3m -> sto - adding the
#=GC RF
line while doing so - then run Hmmbuild with the converted Stockholm file.However, in our convertor, the
#=GC RF
line is generated trivially:x
in every position the first sequence of the a3m (i.e. the query sequence in the previous MSA search step) has an amino acid.
in every column where the first sequence of the a3m has a gap (-
).I realise we are seeing this issue most likely because we are not specifying the reference annotation when using a3m. Is there a way to specify the reference annotation in the a3m case? If not, would it be possible to add that option so that the conversion is not needed?
Reproduction
Input a3m created by removing all columns where the query had a gap and also removing deletions (lower-case letters) in non-query sequences so that everything has exactly the same length (
input.a3m
):Converted input Stockholm (
input.sto
):We then run the following pairs of commands:
hmmbuild --amino output_a3m.hmm input.a3m
hmmbuild --amino output_sto.hmm input.sto
hmmbuild --amino --fast output_fast_a3m.hmm input.a3m
hmmbuild --amino --fast output_fast_sto.hmm input.sto
hmmbuild --amino --hand output_hand_a3m.hmm input.a3m
- fails with "Error: build failed: Alignment input has no reference annotation line"hmmbuild --amino --hand output_hand_sto.hmm input.sto
- worksThe respective pairs for cases 1 and 2 don't match. E.g. running
diff --side-by-side --color --suppress-common-lines output_{a3m,sto}.hmm
will produce:Happy to provide extra details!