BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs

def parse_hmmer_output(self, filename, busco_query):

    processed_genes = []
    matched_genes = []
    records = defaultdict(list)

    data = []
    with open(filename, "r") as f:
        lines = f.readlines()

    for line in lines:
        if line.startswith("#"): # 跳过注释行
            continue
        else:
            values = line.split() # 把每一行拆成字段
            description = "".join(values[22:])
            row_data = tuple(values[:22]) + (description,)
            data.append(row_data)

    arr = np.array(
        data,
        dtype=[
            ("target_name", "U500"),
            ("target_accession", "U500"),
            ("tlen", "i4"),
            ("query_name", "U500"),
            ("query_accession", "U500"),
            ("qlen", "i4"),
            ("eval", "f8"),
            ("score", "f8"),
            ("bias", "f8"),
            ("dom_num", "i4"),
            ("ndom", "i4"),
            ("c-eval", "f8"),
            ("i-eval", "f8"),
            ("dom_score", "f8"),
            ("dom_bias", "f8"),
            ("hmm_coord_from", "i4"),
            ("hmm_coord_to", "i4"),
            ("ali_coord_from", "i4"),
            ("ali_coord_to", "i4"),
            ("env_coord_from", "i4"),
            ("env_coord_to", "i4"),
            ("acc", "f8"),
            ("description", "U500"),
        ],
    )

    target_hits_sorted, hmm_match_lengths = self.sort_hmmer_results(arr) # 排好序的基因ID列表以及每个基因的HMM匹配长度

    for target in target_hits_sorted:
        all_domain_hits = arr[arr["target_name"] == target]
        gene_id = self.get_gene_id(target)
        score = all_domain_hits[0]["score"]

        if not self.mode == "proteins" and not target in processed_genes:
            if self._check_overlap(matched_genes, gene_id):
                continue

        # Extract frame information (present in transcriptome mode)
        frame = (
            all_domain_hits[0]["description"]
            if "frame" in all_domain_hits[0]["description"]
            else None
        )
        # Store bitscore matches for each gene match. If match below cutoff, discard.
        if score < float(self.cutoff_dict[busco_query]["score"]):
            continue

        hmm_coords = [
            (hit["hmm_coord_from"], hit["hmm_coord_to"]) for hit in all_domain_hits
        ]
        env_coords = [
            (hit["env_coord_from"], hit["env_coord_to"]) for hit in all_domain_hits
        ]

        records[gene_id].append(
            {
                "hmm_matched_length": hmm_match_lengths[target],
                "hmm_profile_length": all_domain_hits[0]["qlen"],
                "hmm_coords": hmm_coords,
                "env_coords": env_coords,
                "score": score,
                "frame": frame,
                "orig gene ID": target,
                "ref gene ID": target,
            }
        )

        if gene_id not in matched_genes:
            matched_genes.append(gene_id)
            processed_genes.append(target)
    return records

def _sort_matches(self, matched_record, busco_query):
	# 元组的初始化...
	if self.config.get("busco_run", "datasets_version") == "odb10":
		# odb10数据集（使用统计学z-score方法）
	else:
	 # odb12及以后数据集（使用百分比方法）
	 if size >= 0.8 * profile_length:
      complete
   else:
      fragment
	return (
            busco_complete,
            busco_vlarge,
            busco_fragment,
            matched_genes_complete,
            matched_genes_vlarge,
            matched_genes_fragment,
        )

is_complete = {
  "B1": {
    "P1": [{"bitscore": 200, "length": 461}]
  },
  "B2": {
    "P3": [{"bitscore": 180, "length": 430}],
    "P4": [{"bitscore": 170, "length": 421}]
  }
}
is_very_large = {
  "B5": {
    "P9": [{"bitscore": 210, "length": 530}]
  }
}
is_fragment = {
  "B3": {
    "P5": [{"bitscore": 130, "length": 181}],
    "P6": [{"bitscore": 108, "length": 210}]
  }
}

is_complete = {
  "B1": {
    "P1": [{"bitscore": 200, "length": 461}]
  },
  "B2": {
    "P3": [{"bitscore": 180, "length": 430}],
    "P4": [{"bitscore": 170, "length": 421}]
  }
}
is_very_large = {
  "B5": {
    "P9": [{"bitscore": 210, "length": 530}]
  }
}
is_fragment = {
  "B3": {
    "P5": [{"bitscore": 130, "length": 181}]
  }
}

single_copy_buscos = {
  "B1": {
    "P1": [{"bitscore": 200, "length": 461}]
  },
  "B5": {
    "P9": [{"bitscore": 210, "length": 530}]
  }
}
multi_copy_buscos = {
  "B2": {
    "P3": [{"bitscore": 180, "length": 430}],
    "P4": [{"bitscore": 170, "length": 421}]
  }
}
fragmented_buscos = {
  "B3": {
    "P5": [{"bitscore": 130, "length": 181}]
  }
}

# 结果做的结构化
def record_results(self):
    self._get_busco_percentages()
    extra = ",E:{}%".format(self.e_percent) if self.e_percent else ""
    self.one_line_summary_raw = (
        "C:{}%[S:{}%,D:{}%],F:{}%,M:{}%,n:{}{}\t{}\n".format(
            self.complete_percent,
            self.s_percent,
            self.d_percent,
            self.f_percent,
            self.missing_percent,
            self.total_buscos,
            extra,
            "   ",
        )
    )
    self.one_line_summary = "Results:\t{}".format(self.one_line_summary_raw)

def _get_busco_percentages(self):
    self.single_copy = len(self.single_copy_buscos)  # 单拷贝 BUSCO 家族数
    self.multi_copy = len(self.multi_copy_buscos)  # 重复 BUSCO 家族数
    self.only_fragments = len(self.fragmented_buscos)  # 碎片化 BUSCO 家族数
    self.total_buscos = len(self.cutoff_dict) # 总 BUSCO 家族数
    self.num_missing = (
        self.total_buscos - self.single_copy - self.multi_copy - self.only_fragments
    ) # 缺失 BUSCO 家族数

    # Get percentage of each kind of BUSCO match
    self.s_percent = abs(round((self.single_copy / self.total_buscos) * 100, 1))
    self.d_percent = abs(round((self.multi_copy / self.total_buscos) * 100, 1))
    self.f_percent = abs(round((self.only_fragments / self.total_buscos) * 100, 1))
    self.complete_percent = abs(
        round(((self.single_copy + self.multi_copy) / self.total_buscos) * 100, 1)
    )
    self.missing_percent = abs(
        round((self.num_missing / self.total_buscos) * 100, 1)
    )

    if self.miniprot_pipeline:
        self._get_stop_codon_percent_and_avg_identity()
    else:
        self.e_percent = 0
        self.avg_identity = None
        self.complete_stop_codon_count = 0

    return