comparison assess_alignment.py @ 10:fb905d0f8201 draft

planemo upload for repository https://github.com/pvanheus/polio_report commit aa90f911e6269aba792c9814c98a659e631400b2-dirty

9:acaaf49e2747

import argparse
import json
import sys
from dataclasses import dataclass
from typing import TextIO


@dataclass
class Sample:

                )
    min_start = min([int(al["leadingGaps"]) for al in msas])
    max_end = max([int(al["leadingGaps"]) + len(al["align"]) for al in msas])
    base_state = ["n"] * len(reference["align"])
    mismatch_bases = {}
    consensus = ''
    for i, reference_base in enumerate(reference["align"]):
        for k, al in enumerate(msas):
            leading_gaps = int(al["leadingGaps"])
            align_len = len(al["align"])
            if leading_gaps < i and (leading_gaps + align_len) > i:
                vp1pos = i - offset
                if vp1only and vp1pos < 0 or vp1pos > length:
                    # skip positions outside of vp1 gene region
                    continue
                al_base = al["align"][i - leading_gaps]
                consensus += al_base
                has_secondary_basecall = False
                if sec_is_conflict:
                    gappedTrace = data["gappedTraces"][k]
                    pos = i - int(gappedTrace["leadingGaps"])
                    # print(len(gappedTrace['basecallPos']), pos, k, len(gappedTrace['basecalls']), gappedTrace['basecallPos'][pos])

        if state == "M":
            # for mismatch store [pos_in_genome, pos_in_vp1, reference_base, sequenced_base]
            mismatch_list.append(
                [i + 1, i - offset + 1, reference["align"][i], mismatch_bases[i]]
            )
    return [conflicts, matches, mismatches, mismatch_list, consensus]


def analyse_trace_quality(json_file: TextIO) -> float:
    data = load_json(json_file)

        "perc_mismatches": percent_mismatches,
        "consensus": best_consensus
    }
    json.dump(info, open(args.output_filename, "w"))

    open(args.consensus_output_filename, "w").write(best_consensus)

10:fb905d0f8201

import argparse
import json
import sys
from dataclasses import dataclass
from textwrap import fill
from typing import TextIO


@dataclass
class Sample:

                )
    min_start = min([int(al["leadingGaps"]) for al in msas])
    max_end = max([int(al["leadingGaps"]) + len(al["align"]) for al in msas])
    base_state = ["n"] * len(reference["align"])
    mismatch_bases = {}
    base_start = 0
    for i, reference_base in enumerate(reference["align"]):
        if vp1only and not base_start:
            if data['gappedConsensus'][i] != '-':
                base_start = i
        for k, al in enumerate(msas):
            leading_gaps = int(al["leadingGaps"])
            align_len = len(al["align"])
            if leading_gaps < i and (leading_gaps + align_len) > i:
                vp1pos = i - offset
                if vp1only and vp1pos < 0 or vp1pos > length:
                    # skip positions outside of vp1 gene region
                    continue
                al_base = al["align"][i - leading_gaps]
                has_secondary_basecall = False
                if sec_is_conflict:
                    gappedTrace = data["gappedTraces"][k]
                    pos = i - int(gappedTrace["leadingGaps"])
                    # print(len(gappedTrace['basecallPos']), pos, k, len(gappedTrace['basecalls']), gappedTrace['basecallPos'][pos])

        if state == "M":
            # for mismatch store [pos_in_genome, pos_in_vp1, reference_base, sequenced_base]
            mismatch_list.append(
                [i + 1, i - offset + 1, reference["align"][i], mismatch_bases[i]]
            )
    if vp1only:
        # we have trim consensus bases before vp1 region
        if base_start <= offset:
            cons_start = offset - base_start
        else:
            cons_start = 0
        consensus = data['gapFreeConsensus'][cons_start:cons_start + length]
    else:
        consensus = data['gapFreeConsensus'] 
    return [conflicts, matches, mismatches, mismatch_list, consensus]


def analyse_trace_quality(json_file: TextIO) -> float:
    data = load_json(json_file)

        "perc_mismatches": percent_mismatches,
        "consensus": best_consensus
    }
    json.dump(info, open(args.output_filename, "w"))

    open(args.consensus_output_filename, "w").write(f'>{args.sample_name}\n' + fill(best_consensus))