This is the first workflow guide and covers the whole of one analysis: from an annotation table to a set of calls, and from a call back down to the genes that produced it.
What gifter answers, and what it does not
A GIFT is a genome-inferred functional trait: a biologically meaningful capability whose genomic support is decided by an explicit, curated model. gifter answers one question per GIFT:
Do the markers observed in this genome support at least one complete curated implementation of this capability?
It does not answer whether the capability is expressed, active, used, or consequential. Keep that separation in mind from the start, because every number later in these tutorials inherits it.
Step 1: get your annotations into three columns
Starting from FASTA? giftag searches the profiles used by gifter and writes
giftag_markers.tsvin the format this guide expects. Its handoff guide shows the exact R call. Review giftag’s searchable-marker ledger before interpreting an unsupported GIFT as evidence of absence.
gifter takes a plain data frame. Whatever tool produced your annotations — kofamscan, dbCAN, eggNOG-mapper, DRAM, bakta — the job is to reshape its output into one row per observed marker:
| column | meaning |
|---|---|
gene_id |
your identifier for the gene. Optional, but it is what makes a call traceable back to a locus |
namespace |
which accession system: KO, EC,
PFAM, TIGRFAM, NCBIFAM,
CAZY
|
accession |
the accession itself |
Long form matters: a gene hit by both a KO and a CAZy family is two rows, not one row with two columns.
Here is an illustrative genome. It stands for an anaerobic gut bacterium that degrades xylan, and it is a fixture chosen to exercise the model rather than annotation output from a sequenced organism.
gut_anaerobe <- data.frame(
gene_id = c(
# xylan backbone cleavage and debranching
"g0101", "g0102", "g0103",
# xylose transport and catabolism
"g0210", "g0211", "g0212", "g0330", "g0331",
# purine biosynthesis
"g0450", "g0451", "g0452", "g0453", "g0454", "g0455", "g0456",
# chemotaxis
"g0620", "g0621", "g0622", "g0623", "g0624", "g0625", "g0626",
# restriction-modification
"g0810", "g0811", "g0812",
# amino acid biosynthesis
"g0910", "g0911", "g0912", "g0913", "g0914",
# housekeeping genes gifter does not use
"g1001", "g1002", "g1003"
),
namespace = c(rep("CAZY", 3), rep("KO", 30)),
accession = c(
"GH11", "GH39", "CE1",
"K10543", "K10544", "K10545", "K01805", "K00854",
"K00764", "K01945", "K00601", "K01952", "K01587", "K01756", "K00602",
"K03412", "K00575", "K03407", "K03408", "K03776", "K03413", "K03410",
"K03427", "K01153", "K01154",
"K00812", "K01915", "K00600", "K00826", "K00259",
"K02945", "K02863", "K01873"
),
stringsAsFactors = FALSE
)
nrow(gut_anaerobe)
#> [1] 33If you only have a bare vector of accessions, gifter will infer the
namespace for KO, EC, Pfam, TIGRFAM, NCBIfam and CAZy patterns.
Supplying namespace explicitly is safer.
NCBIFAM is the namespace for the profile library PGAP,
AMRFinderPlus, bakta and InterProScan report, and its accessions carry a
release version suffix: NF040708.3,
TIGR04545.1. Keep the suffix. A bare
TIGR03948 is read as the unversioned TIGRFAM
identity instead, which is a different marker, and a bare
NF040708 matches nothing. A KEGG-only annotation produces
no NCBIfam accession at all; every component that accepts one also
accepts a KO, with one deliberate exception —
siroheme_to_heme_b needs NF040708.3 and
NF040707.3 because no KEGG accession separates the two
subunits of its first enzyme.
Two questions gifter asks before it evaluates
Two input mistakes produce a result that is well formed, internally
consistent, and about something other than what you asked. Neither is
visible in the calls themselves, so evaluate_gifts()
settles both before evaluating rather than reporting them
afterwards.
The first is the gene column. If your table has no
gene_id column, evaluate_gifts() proposes the
first column that is neither namespace nor
accession and waits for your approval; without an answer —
in a script, a knitted document, a batch job — it fails rather than
choosing. Adopting the wrong column would leave every call intact while
attaching the evidence to the wrong loci, which is worse than an error.
Answer it in advance:
# The gene column under its own name.
evaluate_gifts(annotations, gene_id = "locus_tag")
# Approve the first column, whatever it is called.
evaluate_gifts(annotations, gene_id = TRUE)
# No gene column: number the markers and let the trace stop at the marker.
evaluate_gifts(annotations, gene_id = FALSE)The second is scale. evaluate_gifts() evaluates
one genome. Markers pooled from several genomes
complete routes that no single genome encodes, so a table carrying more
than max_genes distinct gene identifiers (5,000 by default,
more than most bacterial genomes annotate) is questioned before it is
treated as one genome. If the input really is one large genome, raise or
disable the threshold with max_genes = Inf. If it is a
collection, split it by genome — the third tutorial builds a community
from separately evaluated genomes.
Step 2: see which markers gifter can use
Before evaluating anything, check how much of your annotation the
reference database actually recognises. map_markers() keeps
every observation, including the ones nothing matched, so marker use is
auditable rather than silent.
mapped <- map_markers(gut_anaerobe)
# One row per (observation, curated component) pair, plus unmatched rows.
table(mapped$matched)
#>
#> FALSE TRUE
#> 3 40The unmatched markers are not a problem — a genome carries thousands of genes and gifter curates a specific catalogue — but you should know which they are:
unique(mapped$accession[!mapped$matched])
#> [1] "K02945" "K02863" "K01873"Those three are ribosomal and tRNA-synthetase orthologues. They are real genes; gifter simply makes no capability claim that rests on them.
When you expect an annotation source to cover a curated marker namespace, inspect unmatched accessions within that namespace:
unmatched <- unique(mapped[!mapped$matched, c("namespace", "accession")])
split(unmatched$accession, unmatched$namespace)
#> $KO
#> [1] "K02945" "K02863" "K01873"Many unmatched annotations are ordinary: gifter curates a specific
set of capabilities, not every function in an annotation file. A run
where expected CAZY eCAMI clusters are consistently
unmatched is different. Their identifiers can be renumbered between
dbCAN releases, so compare your annotation release with the pinned
releases in SOURCES.md
before reading those unsupported GIFTs as absences. The same check
applies to versioned NCBIFAM accessions. Keep the release
identity with your marker table so the audit remains reproducible.
Step 3: evaluate
result <- evaluate_gifts(gut_anaerobe)
result
#> <gifter_genome>
#> # A tibble: 163 × 6
#> gift_id gift_type complete number_of_complete_i…¹ best_implementation
#> <chr> <chr> <lgl> <int> <chr>
#> 1 acetate_interc… metabolic FALSE 0 ACETATE_PTA_ACKA
#> 2 acetoin_format… metabolic FALSE 0 ACETOIN_ALS_BUDA
#> 3 acetyl_coa_to_… metabolic FALSE 0 TCA1_CS_ACN
#> 4 adenylate_bios… metabolic FALSE 0 AMP_ADENYLOSUCCINA…
#> 5 alanine_biosyn… metabolic TRUE 1 ALA_DEHYDROGENASE
#> 6 allantoin_degr… metabolic FALSE 0 ALLANTOIN_ALLANTOI…
#> 7 ammonium_assim… metabolic FALSE 0 AMMONIUM_GDH_NAD
#> 8 anthranilate_d… metabolic FALSE 0 ANTHRANILATE_TO_CA…
#> 9 arabinose_degr… metabolic FALSE 0 ARA_ISOMERASE
#> 10 arabinose_upta… metabolic FALSE 0 ARA_UPTAKE_ABC
#> # ℹ 153 more rows
#> # ℹ abbreviated name: ¹number_of_complete_implementations
#> # ℹ 1 more variable: minimum_missing_requirements <int>Every curated GIFT is evaluated, so the result always has the same number of rows. That matters later: a genome’s negative calls are part of the answer, not an absence of one.
nrow(result$gifts)
#> [1] 163The calls that came out complete
complete <- result$gifts[result$gifts$complete, ]
complete[, c("gift_id", "gift_type", "evidence_confidence")]
#> # A tibble: 10 × 3
#> gift_id gift_type evidence_confidence
#> <chr> <chr> <chr>
#> 1 alanine_biosynthesis metabolic curated
#> 2 aspartate_biosynthesis metabolic curated
#> 3 chemotaxis_signal_transduction regulatory curated
#> 4 glutamine_biosynthesis metabolic curated
#> 5 glycine_biosynthesis metabolic curated
#> 6 type_i_restriction_modification defense curated
#> 7 valine_biosynthesis metabolic curated
#> 8 xylan_degradation metabolic ambiguous
#> 9 xylose_degradation_isomerase metabolic curated
#> 10 xylose_uptake_abc metabolic putativeRead the gift_type column. It tells you which
completeness model produced the call, and therefore what the call is
allowed to mean:
- metabolic — at least one complete curated enzymatic route between the declared boundary molecules;
- regulatory — at least one complete curated circuit that senses a signal and executes a response;
- defense — at least one complete curated defense mechanism;
- structural — at least one complete curated architecture of a cellular structure.
Note evidence_confidence on
xylan_degradation: ambiguous. The call
rests on CAZy family-level accessions, and a glycoside hydrolase family
is polyspecific. gifter reports the call at the confidence its weakest
accepted marker allows rather than presenting it as equal to a call
built on characterised orthologues.
Step 4: read the incomplete calls, because they are informative
An incomplete call is not a blank. gifter reports the closest implementation and exactly what it was missing, so you can tell a capability that is one step short from one that is entirely absent.
gifts <- result$gifts
near <- gifts[!gifts$complete & gifts$minimum_missing_requirements == 1, ]
near <- near[order(-near$completeness_score), ]
head(near[, c("gift_id", "gift_type", "completeness_score")], 5)
#> # A tibble: 5 × 3
#> gift_id gift_type completeness_score
#> <chr> <chr> <dbl>
#> 1 purine_core_biosynthesis metabolic 0.9
#> 2 aspartate_chemoreception regulatory 0.75
#> 3 serine_chemoreception regulatory 0.75
#> 4 adenylate_biosynthesis metabolic 0.5
#> 5 acetoin_formation metabolic 0Take the top one:
purine <- gifts[gifts$gift_id == "purine_core_biosynthesis", ]
purine$best_route
#> [1] "PCR_FOLATE_DIRECT_FOLATE"
purine$minimum_missing_reactions
#> [1] 1
unlist(purine$missing_reactions_best_route)
#> [1] "RHEA:23032"
length(unlist(purine$supporting_reactions))
#> [1] 9Nine of the ten reactions of the closest route are supported and one
is not. That is a specific, checkable claim — go and look at whether
RHEA:23032 is genuinely absent or whether your annotation
missed it — and it is a very different situation from a capability with
no evidence at all.
Note what gifter did not do: it did not call the GIFT complete because 90% of its genes were present. Completeness comes from route logic. The score is a secondary summary, never the decision.
Step 5: the specificity lesson
Look at these two regulatory calls together:
regulatory <- gifts[gifts$gift_type == "regulatory", ]
regulatory[, c("gift_id", "complete", "minimum_missing_requirements")]
#> # A tibble: 4 × 3
#> gift_id complete minimum_missing_requirements
#> <chr> <lgl> <int>
#> 1 aspartate_chemoreception FALSE 1
#> 2 chemotaxis_signal_transduction TRUE 0
#> 3 phosphate_starvation_response FALSE 2
#> 4 serine_chemoreception FALSE 1The genome completes chemotaxis_signal_transduction but
not aspartate_chemoreception. What is missing from the
second?
aspartate <- gifts[gifts$gift_id == "aspartate_chemoreception", ]
unlist(aspartate$missing_requirements)
#> [1] "RF_ASPARTATE_RECEPTION"The genome carries a generic methyl-accepting chemoreceptor. That accession is assigned to dozens of receptors in a single genome, so it supports the claim “this organism has a chemosensory pathway” and cannot support “this organism senses aspartate”. A curated marker anchored on the characterised aspartate receptor would be needed for the second.
This is a rule, not an accident: the specificity of a claim never exceeds the specificity of its evidence. Widening the generic marker to make the specific GIFT fire would damage every other trait that marker also matches.
The same logic runs the other way between types:
gifts[gifts$gift_id == "flagellar_apparatus",
c("complete", "minimum_missing_requirements")]
#> # A tibble: 1 × 2
#> complete minimum_missing_requirements
#> <lgl> <int>
#> 1 FALSE 9This genome has a complete chemosensory pathway and no flagellum. Fusing the two into a “motility” trait would have hidden exactly that.
Step 6: trace a call back to genes
This is the step that makes a call reviewable.
trace_gift() walks the whole evidence chain down to your
own gene identifiers.
trace <- trace_gift(result, "xylose_degradation_isomerase")
trace[, c("reaction_id", "system_id", "component_id", "accession", "gene_id")]
#> # A tibble: 2 × 5
#> reaction_id system_id component_id accession gene_id
#> <chr> <chr> <chr> <chr> <chr>
#> 1 RHEA:22816 SYS_22816_XYLA COMP_22816_CATALYTIC K01805 g0330
#> 2 RHEA:10964 SYS_10964_XYLB COMP_10964_CATALYTIC K00854 g0331The chain for a metabolic GIFT runs route → reaction → enzyme system → component → marker → gene. For a machinery GIFT it runs implementation → function → system → component → marker → gene, under that type’s own names:
defense_trace <- trace_gift(result, "type_i_restriction_modification")
head(defense_trace[, c("function_id", "system_id", "accession", "gene_id")])
#> # A tibble: 3 × 4
#> function_id system_id accession gene_id
#> <chr> <chr> <chr> <chr>
#> 1 DF_TYPE_I_RM_SPECIFICITY SYS_DF_HSDS K01154 g0812
#> 2 DF_TYPE_I_RM_MODIFICATION SYS_DF_HSDM K03427 g0810
#> 3 DF_TYPE_I_RM_RESTRICTION SYS_DF_HSDR K01153 g0811If a reviewer disagrees with a call, this is where the disagreement becomes concrete: a specific marker accepted for a specific component.
Step 7: how the calls connect
Metabolic GIFTs compose through their declared boundary molecules — anchors — and nothing else. Internal reaction intermediates never create a connection.
edges <- gift_graph()
touching <- edges$from_gift %in% complete$gift_id |
edges$to_gift %in% complete$gift_id
edges[touching, c("from_gift", "shared_anchor", "to_gift", "edge_quality")]
#> # A tibble: 37 × 4
#> from_gift shared_anchor to_gift edge_quality
#> <chr> <chr> <chr> <chr>
#> 1 ammonium_assimilation GLUTAMATE glutamine_biosynthes… exact
#> 2 arabinoxylan_debranching XYLAN xylan_degradation exact
#> 3 arginine_deiminase_pathway AMMONIUM glutamine_biosynthes… exact
#> 4 aspartate_biosynthesis ASPARTATE asparagine_biosynthe… exact
#> 5 aspartate_biosynthesis ASPARTATE aspartate_semialdehy… exact
#> 6 aspartate_biosynthesis ASPARTATE pantothenate_biosynt… exact
#> 7 aspartate_biosynthesis ASPARTATE quinolinate_biosynth… exact
#> 8 citrate_fermentation PYRUVATE alanine_biosynthesis exact
#> 9 cysteine_degradation_sulfide AMMONIUM glutamine_biosynthes… exact
#> 10 cysteine_degradation_sulfide PYRUVATE alanine_biosynthesis exact
#> # ℹ 27 more rowsedge_quality is worth reading. exact means
both GIFTs declared the same anchor. compartment_inexact
means they declared the same molecule but one side left its compartment
unresolved, so the link assumes a location curation did not license.
Circular metabolism is derived from this graph rather than curated:
evaluate_gift_cycles(result)
#> # A tibble: 2 × 6
#> cycle_id named_cycle cycle_length supported status broken_at
#> <int> <chr> <int> <int> <chr> <chr>
#> 1 1 glyoxylate_cycle 4 0 absent acetyl_coa…
#> 2 2 citric_acid_cycle_oxidative 5 0 absent acetyl_coa…status is assembled from calls that were already made.
Whether a genome closes a cycle never changes whether any member GIFT is
complete.
What a positive call still does not mean
Worth restating now that you have some:
- it does not mean the machinery is expressed or the pathway carries flux;
- it does not mean the substrate is available in the organism’s environment;
-
type_i_restriction_modificationdoes not mean the genome is phage-resistant; -
chemotaxis_signal_transductiondoes not mean the organism is motile — this genome is the proof, since it has no flagellum.