Negative junction_length for reverse_complement sequence

[Zjianglin]

Hi developer, nice tool!

I have encountered a confused question about 'is_reverse_complement': True sequnce, whose junction_length attribute is negative. Is this a bug? (caculated from junction_end - junction_start). I think for later typically analysis procedures, suchsequene shoule be transformed to its 'reverse_complement' status to keep the 5'->3' direction. Howevere, the groudtruth of the coordinates for domains (such as v_region (v_sequence_start, v_sequence_end), 'junction_region'(junction_start, junction_end)) seems to be numbered according to raw sequence, how these coordinates shoule be properlly processed?

My procedure as below:

result = Experiment.on("human_tcrb").sequence( with_reverse_complement(0.3)).mutate(rate(0.02, 0.08))

for seq in result:
    ...:     if seq['productive'] and seq['is_reverse_complement']:
    ...:         print(seq)
    ...:         break
    ...: 
{'sequence': 'ccagcTctgagagccgggtcccggcgccgaagtactgaatgtttttggctAGcGttCcacagacatacacggccgagtccccctCctctgcaggctggatcttcagagtagagaatgatccatcaggcctctcGgctgagaatcActctCtgggcatacctgcGttatctataatagatcggctccggaagtaaGtcagAaattccaCtccctgcGcgGagCtatctctgtaccagaaaaggaaCttgtggccaaaaattggctcacatctcagagtcactgtttgtcccatttctgtcacctcatActtgggtgactggataatgccagcatc', 'germline_alignment': 'ggtcgtgactctcggcccagggccgcggcttcatgacttacaaaaaccgaNNgNaacgtgtctgtatgtgccggctcagggggacgagacgtccgacctagaagtctcatctcttactaggtagtccggagagtcgactcttagcgagacacccgtatggacgtaatagatattatctagccgaggccttcattgagtcgttaaggtcagggacgtgcttccatagagacatggtcttttcctttaacaccggtttttaaccgagtgtagagtctcagtgacaaacagggtaaagacagtggagtacgaacccactgacctattacggtcgtag', 'v_call': 'TRBV20/OR9-2*01', 'd_call': 'Short-D', 'j_call': 'TRBJ2-6*01,TRBJ2-7*01,TRBJ2-7*02', 'c_call': '', 'v_sequence_start': 54, 'v_sequence_end': 334, 'v_germline_start': 0, 'v_germline_end': 280, 'd_sequence_start': 52, 'd_sequence_end': 53, 'd_germline_start': 0, 'd_germline_end': 1, 'j_sequence_start': 0, 'j_sequence_end': 50, 'j_germline_start': 0, 'j_germline_end': 50, 'junction_start': 60, 'junction_end': 33, 'junction_nt': '', 'junction_aa': '', 'junction_length': -27, 'v_trim_5': 0, 'v_trim_3': 10, 'd_trim_5': 0, 'd_trim_3': 15, 'j_trim_5': 0, 'j_trim_3': 0, 'np1_region': 'G', 'np1_length': 1, 'np2_region': 'AG', 'np2_length': 2, 'mutation_rate': 0.045317, 'n_mutations': 15, 'mutations': '5:t>T,56:c>C,84:c>C,133:t>G,144:c>A,149:c>C,163:t>G,194:g>G,199:g>A,207:c>C,215:t>G,218:t>G,221:c>C,244:t>C,306:c>A', 'n_sequencing_errors': 0, 'sequencing_errors': '', 'n_pcr_errors': 0, 'pcr_errors': '', 'productive': True, 'stop_codon': False, 'vj_in_frame': True, 'note': '', 'is_reverse_complement': True, 'is_contaminant': False, 'sequence_length': 334}

Above sequence has a negative length: 'junction_start': 60, 'junction_end': 33, 'junction_nt': '', 'junction_aa': '', 'junction_length': -27. And this is a productive and vj_in_frame sequence, why its junction_aa and junction_nt is empty?

BTW, I want to use GenAIRR to simulate tens of human TCRB sequences. currently, I used following settings and procedures, is it correct?. Do you have any suggestion for my parameters settings and/or simulation procedures? Thanks.

LOCUS_TO_CONFIG = {
    'IGH': 'HUMAN_IGH_IMGT',
    'IGK': 'HUMAN_IGK_IMGT',
    'IGL': 'HUMAN_IGL_IMGT',
    'TRA': 'HUMAN_TCRA_IMGT',
    'TRB': 'HUMAN_TCRB_IMGT',
}
LOCUS_SIM_PARAMS = {
    'IGH': {'mutation': (0.01, 0.05), 'pcr_cycles': 30, 'five_prime_loss': 24, 'three_prime_loss': 18, 'indel_prob': 0.0025, 'ns_prob': 0.0025},
    'IGK': {'mutation': (0.002, 0.02), 'pcr_cycles': 25, 'five_prime_loss': 20, 'three_prime_loss': 14, 'indel_prob': 0.0015, 'ns_prob': 0.0020},
    'IGL': {'mutation': (0.002, 0.02), 'pcr_cycles': 25, 'five_prime_loss': 20, 'three_prime_loss': 14, 'indel_prob': 0.0015, 'ns_prob': 0.0020},
    'TRA': {'mutation': (0.0, 0.0005), 'pcr_cycles': 25, 'five_prime_loss': 18, 'three_prime_loss': 12, 'indel_prob': 0.0008, 'ns_prob': 0.0015},
    'TRB': {'mutation': (0.0, 0.0005), 'pcr_cycles': 25, 'five_prime_loss': 18, 'three_prime_loss': 12, 'indel_prob': 0.0008, 'ns_prob': 0.0015},
}
QUALITY_PROFILE = {'base': 0.001, 'peak': 0.02}

def _import_genairr():
    from GenAIRR import Experiment
    from GenAIRR.ops import rate, with_3prime_loss, with_5prime_loss, with_indels, with_ns, with_pcr, with_quality_profile, with_reverse_complement
    return {
        'Experiment': Experiment,
        'rate': rate,
        'with_3prime_loss': with_3prime_loss,
        'with_5prime_loss': with_5prime_loss,
        'with_indels': with_indels,
        'with_ns': with_ns,
        'with_pcr': with_pcr,
        'with_quality_profile': with_quality_profile,
        'with_reverse_complement': with_reverse_complement,
    }


def _experiment_for_locus(locus: str, productive: bool, seed: int):
    ga = _import_genairr()
    Experiment = ga['Experiment']
    cfg = LOCUS_TO_CONFIG[locus]
    exp = Experiment.on(cfg)
    params = LOCUS_SIM_PARAMS[locus]
    mutation_low, mutation_high = params['mutation']
    exp = exp.mutate(ga['rate'](mutation_low, mutation_high))
    exp = exp.prepare(ga['with_pcr'](error_rate=1e-4, cycles=params['pcr_cycles']))
    exp = exp.sequence(
        ga['with_5prime_loss'](min_remove=0, max_remove=params['five_prime_loss']),
        ga['with_3prime_loss'](min_remove=0, max_remove=params['three_prime_loss']),
        ga['with_quality_profile'](base=QUALITY_PROFILE['base'], peak=QUALITY_PROFILE['peak']),
        ga['with_reverse_complement'](prob=0.1),
    )
    exp = exp.observe(
        ga['with_indels'](prob=params['indel_prob']),
        ga['with_ns'](prob=params['ns_prob']),
    )
    sim = exp.compile(seed=seed, productive=productive)
    return sim

[MuteJester]

Hi @Zjianglin, thanks for the detailed report and apologies for the slow reply!

Short answer: the negative junction_length on reverse-complemented records was a real bug, and it's fixed in v2.2.0 (just released on PyPI). Your code is using the old v1.x DSL (.sequence(with_*), .observe(with_*), GenAIRR.ops), which the v2.0.0 release replaced with a flatter, more readable method chain. The migration is straightforward.

What changed and why your bug is fixed

In v2.2.0 the reverse-complement projection layer recomputes every coordinate (V/D/J windows, junction, NP regions) against the post-rev-comp sequence, so junction_end > junction_start always holds and junction/junction_aa are populated correctly. I verified locally with random_strand_orientation(prob=0.1) on HUMAN_TCRB_IMGT (100 records, 3 of them rev-comp): zero junction_length mismatches, zero empty junction_aa on productive records.

The AIRR-field renames you'll see going from your record to v2.2.0:

Old (v1.x)	New (v2.2.0)
is_reverse_complement	rev_comp
junction_nt	junction
n_sequencing_errors	n_quality_errors
np1_region / np2_region	np1 / np2
mutations (string)	counters: n_mutations, n_v_mutations, n_d_mutations, n_j_mutations, n_np_mutations

Your pipeline in the v2.2.0 DSL

The new DSL is one method chain (no .sequence/.observe/.prepare separation). Two important things up front:

• TCRs (TRA/TRB/TRG/TRD) don't undergo somatic hypermutation. T-cells lack the AID enzyme that drives SHM in B-cells. v2.2.0 enforces this: calling .mutate() on a TCR cartridge raises a ValueError telling you why. Your original mutation: (0.0, 0.0005) for TRB was effectively zero-rate SHM anyway, so the new behaviour matches your intent — just drop .mutate() entirely on TCR loci.
• rate= takes a single float, not a tuple. For variable mutation amounts per record use count=(min, max) (absolute counts).

Your TRB pipeline:

import GenAIRR as ga

result = (
    ga.Experiment.on("HUMAN_TCRB_IMGT")
      .allow_curatable_refdata()          # tolerate IMGT pseudogene anchor codons
      .recombine()
      .productive_only()
      # no .mutate() — TCRs don't have SHM (engine raises if you try)
      .pcr_amplify(rate=1e-4)             # per-record substitution rate
      .polymerase_indels(count=(0, 1))    # 0-1 indels per record
      .ambiguous_base_calls(count=(0, 1)) # 0-1 N's per record
      .end_loss_5prime(length=(0, 18))
      .end_loss_3prime(length=(0, 12))
      .random_strand_orientation(prob=0.1)
      .run_records(n=10_000, seed=42)
)
result.to_tsv("trb_sim.tsv")

I ran this locally on HUMAN_TCRB_IMGT (100 records, seed=42) and it produces the expected output: positive junction_length, populated junction_aa on productive records, correct rev_comp stamping, zero coordinate mismatches.

For your BCR loci (IGH/IGK/IGL), SHM is real, so .mutate() belongs:

result = (
    ga.Experiment.on("HUMAN_IGH_OGRDB")   # bundled name; HUMAN_IGH_IMGT also works
      .recombine()
      .productive_only()
      .mutate(model="s5f", rate=0.03)     # S5F-targeted SHM at 3% per base
      .pcr_amplify(rate=1e-4)
      .polymerase_indels(count=(0, 1))
      .ambiguous_base_calls(count=(0, 1))
      .end_loss_5prime(length=(0, 24))
      .end_loss_3prime(length=(0, 18))
      .random_strand_orientation(prob=0.1)
      .run_records(n=10_000, seed=42)
)

Parameter notes:

• pcr_amplify(rate=...) — the new API uses one per-record error rate, not a (rate, cycles) pair. For a typical 25–30 cycle amplification with 1e-4 per-base error, a rate of 1e-4 to 5e-4 is a reasonable starting point.
• polymerase_indels(count=(min, max)) and ambiguous_base_calls(count=(min, max)) — both take per-record (min, max) counts now rather than per-base probabilities. For ~350 bp reads, (0, 1) is roughly equivalent to the prob=0.0015 you had.
• end_loss_5prime(length=...) / end_loss_3prime(length=...) — replace with_5prime_loss / with_3prime_loss. length is (min, max) bases removed per record.
• random_strand_orientation(prob=...) — direct replacement for with_reverse_complement.
• .allow_curatable_refdata() on TCRB — the bundled IMGT TRB catalogue has a handful of alleles whose anchor codons don't strictly satisfy the productive-junction rules (pseudogenes / ORF alleles). This opt-in tells the engine to curate them at compile time instead of erroring. The OGRDB bundles (HUMAN_IGH_OGRDB) are pre-curated so this isn't needed there.

Where to look for the rest

• Quick start: https://mutejester.github.io/GenAIRR/getting-started/quick-start.html
• The Experiment builder (the full DSL): https://mutejester.github.io/GenAIRR/guides/experiment-builder.html
• Corruption + sequencing artefacts (PCR / indels / N / end-loss / strand): https://mutejester.github.io/GenAIRR/guides/corruption-sequencing.html
• AIRR record field catalogue: https://mutejester.github.io/GenAIRR/concepts/airr-record.html
• SHM targeting (for BCR): https://mutejester.github.io/GenAIRR/guides/shm-targeting.html

pip install -U GenAIRR gets you v2.2.0. Let me know if anything still misbehaves on your data happy to keep iterating.

Closing the issue once you've confirmed the fix; please reopen (or open a new one) if v2.2.0 still produces negative junction lengths.