def reads_to_variant_sequences(
        variant,
        reads,
        preferred_sequence_length,
        min_alt_rna_reads=MIN_ALT_RNA_READS,
        min_variant_sequence_coverage=MIN_VARIANT_SEQUENCE_COVERAGE,
        variant_sequence_assembly=VARIANT_SEQUENCE_ASSEMBLY):
    """
    Collapse variant-supporting RNA reads into consensus sequences of
    approximately the preferred length (may differ at the ends of transcripts),
    filter consensus sequences by length and number of supporting RNA reads.

    Parameters
    ----------
    variant : varcode.Variant

    reads : list of AlleleRead objects
        Should all support the same variant allele nucleotides.

    preferred_sequence_length : int
        Total number of nucleotides in the assembled sequences, including
        variant nucleotides.

    min_alt_rna_reads : int
        Drop sequences from loci which lack at least this number of reads
        that agree with the variant allele.

    min_variant_sequence_coverage : int
        Drop sequences which don't at least have this number of reads
        covering each cDNA position.

    variant_sequence_assembly : bool
        Construct variant sequences by merging overlapping reads. If False
        then variant sequences must be fully spanned by cDNA reads.

    Returns a collection of VariantSequence objects
    """
    # just in case variant_reads is a generator, convert it to a list
    variant_reads = list(filter_non_alt_reads_for_variant(variant, reads))

    if len(variant_reads) < min_alt_rna_reads:
        logger.info(
            "Skipping %s because only %d alt RNA reads (min=%d)",
            variant,
            len(variant_reads),
            min_alt_rna_reads)
        return []
    if len(variant_reads) == 0:
        return []

    alt_seq = get_single_allele_from_reads(variant_reads)

    # the number of context nucleotides on either side of the variant
    # is half the desired length (minus the number of variant nucleotides)
    n_alt_nucleotides = len(alt_seq)

    n_surrounding_nucleotides = preferred_sequence_length - n_alt_nucleotides
    max_nucleotides_after_variant = n_surrounding_nucleotides // 2

    # if the number of nucleotides we need isn't divisible by 2 then
    # prefer to have one more *before* the variant since we need the
    # prefix sequence to match against reference transcripts
    max_nucleotides_before_variant = (
        n_surrounding_nucleotides - max_nucleotides_after_variant)

    variant_sequences = initial_variant_sequences_from_reads(
        variant_reads=variant_reads,
        max_nucleotides_before_variant=max_nucleotides_before_variant,
        max_nucleotides_after_variant=max_nucleotides_after_variant)

    logger.info(
        "Initial pool of %d variant sequences (min length=%d, max length=%d)",
        len(variant_sequences),
        min(len(s) for s in variant_sequences),
        max(len(s) for s in variant_sequences))

    if variant_sequence_assembly:
        # this is a tricky parameter to set correctly:
        # by how many bases should two sequences overlap before
        # we merge, currently defaulting to either half the non-variant
        # nucleotides or 30 (whichever is smaller)
        variant_sequences = iterative_overlap_assembly(
            variant_sequences,
            min_overlap_size=min(
                MIN_VARIANT_SEQUENCE_ASSEMBLY_OVERLAP_SIZE,
                n_surrounding_nucleotides // 2))

    if variant_sequences:
        logger.info(
            "After overlap assembly: %d variant sequences (min length=%d, max length=%d)",
            len(variant_sequences),
            min(len(s) for s in variant_sequences),
            max(len(s) for s in variant_sequences))
    else:
        logger.info("After overlap assembly: 0 variant sequences")
        return []

    variant_sequences = filter_variant_sequences(
        variant_sequences=variant_sequences,
        preferred_sequence_length=preferred_sequence_length,
        min_variant_sequence_coverage=min_variant_sequence_coverage)

    if variant_sequences:
        logger.info(
            ("After coverage & length filtering: %d variant sequences "
             "(min length=%d, max length=%d)"),
            len(variant_sequences),
            min(len(s) for s in variant_sequences),
            max(len(s) for s in variant_sequences))
    else:
        logger.info("After coverage & length filtering: 0 variant sequences")
        return []

    # sort VariantSequence objects by decreasing order of supporting read
    # counts
    variant_sequences.sort(key=lambda vs: -len(vs.reads))
    return variant_sequences
def reads_to_variant_sequences ( variant , reads , preferred_sequence_length , min_alt_rna_reads = MIN_ALT_RNA_READS , min_variant_sequence_coverage = MIN_VARIANT_SEQUENCE_COVERAGE , variant_sequence_assembly = VARIANT_SEQUENCE_ASSEMBLY ) : variant_reads = list ( filter_non_alt_reads_for_variant ( variant , reads ) ) if len ( variant_reads ) < min_alt_rna_reads : logger . info ( "Skipping %s because only %d alt RNA reads (min=%d)" , variant , len ( variant_reads ) , min_alt_rna_reads ) return [ ] if len ( variant_reads ) == 0 : return [ ] alt_seq = get_single_allele_from_reads ( variant_reads ) n_alt_nucleotides = len ( alt_seq ) n_surrounding_nucleotides = preferred_sequence_length - n_alt_nucleotides max_nucleotides_after_variant = n_surrounding_nucleotides // 2 max_nucleotides_before_variant = ( n_surrounding_nucleotides - max_nucleotides_after_variant ) variant_sequences = initial_variant_sequences_from_reads ( variant_reads = variant_reads , max_nucleotides_before_variant = max_nucleotides_before_variant , max_nucleotides_after_variant = max_nucleotides_after_variant ) logger . info ( "Initial pool of %d variant sequences (min length=%d, max length=%d)" , len ( variant_sequences ) , min ( len ( s ) for s in variant_sequences ) , max ( len ( s ) for s in variant_sequences ) ) if variant_sequence_assembly : variant_sequences = iterative_overlap_assembly ( variant_sequences , min_overlap_size = min ( MIN_VARIANT_SEQUENCE_ASSEMBLY_OVERLAP_SIZE , n_surrounding_nucleotides // 2 ) ) if variant_sequences : logger . info ( "After overlap assembly: %d variant sequences (min length=%d, max length=%d)" , len ( variant_sequences ) , min ( len ( s ) for s in variant_sequences ) , max ( len ( s ) for s in variant_sequences ) ) else : logger . info ( "After overlap assembly: 0 variant sequences" ) return [ ] variant_sequences = filter_variant_sequences ( variant_sequences = variant_sequences , preferred_sequence_length = preferred_sequence_length , min_variant_sequence_coverage = min_variant_sequence_coverage ) if variant_sequences : logger . info ( ( "After coverage & length filtering: %d variant sequences " "(min length=%d, max length=%d)" ) , len ( variant_sequences ) , min ( len ( s ) for s in variant_sequences ) , max ( len ( s ) for s in variant_sequences ) ) else : logger . info ( "After coverage & length filtering: 0 variant sequences" ) return [ ] variant_sequences . sort ( key = lambda vs : - len ( vs . reads ) ) return variant_sequences
Collapse variant-supporting RNA reads into consensus sequences of
    approximately the preferred length (may differ at the ends of transcripts),
    filter consensus sequences by length and number of supporting RNA reads.

    Parameters
    ----------
    variant : varcode.Variant

    reads : list of AlleleRead objects
        Should all support the same variant allele nucleotides.

    preferred_sequence_length : int
        Total number of nucleotides in the assembled sequences, including
        variant nucleotides.

    min_alt_rna_reads : int
        Drop sequences from loci which lack at least this number of reads
        that agree with the variant allele.

    min_variant_sequence_coverage : int
        Drop sequences which don't at least have this number of reads
        covering each cDNA position.

    variant_sequence_assembly : bool
        Construct variant sequences by merging overlapping reads. If False
        then variant sequences must be fully spanned by cDNA reads.

    Returns a collection of VariantSequence objects