注释版本基因对照表v3（终版）

介绍

1. NCC版本的基因组注释已经不需要了，而且最后的筛选似乎还残留了一些冗余的结果，从一致性来看并不是吾辈留下来的那部分，应该是脚本存在问题，因此需要更新一下筛选流程
2. NC和BA的比对，将两个版本注释得到的所有蛋白序列进行比对，筛选；存在问题的基因也纳入比较，表明，比对寻找BA版本的基因；同时对比较的结果，用脚本判断两个版本基因的位置是否重叠，进一步筛选。
3. 如果有后续需要补充的也放在这里
4. 表一共有三张，NC-NCC-BA（网站版）-BA（最终版）的对应表
5. NC-BA的blast比对结果
6. BA（网站版）-BA（最终版）的对应表
   1. 应为BA和NC的比对结果只有18953个结果，而BA的注释一共有36173个注释的基因
7. 后续如果有补充的话也都放在这篇文档中

流程

blastp比对和筛选（v2）

1. 先把两个版本的注释文件按照染色体分割成11个注释子文件
   1. 比对的时候，将对应的染色体上的蛋白进行比对。因为NC版本的基因编号完全看不出来是哪个染色体上的，跨染色体的比对可能会产生一些错误（比如chr1上的基因比对到chr11上的基因）
   2. 拆分两个版本的注释文件

awk '{if ($1 ~ /^Chr[0-9]+$/) print > $1"_nc.gff"}' Tgra.chr.gff
awk '{if ($1 ~ /^Chr[0-9]+$/) print > $1"_ba.gff3"}' female_ba_last.gff3

1. 提取两个注释文件的蛋白序列
   1. NC的

for gff_file in Chr{1..11}_nc.gff; do
    # 使用 gffread 从 GFF 文件中提取对应的 CDS 序列并翻译成蛋白质序列
    gffread "$gff_file" -g Tgra.chr.fa -y "${gff_file%.gff}.pep"
done

1. BA的

gffread -g female_chr_rename.fna Chr01_ba.gff3 -y Chr01_ba.pep
......

1. 创造11个文件夹，把需要比对的蛋白序列（nc和ba）丢到对应的文件夹当中
2. 建库

makeblastdb -in *ba.pep -dbtype prot -out ba_db

1. blastp比对
   1. 重新设定blastp比对的参数如下
      1. evalue设定为：1e-10（默认）
      2. -qcov_hsp_perc 70：设置查询序列的覆盖度为70%，意味着每个高得分区段（HSP）与查询序列的覆盖度必须大于或等于70%才能被返回

nohup blastp \
  -query *nc.pep \
  -db ba_db \
  -outfmt 6 \
  -out blast_results.tsv \
  -evalue 1e-10 \
  -num_threads 20 \
  -qcov_hsp_perc 70 \
  & 

1. 提取NC版注释的mrna的id
   1. 因为蛋白序列的id和mrna的id一样

awk '$3 == "mRNA" { 
    split($9, arr, ";"); 
    for (i in arr) {
        if (arr[i] ~ /^ID=/) {
            split(arr[i], id, "=");
            print id[2] "\t" $4 "\t" $5 "\t" $7
        }
    }
}' Chr1_nc.gff > Chr1_nc.id

1. 提取BA版本注释的mrna的id

awk '$3 == "mRNA" { 
    split($9, arr, ";"); 
    for (i in arr) {
        if (arr[i] ~ /^ID/) {
            # 提取ID并去除ID=前缀
            split(arr[i], id, "="); 
            print id[2] "\t" $4 "\t" $5 "\t" $7
        }
    }
}' Chr01_ba.gff3 > Chr01_ba.id

1. 筛选位置上存在重叠的基因

python gene_overlap.py -a Chr1_nc.id -b Chr01_ba.id -c blast_results.tsv -o filtered_blast_results.tsv

1. 脚本筛选
   1. 一致性大于70的
   2. 保留两个注释中基因id唯一对应的

python filter_blast_results2.py -i filtered_blast_results.tsv -o 01blast.tsv

制作不同注释版本基因id对应表

导入NC，sun，NCC，BA（last）的基因id，使用函数映射完成对应整合

附件

filter_blast_results2.py （v2）

1. 读取blast比对结果第1列内容（搜索序列的id），查找重复；找到id重复的内容时，对两个id所在行进行判定，判定条件如下
2. 判定最优比对结果的条件：identity、E_value、mismatches、gap_openings，这些条件依次进行，降序。
3. 完成第一轮筛选，开始第二轮筛选
4. 读取blast比对结果第2列内容（建库序列的id），查找重复；找到id重复的内容时，对两个id所在行进行判定，判定条件如上；（判定是否存在重复时，去掉后面-mRNA（数字）后判定）
5. 开始第三轮筛选
6. 读取第3列内容（一致性），筛选保留一致性>=70的结果
7. 完成筛选
8. 参数i，指定输入文件

import argparse
from collections import defaultdict
import re  # 用于正则表达式匹配

# 解析命令行参数
def parse_arguments():
    parser = argparse.ArgumentParser(description="Filter BLAST results based on best matches and identity consistency.")
    parser.add_argument("-i", "--input_file", required=True, help="Input BLAST result file.")
    parser.add_argument("-o", "--output_file", required=True, help="Output file for filtered results.")
    return parser.parse_args()

# 解析BLAST比对结果文件
def parse_blast_results(file_path):
    blast_results = []
    with open(file_path, "r") as f:
        for line in f:
            columns = line.strip().split("\t")
            blast_results.append({
                "query_id": columns[0],
                "subject_id": columns[1],
                "identity": float(columns[2]),
                "length": int(columns[3]),
                "mismatches": int(columns[4]),
                "gap_openings": int(columns[5]),
                "q_start": int(columns[6]),
                "q_end": int(columns[7]),
                "s_start": int(columns[8]),
                "s_end": int(columns[9]),
                "e_value": float(columns[10]),
                "bit_score": float(columns[11]),
                "consistency": float(columns[2])  # Assuming the 3rd column is consistency
            })
    return blast_results

# 通过正则表达式去掉 subject_id 中的 "-mRNA数字" 部分
def clean_subject_id(subject_id):
    return re.sub(r"-mRNA\d+", "", subject_id)

# 筛选最优比对
def select_best_matches_by_query(blast_results):
    best_matches = defaultdict(list)

    # 第一轮筛选：按查询序列ID（query_id）分组
    for result in blast_results:
        best_matches[result["query_id"]].append(result)

    filtered_results = []

    for query_id, results in best_matches.items():
        # 对每一组结果按Identity, E_value, Mismatches, Gap_openings进行排序
        sorted_results = sorted(
            results,
            key=lambda x: (-x["identity"], x["e_value"], x["mismatches"], x["gap_openings"])
        )
        # 选择最优结果
        filtered_results.append(sorted_results[0])

    return filtered_results

# 第二轮筛选：按比对序列ID（subject_id）分组，去除"-mRNA数字"部分
def select_best_matches_by_subject(blast_results):
    best_matches = defaultdict(list)

    # 第二轮筛选：按目标序列ID（subject_id）分组
    for result in blast_results:
        cleaned_subject_id = clean_subject_id(result["subject_id"])
        result["cleaned_subject_id"] = cleaned_subject_id  # 给结果添加清理后的subject_id
        best_matches[cleaned_subject_id].append(result)

    filtered_results = []

    for subject_id, results in best_matches.items():
        # 对每一组结果按Identity, E_value, Mismatches, Gap_openings进行排序
        sorted_results = sorted(
            results,
            key=lambda x: (-x["identity"], x["e_value"], x["mismatches"], x["gap_openings"])
        )
        # 选择最优结果
        filtered_results.append(sorted_results[0])

    return filtered_results

# 第三轮筛选：根据一致性（Identity）筛选
def filter_by_consistency(blast_results, min_consistency=70.0):
    return [result for result in blast_results if result["consistency"] >= min_consistency]

# 写入过滤后的结果
def write_filtered_results(filtered_results, output_file):
    with open(output_file, "w") as f:
        for result in filtered_results:
            f.write("\t".join([str(result[key]) for key in result if key != "cleaned_subject_id"]) + "\n")

# 主函数
def main():
    args = parse_arguments()

    # 解析BLAST结果文件
    blast_results = parse_blast_results(args.input_file)

    # 第一轮筛选：根据查询序列ID筛选最优比对
    best_by_query = select_best_matches_by_query(blast_results)

    # 第二轮筛选：根据目标序列ID筛选最优比对，并去掉-mRNA数字
    best_by_subject = select_best_matches_by_subject(best_by_query)

    # 第三轮筛选：根据一致性筛选
    final_filtered_results = filter_by_consistency(best_by_subject)

    # 写入筛选后的结果
    write_filtered_results(final_filtered_results, args.output_file)

if __name__ == "__main__":
    main()

gene_overlap.py

判定条件

1. 是否在同一条正负链上，如果不在同一条链上则丢掉
2. 比对的基因（gene1和gene2好了）：gene2起始<=gene1起始/终止<=gene2起始/终止，简单描述就是这么一回事：基因1任意区域是否落在基因2的区域内，或者两个基因完全重叠

import argparse

# 用于解析命令行参数
def parse_arguments():
    parser = argparse.ArgumentParser(description="Filter BLAST results based on gene overlap and strand.")
    parser.add_argument("-a", "--file_a", required=True, help="File containing gene1 info (ID, start, end, strand).")
    parser.add_argument("-b", "--file_b", required=True, help="File containing gene2 info (ID, start, end, strand).")
    parser.add_argument("-c", "--blast_results", required=True, help="File containing BLAST results.")
    parser.add_argument("-o", "--output", required=True, help="Output file for filtered results.")
    return parser.parse_args()

# 读取ID文件，返回ID对应的起始位置、终止位置和正负链信息
def load_gene_info(file_path):
    gene_info = {}
    with open(file_path, "r") as f:
        for line in f:
            parts = line.strip().split()
            gene_id = parts[0]
            start = int(parts[1])
            end = int(parts[2])
            strand = parts[3]
            gene_info[gene_id] = (start, end, strand)
    return gene_info

# 判断两个基因是否在同一链上
def same_strand(gene1, gene2):
    start1, end1, strand1 = gene1
    start2, end2, strand2 = gene2
    # 如果两个基因的链不同，直接返回False
    return strand1 == strand2

# 判断两个基因是否有重叠
def is_overlap(gene1, gene2):
    start1, end1, _ = gene1
    start2, end2, _ = gene2

    # 判断基因1的起始位置或终止位置是否在基因2的范围内
    if (start2 <= start1 <= end2) or (start2 <= end1 <= end2) or (start1 == start2 and end1 == end2):
        return True
    return False

# 处理BLAST结果并筛选符合条件的比对
def process_blast_results(blast_file, gene_info_a, gene_info_b, output_file):
    with open(blast_file, "r") as f_in, open(output_file, "w") as f_out:
        for line in f_in:
            columns = line.strip().split("\t")
            gene1_id = columns[0]
            gene2_id = columns[1]

            # 获取基因ID对应的信息
            if gene1_id in gene_info_a and gene2_id in gene_info_b:
                gene1 = gene_info_a[gene1_id]
                gene2 = gene_info_b[gene2_id]

                # 先判断基因是否在同一链上，不在同一链的直接跳过
                if not same_strand(gene1, gene2):
                    continue

                # 判断基因是否重叠
                if is_overlap(gene1, gene2):
                    f_out.write(line)

# 主函数
def main():
    args = parse_arguments()

    # 读取两个ID文件
    gene_info_a = load_gene_info(args.file_a)
    gene_info_b = load_gene_info(args.file_b)

    # 处理BLAST比对结果
    process_blast_results(args.blast_results, gene_info_a, gene_info_b, args.output)

if __name__ == "__main__":
    main()