phasing.py 文件源码

def generate_new_hsets(self, haplotype_set: HaplotypeSet,
                           possible_paths: List[Tuple[OrientedDNASegment]],
                           relevant_reads: RelevantReads
                           ) -> Iterable[HaplotypeSet]:
        """For a given haplotype set, generate all possible extensions at the
        current bubble. Only yield the new haplotype sets that have a relative
        likelihood above a given threshold."""

        if callable(self.threshold):
            threshold = self.threshold(len(relevant_reads))
        else:
            threshold = self.threshold

        if threshold == 0.0:
            threshold = float('-inf')
        else:
            threshold = math.log10(threshold)

        if self.start_of_block:
            # For the first bubble the order does not matter, as a permutation
            # in a different order will in the end result in the same haplotype
            # set.
            extension_iter = iter(combinations_with_replacement(possible_paths,
                                                                self.ploidy))
        else:
            # Otherwise all possible k-tuples of possible paths, because now
            # order does matter
            extension_iter = iter(product(possible_paths, repeat=self.ploidy))

        num_possible_sets = len(possible_paths)**self.ploidy

        for extension in extension_iter:
            ext_read_sets = []
            # Get graph reads of the extension
            for hap_ext in extension:
                # We index with [1:-1] to ignore the entrance and exit of the
                # bubble
                ext_read_sets.append(set(self.get_all_reads(hap_ext[1:-1])))

            rl = self.calculate_rl(haplotype_set, extension, ext_read_sets,
                                   relevant_reads, num_possible_sets)

            if rl >= threshold:
                new_set = haplotype_set.extend(extension, ext_read_sets)
                new_set.log_rl = rl
                yield new_set