在浩瀚基因数据库里寻找一个线索

ATGCTAGCTAGCGTAC

MALWMRLLPLLALLAL

ACGTTGACCTGA

TTACGTTGACCTGACC

Query:       ACGTTGACCTGA
Subject:   TTACGTTGACCTGACC
             ||||||||||||

ACGTTGACCTGA

ACG
CGT
GTT
TTG
TGA
GAC
ACC
CCT
CTG
TGA

{
    "ACG": [0],
    "CGT": [1],
    "GTT": [2],
    "TGA": [4, 9],
}

ACG

index.get("ACG")

TTACGTTGACCTGACC

ACG

ACG

Query:    ACG
Subject:  ACG

Query:    ACGT
Subject:  ACGT

Query:    ACGTT
Subject:  ACGTT

Query:    ACGTTGACCTGA
Subject:  ACGTTGACCTGA

Query:    ACGTTGACCTGA
Subject:  ACGTCGACATGA

ACG

Query:    A L V G
Subject:  A I V G

Query:    A C G T T G A C C T G A
Subject:  A C G T C G A C C T G A

Query:    ACGTTGACCTGA
          |||| |||||||
Subject:  ACGTCGACCTGA

ACGTTGACCTGA

AAAA

CCCC

AAAAAAAAAAAAAAAAAAAA

ATATATATATATATATATAT

AAAAAAAAAAAA

ATATATATATAT

AAAAAAAAAAAA

AAAAAA

AAA

AAA

AAAAAAAAAAAA

AACAACAACAAC

AAA

from collections import defaultdict


def build_word_index(query, word_size):
    """
    为查询序列建立短词 -> 出现位置列表。

    例如：
        ACGTT
        w = 3

    得到：
        ACG -> [0]
        CGT -> [1]
        GTT -> [2]
    """

    if not (1 <= word_size <= len(query)):
        raise ValueError("短词长度不合法")

    index = defaultdict(list)

    for i in range(len(query) - word_size + 1):
        word = query[i:i + word_size]
        index[word].append(i)

    return dict(index)


query = "ACGTTGACCTGA"

index = build_word_index(query, 3)

for word, positions in sorted(index.items()):
    print(word, positions)

def extend_one_side(
    query,
    subject,
    query_start,
    subject_start,
    word_size,
    direction,
    x_drop,
    match=2,
    mismatch=-1,
):
    """
    从已知的完全匹配种子向单侧延伸。

    direction:
        -1 -> 向左
         1 -> 向右

    返回：
        (最佳额外得分, 最佳额外延伸字符数)
    """

    current = 0
    best = 0

    steps = 0
    best_steps = 0

    while True:
        if direction == -1:
            qi = query_start - 1 - steps
            sj = subject_start - 1 - steps
        else:
            qi = query_start + word_size + steps
            sj = subject_start + word_size + steps

        if (
            qi < 0 or
            sj < 0 or
            qi >= len(query) or
            sj >= len(subject)
        ):
            break

        delta = (
            match
            if query[qi] == subject[sj]
            else mismatch
        )

        current += delta
        steps += 1

        if current > best:
            best = current
            best_steps = steps

        if best - current > x_drop:
            break

    return best, best_steps

best
best_steps

def mini_blast(
    query,
    database,
    word_size=3,
    x_drop=3,
    min_score=10,
):
    """
    教学版BLAST式DNA搜索。

    参数：
        query:
            查询DNA序列，例如ACGTTGACCTGA

        database:
            字典 {记录名: DNA序列}

        word_size:
            完全匹配种子长度

        x_drop:
            延伸分数允许下降的最大幅度

        min_score:
            最小报告得分

    只生成无Gap局部比对。
    不计算E-value。
    """

    query = query.upper()

    if not query or set(query) - set("ACGT"):
        raise ValueError("Query必须由ACGT组成")

    if not 1 <= word_size <= len(query):
        raise ValueError("不合法的word_size")

    if x_drop < 0:
        raise ValueError("x_drop不能为负")

    index = build_word_index(query, word_size)

    all_results = []

    for name, subject in database.items():
        subject = subject.upper()

        if set(subject) - set("ACGT"):
            raise ValueError(
                f"{name}: 只支持ACGT"
            )

        candidates = {}

        seed_count = 0

        for j in range(
            max(0, len(subject) - word_size + 1)
        ):
            word = subject[j:j + word_size]

            query_positions = index.get(word, [])

            for i in query_positions:
                seed_count += 1

                # 向左延伸
                left_score, left_steps = extend_one_side(
                    query,
                    subject,
                    i,
                    j,
                    word_size,
                    direction=-1,
                    x_drop=x_drop,
                )

                # 向右延伸
                right_score, right_steps = extend_one_side(
                    query,
                    subject,
                    i,
                    j,
                    word_size,
                    direction=1,
                    x_drop=x_drop,
                )

                q_start = i - left_steps
                q_end = i + word_size + right_steps

                s_start = j - left_steps
                s_end = j + word_size + right_steps

                score = (
                    2 * word_size
                    + left_score
                    + right_score
                )

                key = (
                    q_start, q_end,
                    s_start, s_end
                )

                candidate = {
                    "name": name,
                    "score": score,
                    "query_range": (q_start, q_end),
                    "subject_range": (s_start, s_end),
                    "query_segment": query[q_start:q_end],
                    "subject_segment": subject[s_start:s_end],
                }

                old = candidates.get(key)

                if old is None or score > old["score"]:
                    candidates[key] = candidate

        hits = sorted(
            (
                candidate
                for candidate in candidates.values()
                if candidate["score"] >= min_score
            ),
            key=lambda x: -x["score"],
        )

        all_results.append({
            "name": name,
            "seed_count": seed_count,
            "hits": hits,
            "best": hits[0] if hits else None,
        })

    all_results.sort(
        key=lambda r: (
            -r["best"]["score"]
            if r["best"] else float("inf")
        )
    )

    return all_results

DATABASE = {
    "sample_A_complete":
        "TTACGTTGACCTGACC",

    "sample_B_substitution":
        "GGACGTCGACCTGAT",

    "sample_C_insertion":
        "TTACGTTTGACCTGAA",

    "sample_D_fragment":
        "CCCCGACCTGAAAA",

    "sample_E_repeat":
        "AAAAAAAACCCCCCCC",

    "sample_F_other":
        "TCCAGATATCTGAC",
}


QUERY = "ACGTTGACCTGA"


results = mini_blast(
    query=QUERY,
    database=DATABASE,
    word_size=3,
    x_drop=3,
    min_score=10,
)


for record in results:
    print()
    print("数据库记录:", record["name"])
    print("种子命中数:", record["seed_count"])

    if record["best"] is None:
        print("没有达到阈值的候选")
        continue

    hit = record["best"]

    print("最佳分数:", hit["score"])
    print(
        "Query位置:",
        hit["query_range"],
    )
    print(
        "Subject位置:",
        hit["subject_range"],
    )

    print(hit["query_segment"])
    print(hit["subject_segment"])


# 全部12个碱基完全匹配：
# 每个匹配+2，所以最高分24。
assert results[0]["best"]["score"] == 24

print()
print("微型BLAST基础实验通过")

word_size=6

x_drop=0

word_size=2

def exact_ungapped_score(
    query,
    subject,
    match=2,
    mismatch=-1,
):
    """
    精确求解两条序列的最高无Gap局部匹配分数。

    沿所有对角线进行线性扫描。
    每个字符位置对最多处理一次。

    时间复杂度 O(m*n)。
    """

    m = len(query)
    n = len(subject)

    best_score = 0
    best_hit = None

    # diagonal = subject_position - query_position
    for diagonal in range(-m + 1, n):

        i = max(0, -diagonal)
        j = i + diagonal

        current = 0

        start_i = i
        start_j = j

        while i < m and j < n:

            delta = (
                match
                if query[i] == subject[j]
                else mismatch
            )

            current += delta

            if current <= 0:
                current = 0

                start_i = i + 1
                start_j = j + 1

            elif current > best_score:
                best_score = current

                best_hit = {
                    "query_range": (
                        start_i, i + 1
                    ),
                    "subject_range": (
                        start_j, j + 1
                    ),
                }

            i += 1
            j += 1

    return best_score, best_hit


query = "ACGTTGACCTGA"
subject = "TTACGTTGACCTGACC"

score, hit = exact_ungapped_score(
    query, subject
)

print("精确无Gap最高得分:", score)
print("最高分区间:", hit)

assert score == 24

Query:   AAAAAAAAAAAA
Subject: AACAACAACAAC

AAA

query = "AAAAAAAAAAAA"
subject = "AACAACAACAAC"

# 先执行教学版BLAST
search = mini_blast(
    query,
    {"test": subject},
    word_size=3,
    x_drop=3,
    min_score=1,
)

seed_hits = search[0]["seed_count"]

# 再执行精确无Gap局部比对
exact_score, exact_hit = exact_ungapped_score(
    query,
    subject,
)

print("种子命中数:", seed_hits)
print("精确最佳分数:", exact_score)

assert seed_hits == 0
assert exact_score > 0

print("成功观察到启发式搜索遗漏")

Query:    A C G T - T G A C
Subject:  A C G T C T G A C

>example_dna_sequence
ACGTTGACCTGA

ACGTTGACCTGA

ACG
CGT
GTT
TTG
GAC