Apply cutoffs for domains

[apcamargo]

This is more of a question than a bug.

My understanding is that when the user set a cutoff, PyHMMER will not report hits that would have hit.included == False. However, the same is not applied at the domain level, so there are cases where the hit is reported because one of the domains passes the cutoffs, but domains that do not satisfy the cutoffs are also reported.

For example:

>P03865 1-315
MATTKLGNTKSASRAINYAEKRAEEKSGLNCDVDYAKSAFKQTRALYGKEDGIQAHTVIQSFKPGEVTPEQCNQLGLELA
EKIAPNHQVAVYTHTDKDHYHNHIVINSVDLETGKKYQSNKKQRDLVKKENDNICREHGLSVTERGIAKMRYTQAEKGIV
FDRDEYSWKDELRDLIENAKTHTSNLETFSEHLEEKGVGVKLRGETISYKPENENKWVRGRTLGSEYEKGAIDHEHERHQ
KQQREPEYADEFKINWDAVEQHTEQLKQRRVERAQETKQAHSKISSRDTRESENQRERAKGNNIRIERGDEGLSR

with open("hmmsearch_plasmid_hallmarks.tsv", "w") as fo:
    fo.write("query\thit\tstart\tend\tevalue\tdomain_score\tdomain_included\n")
    
    for p in Path(".").glob("*.faa"):
        with pyhmmer.easel.SequenceFile(p, digital=True) as seq_file:
            seqs = seq_file.read_block()

        with pyhmmer.plan7.HMMFile("Pfam-A.hmm") as hmm_file:
            for hits in pyhmmer.hmmsearch(hmm_file, seqs, bit_cutoffs="gathering"):
                for hit in hits:
                    for domain in hit.domains:
                        fout.write(
                            f"{hits.query_accession.decode()}\t"
                            f"{hit.name.decode()}\t"
                            f"{domain.env_from}\t"
                            f"{domain.env_to}\t"
                            f"{hit.evalue:.2E}\t"
                            f"{domain.score:.2f}\t"
                            f"{domain.included}\n"
                        )

Generates this output:

query        hit      start   end   evalue     domain_score   domain_included
PF03432.18   P03865   8       239   3.25E-90   288.01         True
PF03432.18   P03865   240     314   3.25E-90   -0.23          False

To be honest, I don't really know what determines if a domain is included or not, given that the cutoffs are applied to the hit. I expected that when a cutoff is applied, domains were also filtered.

[althonos]

A bit of background before I give you the easy answer 😁

So there are two kind of thresholds in HMMER, reporting and inclusion. When you use bitscore cutoffs with the pipeline, you actually change the inclusion thresholds, not the reporting thresholds. That's why you sometimes get hits that are not "included", but all hits should have hit.reported == True (otherwise they are not kept in the TopHits, which as its name suggests, only keep the top hits, not all hits).

The reason why there are different thresholds for hits and domains is to account for multi-domain hits, where sometimes a hit passes the inclusion thresholds because the combination of domains do, however individual domains do not pass the threshold (this is reported in HMMER's output, sometimes you'll see "no domain above threshold" or something like that).

Now in your example, to only iterate on included hits, and for each hit on included domains, you can use the included property to avoid a cascade of if/else:

for hit in hits.included:
    for domain in hit.domains.included:
        # ...

[apcamargo]

Thanks for the explanation, Martin! That makes sense.

However, I've seem at least a case where the hit was reported, but it had a single domain that wasn't. If having multiple domains is what makes hit with unreported domains being reported, how can a hit with a single unreported domain be reported?

Here's an example:

>A0A010SHV5 1-984
MLDITTISRQSLGKVVSYYADGADDYYAKDGGAMQWQGAGAEALGLSGEVEQARFRELLDGRISDSTKLMRTVKEADGKV
VSKERLGYDLTFSAPKGVSLQALVHGDASIIEAHDKAVAAAIREAERLSQARITVNKKTGTENTNNLVVAKFRHETSRAL
DPDLHTHAFVLNMTQRSDGEWRALKNDGVFNSSMFLGNVYKAELARELEKAGFQLRYERNGTFDLAHFSDEQIREFSSRS
QQIEAALAAKGLDRSTASYAEKNQAALATRDKKQGGIDREELRQVWLERSRALGIDYHSREWAGVGADAQGGRERNSAAT
PQIEKPLEYRADQVIEFAIKSLTERQAVIGQKELMDTALRHGYGALTIDDVRAGIERRAASGHLIREEPLYASQNPADGK
KGKAAEKARQEAPQLSRKEWVANLVRAGKSRSEAARLVDEGIRTGRLRQGENRFTTHIAQKREREVLQIERMGRGTVEPR
ISKEAAEAFLADRGLKAEQQASVMRIARTQNQFIGVQGFAGVGKSYMTVAAKDLLEANGYRVTSLAPYGSQVKALQAEGL
EARTLQSFLKARDKKIDSNTVVFIDEAGVIPARQMHEAMKTIEAAGARVVFLGDVSQTKAIEAGKPFEQLMKAGMETSRL
TDIQRQKDPQLLEAVKLAAEGKAKQSLPLVNEIHEIKEDGARYQAIVDAYASMPKAERDQALIITGTNASRIEINEGVRE
ALGLKGQGAEYPLLNRLDTTQAERRHSKYYGKGSIVVPEVDYKNGLQRGVQYVVLDTGPGNKLTVRGPEGETLQFTPARC
TKLSVYSVERTELAPGDQVKITRNDAEKDLANGDRFVVKEIHSDRVVLEGDKGRQVELDTASAMFVGLAYASTVHSAQGL
TCDKVLINLETQSRTTAKDVYYVAISRARHAAEIFTDNRQKLGAAISRLNAKAGALDIKQLQRHALERKGHDQAGKQKDA
AQKQQQGQQLPTKQPKEKGRAFGL

print("pfam\tprotein\tdomain_number\thit_included\tdomain_included")
with pyhmmer.easel.SequenceFile("test.faa", digital=True) as seq_file:
    seqs = seq_file.read_block()
    with pyhmmer.plan7.HMMFile("Pfam-A.hmm") as hmm_file:
        for hits in pyhmmer.hmmsearch(hmm_file, seqs, bit_cutoffs="gathering"):
            for hit in hits:
                for i, domain in enumerate(hit.domains, 1):
                    print(
                        f"{hit.name.decode()}\t{hits.query_accession.decode()}\t"
                        f"{i}\t{hit.included}\t{domain.included}"
                    )

Output:

pfam          protein       domain_number    hit_included    domain_included
A0A010SHV5    PF13245.10                1            True               True
A0A010SHV5    PF13245.10                2            True              False
A0A010SHV5    PF13604.10                1            True               True
A0A010SHV5    PF08751.15                1            True               True
A0A010SHV5    PF08751.15                2            True              False
A0A010SHV5    PF13538.10                1            True              False
A0A010SHV5    PF01443.22                1            True               True
A0A010SHV5    PF01443.22                2            True               True

As you can see, PF13538.10 has a single domain that was not included, but the hit is. Why is that?

[raufs]

Also interested here!

Maybe useful - the domain.score (26.71) for PF1538.10 seems to be slightly deflated relative to the hit.score (27.91), with the gathering threshold for the profile being 27.4.

This behavior is consistent with using HMMER3 though, where similarly domain.score of 26.7 and the full sequence score is 27.9.

hmmsearch -E 10.0 --domtblout test_hmmer3.txt Pfam-A.hmm test.faa

Results in test_hmmer3.txt:

#                                                                            --- full sequence --- -------------- this domain -------------   hmm coord   ali coord   env coord
# target name        accession   tlen query name           accession   qlen   E-value  score  bias   #  of  c-Evalue  i-Evalue  score  bias  from    to  from    to  from    to  acc description of target
#------------------- ---------- ----- -------------------- ---------- ----- --------- ------ ----- --- --- --------- --------- ------ ----- ----- ----- ----- ----- ----- ----- ---- ---------------------
A0A010SHV5           -            984 AAA_19               PF13245.10   135   3.4e-12   32.9   0.8   1   2     1e-11     1e-11   31.3   0.3     8   130   508   618   499   623 0.89 1-984
A0A010SHV5           -            984 AAA_19               PF13245.10   135   3.4e-12   32.9   0.8   2   2      0.79      0.79   -3.9   0.0    32    62   687   718   684   731 0.70 1-984
A0A010SHV5           -            984 AAA_30               PF13604.10   191   1.5e-48  151.2   5.0   1   1   4.3e-48   4.3e-48  149.7   5.0     2   187   495   677   494   680 0.95 1-984
A0A010SHV5           -            984 DUF2563              PF10817.12   104   1.1e-05   12.3   1.8   1   2      0.74      0.74   -3.2   0.1    40    59   106   125    99   142 0.57 1-984
A0A010SHV5           -            984 DUF2563              PF10817.12   104   1.1e-05   12.3   1.8   2   2     2e-05     2e-05   11.4   0.3    35    97   861   928   850   935 0.79 1-984
A0A010SHV5           -            984 Herpes_Helicase      PF02689.18   808   1.9e-07   15.4   0.1   1   2     0.036     0.036   -2.0   0.0   171   206   556   591   549   597 0.71 1-984
A0A010SHV5           -            984 Herpes_Helicase      PF02689.18   808   1.9e-07   15.4   0.1   2   2   5.6e-07   5.6e-07   13.9   0.1   731   777   870   917   847   927 0.83 1-984
A0A010SHV5           -            984 PhoH                 PF02562.20   205   1.9e-09   23.4   1.3   1   3   2.2e-05   2.2e-05   10.2   0.1     6    39   497   531   493   543 0.77 1-984
A0A010SHV5           -            984 PhoH                 PF02562.20   205   1.9e-09   23.4   1.3   2   3   1.1e-05   1.1e-05   11.2   0.1   120   162   579   622   561   646 0.85 1-984
A0A010SHV5           -            984 PhoH                 PF02562.20   205   1.9e-09   23.4   1.3   3   3      0.28      0.28   -3.3   0.0    81   126   731   776   718   777 0.83 1-984
A0A010SHV5           -            984 PIF1                 PF05970.18   223   2.3e-06   13.4   0.3   1   2     2e-05     2e-05   10.4   0.1     5    66   498   557   494   587 0.74 1-984
A0A010SHV5           -            984 PIF1                 PF05970.18   223   2.3e-06   13.4   0.3   2   2     0.023     0.023    0.4   0.0   131   174   597   640   590   675 0.73 1-984
A0A010SHV5           -            984 SH3_13               PF18335.5     65   4.5e-07   15.9   0.0   1   1   1.4e-06   1.4e-06   14.3   0.0    23    57   809   844   794   852 0.79 1-984
A0A010SHV5           -            984 TrwC                 PF08751.15   283   9.4e-91  290.3   5.0   1   2   9.4e-91   9.4e-91  290.3   5.0     3   282    14   292    12   293 0.95 1-984
A0A010SHV5           -            984 TrwC                 PF08751.15   283   9.4e-91  290.3   5.0   2   2     0.072     0.072   -1.3   1.3   194   231   915   953   891   981 0.55 1-984
A0A010SHV5           -            984 UvrD_C_2             PF13538.10    52   8.5e-11   27.9   0.3   1   1     2e-10     2e-10   26.7   0.3     1    48   868   911   868   914 0.89 1-984
A0A010SHV5           -            984 Viral_helicase1      PF01443.22   228   2.4e-17   49.6   0.5   1   2     1e-08     1e-08   21.4   0.0     3   102   516   623   514   690 0.80 1-984
A0A010SHV5           -            984 Viral_helicase1      PF01443.22   228   2.4e-17   49.6   0.5   2   2   4.1e-10   4.1e-10   25.9   0.2   160   228   852   916   791   916 0.62 1-984

[althonos]

Hi both, there may have been a bug in TopHits.merge with bitscore cutoffs in particular which caused flags to be reset incorrectly. If the issue persists, I'm still not sure what caused it and I will have a second look.

[raufs]

Thank you Martin! There might not be a bug here - or I think PyHMMER acts the same as HMMER at least. The --tblout and --domtblout arguments function more akin to different modes rather than simple output formats when score thresholds are used. Specifically, I think that while Pfam pHMMs have two scores for thresholding type (type being gathering, noise, trusted, and values corresponding to reporting/inclusion at the per-sequence and per-domain levels, respectively) - because the domain-specific score and full-sequence score can differ - what is shown in the output will depend on which "output" mode is requested and which cutoff is applied. This seems to be the case based on basic testing using hmmsearch (HMMER v3.4) and the example Antonio shared. However, from pg 100 of the HMMER 3.4 (Aug 2023) documentation, it might be that both should be getting applied even when --tblout is used:

Each thresholding option has two scores: the per-sequence threshold and the per-domain threshold . These act as if -T --incT --domT --incdomT has been applied specifically using each model’s curated thresholds.

and

--cut_ga: Use the GA (gathering) bit scores in the model to set per-sequence (GA1) and per-domain (GA2) reporting and inclusion thresholds. GA thresholds are generally considered to be the reliable curated thresholds defining family membership; for example, in Pfam, these thresholds define what gets included in Pfam Full alignments based on searches with Pfam Seed models.

Practically, given the current setup, I think that if you are annotating using domain-level pHMMs, such as Pfam domains, it will often probably make the most sense to use the --domtblout option in hmmsearch/hmmscan or go off the hit.domains.included listing with PyHMMER. This is generally more conservative since I think sequence scores are larger than or equal to domain scores and GA1 is often equal to GA2. From looking at Pfam-A v38.0, it seems that for many of the 25,545 pHMMs (98.1%), GA1 == GA2, but for 476 pHMMs GA1 > GA2 and for 6 pHMMs GA2 > GA1.

Sorry for all the edits here, tried my best to use the same wording as hmmer/pyhmmer, look forward to hearing both your thoughts/perspectives here.