RIGHOR

本包基于 IGoR，旨在学习 V(D)J 重组的模型。

它可以

• 生成序列
• 评估序列（推断最可能的重组场景）
• 计算 "pgen"

可能更容易使用配套的 Python 包（pip install righor），但直接在 Rust 中工作也应可行。

如何使用 Python 包

加载模型

import righor
import matplotlib.pyplot as plt
import seaborn
import pandas as pd
from tqdm.notebook import tqdm
from collections import Counter
import numpy as np


igor_model = righor.load_model("human", "trb")

# alternatively, you can load a model from igor files
# igor_model = righor.load_model_from_files(params.txt, marginals.txt, anchor_v.csv, anchor_j.csv)

快速生成序列

# Create a generator object
generator = igor_model.generator(seed=42) # or igor_model.generator() to run it without a seed

# Generate 10'000 functional sequences (not out-of-frame, no stop codons, right boundaries)
for _ in tqdm(range(10000)):
    # generate_without_errors ignore Igor error model, use "generate" if this is needed
    sequence = generator.generate_without_errors(functional=True)
    if "IGH" in sequence.cdr3_aa:
        print("TRB CDR3 containing \"IGH\":", sequence.cdr3_aa)

# Generate one sequence with a particular V/J genes family
V_genes = righor.genes_matching("TRBV5", igor_model) # return all the V genes that match TRBV5
J_genes = righor.genes_matching("TRBJ", igor_model) # all the J genes
generator = igor_model.generator(seed=42, available_v=V_genes, available_j=J_genes)
generation_result = generator.generate_without_errors(functional=True)
print("Result:")
print(generation_result)
print("Explicit recombination event:")
print(generation_result.recombination_event)

评估给定的序列

my_sequence = "ACCCTCCAGTCTGCCAGGCCCTCACATACCTCTCAGTACCTCTGTGCCAGCAGTGAGGACAGGGACGTCACTGAAGCTTTCTTTGGACAAGGCACC"

# first align the sequence
align_params = righor.AlignmentParameters() # default alignment parameters
aligned_sequence = igor_model.align_sequence(my_sequence, align_params)

# we can also align a sequence from a CDR3 and a list of V-genes and J-genes (much faster)
# v_genes = righor.genes_matching("TRBV1", igor_model)
# j_genes = righor.genes_matching("TRBJ1", igor_model)
# igor_model.align_cdr3('TGTGTGAGAGATATTGTAGTAGTACCAGCTGCTAACCGCTTTCCTTCTTACTACTACTACTACTACATGGACGTCTGG', v_genes, j_genes)

# then evaluate it
infer_params = righor.InferenceParameters() # default inference parameters
result_inference = igor_model.evaluate(aligned_sequence, infer_params)

# Most likely scenario
best_event = result_inference.best_event

print(f"Probability that this specific event chain created the sequence: {best_event.likelihood / result_inference.likelihood:.2f}.")
print(f"Reconstructed sequence (without errors):", best_event.reconstructed_sequence)
print(f"Pgen: {result_inference.pgen:.1e}")

推断模型

# here we just generate the sequences needed
generator = igor_model.generator()
example_seq = generator.generate(False)
sequences = [generator.generate(False).full_seq for _ in range(1000)]

# define parameters for the alignment and the inference
align_params = righor.AlignmentParameters()
align_params.left_v_cutoff = 40
infer_params = righor.InferenceParameters()

# generate an uniform model as a starting point
# (it's generally *much* faster to start from an already inferred model)
model = igor_model.copy()
model.p_ins_vd = np.ones(model.p_ins_vd.shape)
model.error_rate = 0

# align multiple sequences at once
aligned_sequences = model.align_all_sequences(sequences, align_params)

# multiple round of expectation-maximization to infer the model
models = {}
model = igor_model.uniform()
model.error_rate = 0
models[0] = model
for ii in tqdm(range(35)):
    models[ii+1] = models[ii].copy()
    models[ii+1].infer(aligned_sequences, infer_params)

额外功能

与 IGoR 的主要区别

• "动态规划" 方法，而不是对所有事件求和，我们首先对事件的总和进行预计算。这意味着我们可以运行它，使用未定义的核苷酸如 N（至少在理论上，我需要添加对这些的完全支持）。
• D 基因对齐不太受限制

限制

• 在运行之前需要去除 V 基因侧的任何引物/末端
• 读取的长度需要足够长，才能完全覆盖 CDR3（即使它特别长）

编程相关

• 有一个用于网页使用的 wasm 版本。
• Python 版本现在在另一个 crate 中。
• 要永久添加模型，请将其添加到 "models.json"。同一类别中的第一个模型是默认模型。每个字段都是一个独立模型。chain 和 species 中的元素始终应该是小写。

当前状态

• 速度还不错（大约 50 个序列/秒）。可能会稍微快一点。我认为某些范围应该被迭代器替换。
• pgen 可以工作，但由于我考虑的 D 基因对齐比 Quentin 多，当 endD 和 startD 真的非常接近时，可能会有一些问题。