# HG changeset patch # User mheinzl # Date 1601831979 0 # Node ID e5953c54cfb53194c8024258eb107c63cafefa58 planemo upload for repository https://github.com/gpovysil/VariantAnalyzerGalaxy/tree/master/tools/variant_analyzer commit ee4a8e6cf290e6c8a4d55f9cd2839d60ab3b11c8 diff -r 000000000000 -r e5953c54cfb5 mut2read.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mut2read.py Sun Oct 04 17:19:39 2020 +0000 @@ -0,0 +1,155 @@ +#!/usr/bin/env python + +"""mut2read.py + +Author -- Gundula Povysil +Contact -- povysil@bioinf.jku.at + +Takes a tabular file with mutations and a BAM file as input and prints +all tags of reads that carry the mutation to a user specified output file. +Creates fastq file of reads of tags with mutation. + +======= ========== ================= ================================ +Version Date Author Description +0.2.1 2019-10-27 Gundula Povysil - +======= ========== ================= ================================ + +USAGE: python mut2read.py DCS_Mutations.tabular DCS.bam Aligned_Families.tabular Interesting_Reads.fastq + tag_count_dict.json +""" + +import argparse +import json +import os +import sys + +import numpy as np +import pysam + + +def make_argparser(): + parser = argparse.ArgumentParser(description='Takes a tabular file with mutations and a BAM file as input and prints all tags of reads that carry the mutation to a user specified output file and creates a fastq file of reads of tags with mutation.') + parser.add_argument('--mutFile', + help='TABULAR file with DCS mutations.') + parser.add_argument('--bamFile', + help='BAM file with aligned DCS reads.') + parser.add_argument('--familiesFile', + help='TABULAR file with aligned families.') + parser.add_argument('--outputFastq', + help='Output FASTQ file of reads with mutations.') + parser.add_argument('--outputJson', + help='Output JSON file to store collected data.') + return parser + + +def mut2read(argv): + parser = make_argparser() + args = parser.parse_args(argv[1:]) + + file1 = args.mutFile + file2 = args.bamFile + file3 = args.familiesFile + outfile = args.outputFastq + json_file = args.outputJson + + if os.path.isfile(file1) is False: + sys.exit("Error: Could not find '{}'".format(file1)) + + if os.path.isfile(file2) is False: + sys.exit("Error: Could not find '{}'".format(file2)) + + if os.path.isfile(file3) is False: + sys.exit("Error: Could not find '{}'".format(file3)) + + # read mut file + with open(file1, 'r') as mut: + mut_array = np.genfromtxt(mut, skip_header=1, delimiter='\t', comments='#', dtype=str) + + # read dcs bam file + # pysam.index(file2) + bam = pysam.AlignmentFile(file2, "rb") + + # get tags + tag_dict = {} + cvrg_dict = {} + + if mut_array.shape == (1,13): + mut_array = mut_array.reshape((1, len(mut_array))) + + for m in range(len(mut_array[:, 0])): + print(str(m + 1) + " of " + str(len(mut_array[:, 0]))) + chrom = mut_array[m, 1] + stop_pos = mut_array[m, 2].astype(int) + chrom_stop_pos = str(chrom) + "#" + str(stop_pos) + ref = mut_array[m, 9] + alt = mut_array[m, 10] + + dcs_len = [] + + for pileupcolumn in bam.pileup(chrom.tostring(), stop_pos - 2, stop_pos, max_depth=100000000): + + if pileupcolumn.reference_pos == stop_pos - 1: + count_alt = 0 + count_ref = 0 + count_indel = 0 + count_n = 0 + count_other = 0 + count_lowq = 0 + print("unfiltered reads=", pileupcolumn.n, "filtered reads=", len(pileupcolumn.pileups), + "difference= ", len(pileupcolumn.pileups) - pileupcolumn.n) + for pileupread in pileupcolumn.pileups: + if not pileupread.is_del and not pileupread.is_refskip: + # query position is None if is_del or is_refskip is set. + nuc = pileupread.alignment.query_sequence[pileupread.query_position] + dcs_len.append(len(pileupread.alignment.query_sequence)) + if nuc == alt: + count_alt += 1 + tag = pileupread.alignment.query_name + if tag in tag_dict: + tag_dict[tag][chrom_stop_pos] = alt + else: + tag_dict[tag] = {} + tag_dict[tag][chrom_stop_pos] = alt + elif nuc == ref: + count_ref += 1 + elif nuc == "N": + count_n += 1 + elif nuc == "lowQ": + count_lowq += 1 + else: + count_other += 1 + else: + count_indel += 1 + dcs_median = np.median(np.array(dcs_len)) + cvrg_dict[chrom_stop_pos] = (count_ref, count_alt, dcs_median) + + print("coverage at pos %s = %s, ref = %s, alt = %s, other bases = %s, N = %s, indel = %s, low quality = %s, median length of DCS = %s\n" % + (pileupcolumn.pos, count_ref + count_alt, count_ref, count_alt, count_other, count_n, + count_indel, count_lowq, dcs_median)) + bam.close() + + with open(json_file, "w") as f: + json.dump((tag_dict, cvrg_dict), f) + + # create fastq from aligned reads + with open(outfile, 'w') as out: + with open(file3, 'r') as families: + for line in families: + line = line.rstrip('\n') + splits = line.split('\t') + tag = splits[0] + + if tag in tag_dict: + str1 = splits[4] + curr_seq = str1.replace("-", "") + str2 = splits[5] + curr_qual = str2.replace(" ", "") + + out.write("@" + splits[0] + "." + splits[1] + "." + splits[2] + "\n") + out.write(curr_seq + "\n") + out.write("+" + "\n") + out.write(curr_qual + "\n") + + +if __name__ == '__main__': + sys.exit(mut2read(sys.argv)) diff -r 000000000000 -r e5953c54cfb5 mut2read.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mut2read.xml Sun Oct 04 17:19:39 2020 +0000 @@ -0,0 +1,65 @@ + + + Extracts all tags that carry a mutation in the duplex consensus sequence (DCS) + + va_macros.xml + + + + + + + + + + + + + + + + + + + + + + + `_. + +**Dataset 3:** Tabular file with reads as produced by the +**Du Novo: Align families** tool of the `Du Novo Analysis Pipeline +`_ + +**Output** + +The output is a json file containing dictonaries of the tags of reads containing mutations +in the DCS and a fastq file of all reads of these tags. + + ]]> + + + diff -r 000000000000 -r e5953c54cfb5 mut2sscs.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mut2sscs.py Sun Oct 04 17:19:39 2020 +0000 @@ -0,0 +1,133 @@ +#!/usr/bin/env python + +"""mut2sscs.py + +Author -- Gundula Povysil +Contact -- povysil@bioinf.jku.at + +Takes a tabular file with mutations from DCS and a BAM file of SSCS as input +and extracts all tags of reads that carry the mutation. +Calculates statistics about number of ab/ba/duplex per mutation. + +======= ========== ================= ================================ +Version Date Author Description +0.2.1 2019-10-27 Gundula Povysil - +======= ========== ================= ================================ + +USAGE: python mut2sscs.py DCS_Mutations.tabular SSCS.bam SSCS_counts.json + +""" + +from __future__ import division + +import argparse +import json +import os +import sys + +import numpy as np +import pysam + + +def make_argparser(): + parser = argparse.ArgumentParser(description='Takes a tabular file with mutations and a BAM file as input and prints all tags of reads that carry the mutation to a user specified output file.') + parser.add_argument('--mutFile', + help='TABULAR file with DCS mutations.') + parser.add_argument('--bamFile', + help='BAM file with aligned SSCS reads.') + parser.add_argument('--outputJson', + help='Output JSON file to store SSCS counts.') + return parser + + +def mut2sscs(argv): + parser = make_argparser() + args = parser.parse_args(argv[1:]) + + file1 = args.mutFile + file2 = args.bamFile + sscs_counts_json = args.outputJson + + if os.path.isfile(file1) is False: + sys.exit("Error: Could not find '{}'".format(file1)) + + if os.path.isfile(file2) is False: + sys.exit("Error: Could not find '{}'".format(file2)) + + # 1. read mut file + with open(file1, 'r') as mut: + mut_array = np.genfromtxt(mut, skip_header=1, delimiter='\t', comments='#', dtype=str) + + # 2 read SSCS bam file + # pysam.index(file2) + bam = pysam.AlignmentFile(file2, "rb") + + # get tags + mut_pos_dict = {} + ref_pos_dict = {} + if mut_array.shape == (1,13): + mut_array = mut_array.reshape((1, len(mut_array))) + + for m in range(0, len(mut_array[:, 0])): + print(str(m + 1) + " of " + str(len(mut_array[:, 0]))) + chrom = mut_array[m, 1] + stop_pos = mut_array[m, 2].astype(int) + chrom_stop_pos = str(chrom) + "#" + str(stop_pos) + ref = mut_array[m, 9] + alt = mut_array[m, 10] + + for pileupcolumn in bam.pileup(chrom.tostring(), stop_pos - 2, stop_pos, max_depth=1000000000): + if pileupcolumn.reference_pos == stop_pos - 1: + count_alt = 0 + count_ref = 0 + count_indel = 0 + print("unfiltered reads=", pileupcolumn.n, "filtered reads=", len(pileupcolumn.pileups), + "difference= ", len(pileupcolumn.pileups) - pileupcolumn.n) + for pileupread in pileupcolumn.pileups: + if not pileupread.is_del and not pileupread.is_refskip: + tag = pileupread.alignment.query_name + abba = tag[-2:] + # query position is None if is_del or is_refskip is set. + if pileupread.alignment.query_sequence[pileupread.query_position] == alt: + count_alt += 1 + if chrom_stop_pos in mut_pos_dict: + if abba in mut_pos_dict[chrom_stop_pos]: + mut_pos_dict[chrom_stop_pos][abba] += 1 + else: + mut_pos_dict[chrom_stop_pos][abba] = 1 + else: + mut_pos_dict[chrom_stop_pos] = {} + mut_pos_dict[chrom_stop_pos][abba] = 1 + elif pileupread.alignment.query_sequence[pileupread.query_position] == ref: + count_ref += 1 + if chrom_stop_pos in ref_pos_dict: + if abba in ref_pos_dict[chrom_stop_pos]: + ref_pos_dict[chrom_stop_pos][abba] += 1 + else: + ref_pos_dict[chrom_stop_pos][abba] = 1 + else: + ref_pos_dict[chrom_stop_pos] = {} + ref_pos_dict[chrom_stop_pos][abba] = 1 + else: + count_indel += 1 + + print("coverage at pos %s = %s, ref = %s, alt = %s, indel = %s,\n" % + (pileupcolumn.pos, count_ref + count_alt, count_ref, count_alt, count_indel)) + + # if mutation is in DCS file but not in SSCS, then set counts to NA + if chrom_stop_pos not in mut_pos_dict.keys(): + mut_pos_dict[chrom_stop_pos] = {} + mut_pos_dict[chrom_stop_pos]["ab"] = 0 + mut_pos_dict[chrom_stop_pos]["ba"] = 0 + ref_pos_dict[chrom_stop_pos] = {} + ref_pos_dict[chrom_stop_pos]["ab"] = 0 + ref_pos_dict[chrom_stop_pos]["ba"] = 0 + bam.close() + + # save counts + with open(sscs_counts_json, "w") as f: + json.dump((mut_pos_dict, ref_pos_dict), f) + + +if __name__ == '__main__': + sys.exit(mut2sscs(sys.argv)) diff -r 000000000000 -r e5953c54cfb5 mut2sscs.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mut2sscs.xml Sun Oct 04 17:19:39 2020 +0000 @@ -0,0 +1,59 @@ + + + Extracts all tags from the single stranded consensus sequence (SSCS) bam file that carry a mutation at the same position a mutation is called in the duplex consensus sequence (DCS) and calculates their frequencies + + va_macros.xml + + + + + + + + + + + + + + + + + + + `_. + +**Dataset 3:** Tabular file with reads as produced by the +**Du Novo: Align families** tool of the `Du Novo Analysis Pipeline +`_ + +**Output** + +The output is a json file containing dictonaries with stats of tags that carry a mutation in the SSCS +at the same position a mutation is called in the DCS. + + ]]> + + + diff -r 000000000000 -r e5953c54cfb5 read2mut.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/read2mut.py Sun Oct 04 17:19:39 2020 +0000 @@ -0,0 +1,953 @@ +#!/usr/bin/env python + +"""read2mut.py + +Author -- Gundula Povysil +Contact -- povysil@bioinf.jku.at + +Looks for reads with mutation at known +positions and calculates frequencies and stats. + +======= ========== ================= ================================ +Version Date Author Description +0.2.1 2019-10-27 Gundula Povysil - +======= ========== ================= ================================ + + +USAGE: python read2mut.py --mutFile DCS_Mutations.tabular --bamFile Interesting_Reads.trim.bam + --inputJson tag_count_dict.json --sscsJson SSCS_counts.json + --outputFile mutant_reads_summary_short_trim.xlsx --thresh 10 --phred 20 --trim 10 --chimera_correction + +""" + +from __future__ import division + +import argparse +import json +import operator +import os +import re +import sys + +import numpy as np +import pysam +import xlsxwriter + + +def make_argparser(): + parser = argparse.ArgumentParser(description='Takes a tabular file with mutations, a BAM file and JSON files as input and prints stats about variants to a user specified output file.') + parser.add_argument('--mutFile', + help='TABULAR file with DCS mutations.') + parser.add_argument('--bamFile', + help='BAM file with aligned raw reads of selected tags (FASTQ created by mut2read.py - trimming with Trimmomatic - alignment with bwa).') + parser.add_argument('--inputJson', + help='JSON file with data collected by mut2read.py.') + parser.add_argument('--sscsJson', + help='JSON file with SSCS counts collected by mut2sscs.py.') + parser.add_argument('--outputFile', + help='Output xlsx file of mutation details.') + parser.add_argument('--thresh', type=int, default=0, + help='Integer threshold for displaying mutations. Only mutations occuring less than thresh times are displayed. Default of 0 displays all.') + parser.add_argument('--phred', type=int, default=20, + help='Integer threshold for Phred score. Only reads higher than this threshold are considered. Default 20.') + parser.add_argument('--trim', type=int, default=10, + help='Integer threshold for assigning mutations at start and end of reads to lower tier. Default 10.') + parser.add_argument('--chimera_correction', action="store_true", + help='Add additional tier for chimeric variants and correct the variant frequencies') + + return parser + + +def safe_div(x, y): + if y == 0: + return None + return x / y + + +def read2mut(argv): + parser = make_argparser() + args = parser.parse_args(argv[1:]) + file1 = args.mutFile + file2 = args.bamFile + json_file = args.inputJson + sscs_json = args.sscsJson + outfile = args.outputFile + thresh = args.thresh + phred_score = args.phred + trim = args.trim + chimera_correction = args.chimera_correction + + if os.path.isfile(file1) is False: + sys.exit("Error: Could not find '{}'".format(file1)) + if os.path.isfile(file2) is False: + sys.exit("Error: Could not find '{}'".format(file2)) + if os.path.isfile(json_file) is False: + sys.exit("Error: Could not find '{}'".format(json_file)) + if thresh < 0: + sys.exit("Error: thresh is '{}', but only non-negative integers allowed".format(thresh)) + if phred_score < 0: + sys.exit("Error: phred is '{}', but only non-negative integers allowed".format(phred_score)) + if trim < 0: + sys.exit("Error: trim is '{}', but only non-negative integers allowed".format(thresh)) + + # 1. read mut file + with open(file1, 'r') as mut: + mut_array = np.genfromtxt(mut, skip_header=1, delimiter='\t', comments='#', dtype=str) + + # 2. load dicts + with open(json_file, "r") as f: + (tag_dict, cvrg_dict) = json.load(f) + + with open(sscs_json, "r") as f: + (mut_pos_dict, ref_pos_dict) = json.load(f) + + # 3. read bam file + # pysam.index(file2) + bam = pysam.AlignmentFile(file2, "rb") + + # 4. create mut_dict + mut_dict = {} + mut_read_pos_dict = {} + mut_read_dict = {} + reads_dict = {} + if mut_array.shape == (1, 13): + mut_array = mut_array.reshape((1, len(mut_array))) + + for m in range(0, len(mut_array[:, 0])): + print(str(m + 1) + " of " + str(len(mut_array[:, 0]))) + # for m in range(0, 5): + chrom = mut_array[m, 1] + stop_pos = mut_array[m, 2].astype(int) + chrom_stop_pos = str(chrom) + "#" + str(stop_pos) + ref = mut_array[m, 9] + alt = mut_array[m, 10] + mut_dict[chrom_stop_pos] = {} + mut_read_pos_dict[chrom_stop_pos] = {} + reads_dict[chrom_stop_pos] = {} + + for pileupcolumn in bam.pileup(chrom.tostring(), stop_pos - 2, stop_pos, max_depth=1000000000): + if pileupcolumn.reference_pos == stop_pos - 1: + count_alt = 0 + count_ref = 0 + count_indel = 0 + count_n = 0 + count_other = 0 + count_lowq = 0 + n = 0 + print("unfiltered reads=", pileupcolumn.n, "filtered reads=", len(pileupcolumn.pileups), + "difference= ", len(pileupcolumn.pileups) - pileupcolumn.n) + for pileupread in pileupcolumn.pileups: + n += 1 + if not pileupread.is_del and not pileupread.is_refskip: + tag = pileupread.alignment.query_name + nuc = pileupread.alignment.query_sequence[pileupread.query_position] + phred = ord(pileupread.alignment.qual[pileupread.query_position]) - 33 + if phred < phred_score: + nuc = "lowQ" + if tag not in mut_dict[chrom_stop_pos]: + mut_dict[chrom_stop_pos][tag] = {} + if nuc in mut_dict[chrom_stop_pos][tag]: + mut_dict[chrom_stop_pos][tag][nuc] += 1 + else: + mut_dict[chrom_stop_pos][tag][nuc] = 1 + if tag not in mut_read_pos_dict[chrom_stop_pos]: + mut_read_pos_dict[chrom_stop_pos][tag] = np.array(pileupread.query_position) + 1 + reads_dict[chrom_stop_pos][tag] = len(pileupread.alignment.query_sequence) + else: + mut_read_pos_dict[chrom_stop_pos][tag] = np.append( + mut_read_pos_dict[chrom_stop_pos][tag], pileupread.query_position + 1) + reads_dict[chrom_stop_pos][tag] = np.append( + reads_dict[chrom_stop_pos][tag], len(pileupread.alignment.query_sequence)) + + if nuc == alt: + count_alt += 1 + if tag not in mut_read_dict: + mut_read_dict[tag] = {} + mut_read_dict[tag][chrom_stop_pos] = alt + else: + mut_read_dict[tag][chrom_stop_pos] = alt + elif nuc == ref: + count_ref += 1 + elif nuc == "N": + count_n += 1 + elif nuc == "lowQ": + count_lowq += 1 + else: + count_other += 1 + else: + count_indel += 1 + + print("coverage at pos %s = %s, ref = %s, alt = %s, other bases = %s, N = %s, indel = %s, low quality = %s\n" % (pileupcolumn.pos, count_ref + count_alt, count_ref, count_alt, count_other, count_n, count_indel, count_lowq)) + + for read in bam.fetch(until_eof=True): + if read.is_unmapped: + pure_tag = read.query_name[:-5] + nuc = "na" + for key in tag_dict[pure_tag].keys(): + if key not in mut_dict: + mut_dict[key] = {} + if read.query_name not in mut_dict[key]: + mut_dict[key][read.query_name] = {} + if nuc in mut_dict[key][read.query_name]: + mut_dict[key][read.query_name][nuc] += 1 + else: + mut_dict[key][read.query_name][nuc] = 1 + bam.close() + + # 5. create pure_tags_dict + pure_tags_dict = {} + for key1, value1 in sorted(mut_dict.items()): + i = np.where(np.array(['#'.join(str(i) for i in z) + for z in zip(mut_array[:, 1], mut_array[:, 2])]) == key1)[0][0] + ref = mut_array[i, 9] + alt = mut_array[i, 10] + pure_tags_dict[key1] = {} + for key2, value2 in sorted(value1.items()): + for key3, value3 in value2.items(): + pure_tag = key2[:-5] + if key3 == alt: + if pure_tag in pure_tags_dict[key1]: + pure_tags_dict[key1][pure_tag] += 1 + else: + pure_tags_dict[key1][pure_tag] = 1 + + # 6. create pure_tags_dict_short with thresh + if thresh > 0: + pure_tags_dict_short = {} + for key, value in sorted(pure_tags_dict.items()): + if len(value) < thresh: + pure_tags_dict_short[key] = value + else: + pure_tags_dict_short = pure_tags_dict + + #whole_array = [] + #for k in pure_tags_dict.values(): + # whole_array.extend(k.keys()) + + # 7. output summary with threshold + workbook = xlsxwriter.Workbook(outfile) + ws1 = workbook.add_worksheet("Results") + ws2 = workbook.add_worksheet("Allele frequencies") + ws3 = workbook.add_worksheet("Tiers") + + format1 = workbook.add_format({'bg_color': '#BCF5A9'}) # green + format2 = workbook.add_format({'bg_color': '#FFC7CE'}) # red + format3 = workbook.add_format({'bg_color': '#FACC2E'}) # yellow + + header_line = ('variant ID', 'tier', 'tag', 'mate', 'read pos.ab', 'read pos.ba', 'read median length.ab', + 'read median length.ba', 'DCS median length', + 'FS.ab', 'FS.ba', 'FSqc.ab', 'FSqc.ba', 'ref.ab', 'ref.ba', 'alt.ab', 'alt.ba', + 'rel. ref.ab', 'rel. ref.ba', 'rel. alt.ab', 'rel. alt.ba', + 'na.ab', 'na.ba', 'lowq.ab', 'lowq.ba', 'trim.ab', 'trim.ba', + 'SSCS alt.ab', 'SSCS alt.ba', 'SSCS ref.ab', 'SSCS ref.ba', + 'other mut', 'chimeric tag') + ws1.write_row(0, 0, header_line) + + counter_tier11 = 0 + counter_tier12 = 0 + counter_tier21 = 0 + counter_tier22 = 0 + counter_tier23 = 0 + counter_tier24 = 0 + counter_tier31 = 0 + counter_tier32 = 0 + counter_tier41 = 0 + counter_tier42 = 0 + + if chimera_correction: + counter_tier43 = 0 + + counter_tier5 = 0 + + row = 1 + tier_dict = {} + for key1, value1 in sorted(mut_dict.items()): + counts_mut = 0 + if key1 in pure_tags_dict_short.keys(): + i = np.where(np.array(['#'.join(str(i) for i in z) + for z in zip(mut_array[:, 1], mut_array[:, 2])]) == key1)[0][0] + ref = mut_array[i, 9] + alt = mut_array[i, 10] + dcs_median = cvrg_dict[key1][2] + whole_array = pure_tags_dict_short[key1].keys() + + tier_dict[key1] = {} + if chimera_correction: + values_tier_dict = [("tier 1.1", 0), ("tier 1.2", 0), ("tier 2.1", 0), ("tier 2.2", 0), ("tier 2.3", 0), ("tier 2.4", 0), ("tier 3.1", 0), ("tier 3.2", 0), ("tier 4.1", 0), ("tier 4.2", 0), ("tier 4.3", 0), ("tier 5", 0)] + else: + values_tier_dict = [("tier 1.1", 0), ("tier 1.2", 0), ("tier 2.1", 0), ("tier 2.2", 0), ("tier 2.3", 0), ("tier 2.4", 0), ("tier 3.1", 0), ("tier 3.2", 0), ("tier 4.1", 0), ("tier 4.2", 0), ("tier 5", 0)] + + for k, v in values_tier_dict: + tier_dict[key1][k] = v + + used_keys = [] + if 'ab' in mut_pos_dict[key1].keys(): + sscs_mut_ab = mut_pos_dict[key1]['ab'] + else: + sscs_mut_ab = 0 + if 'ba' in mut_pos_dict[key1].keys(): + sscs_mut_ba = mut_pos_dict[key1]['ba'] + else: + sscs_mut_ba = 0 + if 'ab' in ref_pos_dict[key1].keys(): + sscs_ref_ab = ref_pos_dict[key1]['ab'] + else: + sscs_ref_ab = 0 + if 'ba' in ref_pos_dict[key1].keys(): + sscs_ref_ba = ref_pos_dict[key1]['ba'] + else: + sscs_ref_ba = 0 + for key2, value2 in sorted(value1.items()): + add_mut14 = "" + add_mut23 = "" + if (key2[:-5] in pure_tags_dict_short[key1].keys()) and (key2[:-5] not in used_keys) and (key1 in tag_dict[key2[:-5]].keys()): + if key2[:-5] + '.ab.1' in mut_dict[key1].keys(): + total1 = sum(mut_dict[key1][key2[:-5] + '.ab.1'].values()) + if 'na' in mut_dict[key1][key2[:-5] + '.ab.1'].keys(): + na1 = mut_dict[key1][key2[:-5] + '.ab.1']['na'] + # na1f = na1/total1 + else: + # na1 = na1f = 0 + na1 = 0 + if 'lowQ' in mut_dict[key1][key2[:-5] + '.ab.1'].keys(): + lowq1 = mut_dict[key1][key2[:-5] + '.ab.1']['lowQ'] + # lowq1f = lowq1 / total1 + else: + # lowq1 = lowq1f = 0 + lowq1 = 0 + if ref in mut_dict[key1][key2[:-5] + '.ab.1'].keys(): + ref1 = mut_dict[key1][key2[:-5] + '.ab.1'][ref] + ref1f = ref1 / (total1 - na1 - lowq1) + else: + ref1 = ref1f = 0 + if alt in mut_dict[key1][key2[:-5] + '.ab.1'].keys(): + alt1 = mut_dict[key1][key2[:-5] + '.ab.1'][alt] + alt1f = alt1 / (total1 - na1 - lowq1) + else: + alt1 = alt1f = 0 + total1new = total1 - na1 - lowq1 + if (key2[:-5] + '.ab.1') in mut_read_dict.keys(): + k1 = mut_read_dict[(key2[:-5] + '.ab.1')].keys() + add_mut1 = len(k1) + if add_mut1 > 1: + for k, v in mut_read_dict[(key2[:-5] + '.ab.1')].items(): + if k != key1: + if len(add_mut14) == 0: + add_mut14 = str(k) + "_" + v + else: + add_mut14 = add_mut14 + ", " + str(k) + "_" + v + else: + k1 = [] + else: + total1 = total1new = na1 = lowq1 = 0 + ref1 = alt1 = ref1f = alt1f = 0 + k1 = [] + + if key2[:-5] + '.ab.2' in mut_dict[key1].keys(): + total2 = sum(mut_dict[key1][key2[:-5] + '.ab.2'].values()) + if 'na' in mut_dict[key1][key2[:-5] + '.ab.2'].keys(): + na2 = mut_dict[key1][key2[:-5] + '.ab.2']['na'] + # na2f = na2 / total2 + else: + # na2 = na2f = 0 + na2 = 0 + if 'lowQ' in mut_dict[key1][key2[:-5] + '.ab.2'].keys(): + lowq2 = mut_dict[key1][key2[:-5] + '.ab.2']['lowQ'] + # lowq2f = lowq2 / total2 + else: + # lowq2 = lowq2f = 0 + lowq2 = 0 + if ref in mut_dict[key1][key2[:-5] + '.ab.2'].keys(): + ref2 = mut_dict[key1][key2[:-5] + '.ab.2'][ref] + ref2f = ref2 / (total2 - na2 - lowq2) + else: + ref2 = ref2f = 0 + if alt in mut_dict[key1][key2[:-5] + '.ab.2'].keys(): + alt2 = mut_dict[key1][key2[:-5] + '.ab.2'][alt] + alt2f = alt2 / (total2 - na2 - lowq2) + else: + alt2 = alt2f = 0 + total2new = total2 - na2 - lowq2 + if (key2[:-5] + '.ab.2') in mut_read_dict.keys(): + k2 = mut_read_dict[(key2[:-5] + '.ab.2')].keys() + add_mut2 = len(k2) + if add_mut2 > 1: + for k, v in mut_read_dict[(key2[:-5] + '.ab.2')].items(): + if k != key1: + if len(add_mut23) == 0: + add_mut23 = str(k) + "_" + v + else: + add_mut23 = add_mut23 + ", " + str(k) + "_" + v + else: + k2 = [] + else: + total2 = total2new = na2 = lowq2 = 0 + ref2 = alt2 = ref2f = alt2f = 0 + k2 = [] + + if key2[:-5] + '.ba.1' in mut_dict[key1].keys(): + total3 = sum(mut_dict[key1][key2[:-5] + '.ba.1'].values()) + if 'na' in mut_dict[key1][key2[:-5] + '.ba.1'].keys(): + na3 = mut_dict[key1][key2[:-5] + '.ba.1']['na'] + # na3f = na3 / total3 + else: + # na3 = na3f = 0 + na3 = 0 + if 'lowQ' in mut_dict[key1][key2[:-5] + '.ba.1'].keys(): + lowq3 = mut_dict[key1][key2[:-5] + '.ba.1']['lowQ'] + # lowq3f = lowq3 / total3 + else: + # lowq3 = lowq3f = 0 + lowq3 = 0 + if ref in mut_dict[key1][key2[:-5] + '.ba.1'].keys(): + ref3 = mut_dict[key1][key2[:-5] + '.ba.1'][ref] + ref3f = ref3 / (total3 - na3 - lowq3) + else: + ref3 = ref3f = 0 + if alt in mut_dict[key1][key2[:-5] + '.ba.1'].keys(): + alt3 = mut_dict[key1][key2[:-5] + '.ba.1'][alt] + alt3f = alt3 / (total3 - na3 - lowq3) + else: + alt3 = alt3f = 0 + total3new = total3 - na3 - lowq3 + if (key2[:-5] + '.ba.1') in mut_read_dict.keys(): + add_mut3 = len(mut_read_dict[(key2[:-5] + '.ba.1')].keys()) + if add_mut3 > 1: + for k, v in mut_read_dict[(key2[:-5] + '.ba.1')].items(): + if k != key1 and k not in k2: + if len(add_mut23) == 0: + add_mut23 = str(k) + "_" + v + else: + add_mut23 = add_mut23 + ", " + str(k) + "_" + v + else: + total3 = total3new = na3 = lowq3 = 0 + ref3 = alt3 = ref3f = alt3f = 0 + + if key2[:-5] + '.ba.2' in mut_dict[key1].keys(): + total4 = sum(mut_dict[key1][key2[:-5] + '.ba.2'].values()) + if 'na' in mut_dict[key1][key2[:-5] + '.ba.2'].keys(): + na4 = mut_dict[key1][key2[:-5] + '.ba.2']['na'] + # na4f = na4 / total4 + else: + # na4 = na4f = 0 + na4 = 0 + if 'lowQ' in mut_dict[key1][key2[:-5] + '.ba.2'].keys(): + lowq4 = mut_dict[key1][key2[:-5] + '.ba.2']['lowQ'] + # lowq4f = lowq4 / total4 + else: + # lowq4 = lowq4f = 0 + lowq4 = 0 + if ref in mut_dict[key1][key2[:-5] + '.ba.2'].keys(): + ref4 = mut_dict[key1][key2[:-5] + '.ba.2'][ref] + ref4f = ref4 / (total4 - na4 - lowq4) + else: + ref4 = ref4f = 0 + if alt in mut_dict[key1][key2[:-5] + '.ba.2'].keys(): + alt4 = mut_dict[key1][key2[:-5] + '.ba.2'][alt] + alt4f = alt4 / (total4 - na4 - lowq4) + else: + alt4 = alt4f = 0 + total4new = total4 - na4 - lowq4 + if (key2[:-5] + '.ba.2') in mut_read_dict.keys(): + add_mut4 = len(mut_read_dict[(key2[:-5] + '.ba.2')].keys()) + if add_mut4 > 1: + for k, v in mut_read_dict[(key2[:-5] + '.ba.2')].items(): + if k != key1 and k not in k1: + if len(add_mut14) == 0: + add_mut14 = str(k) + "_" + v + else: + add_mut14 = add_mut14 + ", " + str(k) + "_" + v + else: + total4 = total4new = na4 = lowq4 = 0 + ref4 = alt4 = ref4f = alt4f = 0 + + read_pos1 = read_pos2 = read_pos3 = read_pos4 = -1 + read_len_median1 = read_len_median2 = read_len_median3 = read_len_median4 = 0 + + if key2[:-5] + '.ab.1' in mut_read_pos_dict[key1].keys(): + read_pos1 = np.median(mut_read_pos_dict[key1][key2[:-5] + '.ab.1']) + read_len_median1 = np.median(reads_dict[key1][key2[:-5] + '.ab.1']) + if key2[:-5] + '.ab.2' in mut_read_pos_dict[key1].keys(): + read_pos2 = np.median(mut_read_pos_dict[key1][key2[:-5] + '.ab.2']) + read_len_median2 = np.median(reads_dict[key1][key2[:-5] + '.ab.2']) + if key2[:-5] + '.ba.1' in mut_read_pos_dict[key1].keys(): + read_pos3 = np.median(mut_read_pos_dict[key1][key2[:-5] + '.ba.1']) + read_len_median3 = np.median(reads_dict[key1][key2[:-5] + '.ba.1']) + if key2[:-5] + '.ba.2' in mut_read_pos_dict[key1].keys(): + read_pos4 = np.median(mut_read_pos_dict[key1][key2[:-5] + '.ba.2']) + read_len_median4 = np.median(reads_dict[key1][key2[:-5] + '.ba.2']) + + used_keys.append(key2[:-5]) + counts_mut += 1 + if (alt1f + alt2f + alt3f + alt4f) > 0.5: + if total1new == 0: + ref1f = alt1f = None + alt1ff = -1 + else: + alt1ff = alt1f + if total2new == 0: + ref2f = alt2f = None + alt2ff = -1 + else: + alt2ff = alt2f + if total3new == 0: + ref3f = alt3f = None + alt3ff = -1 + else: + alt3ff = alt3f + if total4new == 0: + ref4f = alt4f = None + alt4ff = -1 + else: + alt4ff = alt4f + + beg1 = beg4 = beg2 = beg3 = 0 + + details1 = (total1, total4, total1new, total4new, ref1, ref4, alt1, alt4, ref1f, ref4f, alt1f, alt4f, na1, na4, lowq1, lowq4, beg1, beg4) + details2 = (total2, total3, total2new, total3new, ref2, ref3, alt2, alt3, ref2f, ref3f, alt2f, alt3f, na2, na3, lowq2, lowq3, beg2, beg3) + + chrom, pos = re.split(r'\#', key1) + var_id = '-'.join([chrom, pos, ref, alt]) + sample_tag = key2[:-5] + array2 = np.unique(whole_array) # remove duplicate sequences to decrease running time + # exclude identical tag from array2, to prevent comparison to itself + same_tag = np.where(array2 == sample_tag) + index_array2 = np.arange(0, len(array2), 1) + index_withoutSame = np.delete(index_array2, same_tag) # delete identical tag from the data + array2 = array2[index_withoutSame] + if len(array2) != 0: # only perform chimera analysis if there is more than 1 variant + array1_half = sample_tag[0:int(len(sample_tag) / 2)] # mate1 part1 + array1_half2 = sample_tag[int(len(sample_tag) / 2):int(len(sample_tag))] # mate1 part 2 + array2_half = np.array([ii[0:int(len(ii) / 2)] for ii in array2]) # mate2 part1 + array2_half2 = np.array([ii[int(len(ii) / 2):int(len(ii))] for ii in array2]) # mate2 part2 + + min_tags_list_zeros = [] + chimera_tags = [] + for mate_b in [False, True]: + i = 0 # counter, only used to see how many HDs of tags were already calculated + if mate_b is False: # HD calculation for all a's + half1_mate1 = array1_half + half2_mate1 = array1_half2 + half1_mate2 = array2_half + half2_mate2 = array2_half2 + elif mate_b is True: # HD calculation for all b's + half1_mate1 = array1_half2 + half2_mate1 = array1_half + half1_mate2 = array2_half2 + half2_mate2 = array2_half + # calculate HD of "a" in the tag to all "a's" or "b" in the tag to all "b's" + dist = np.array([sum(map(operator.ne, half1_mate1, c)) for c in half1_mate2]) + min_index = np.where(dist == dist.min()) # get index of min HD + # get all "b's" of the tag or all "a's" of the tag with minimum HD + min_tag_half2 = half2_mate2[min_index] + min_tag_array2 = array2[min_index] # get whole tag with min HD + min_value = dist.min() + # calculate HD of "b" to all "b's" or "a" to all "a's" + dist_second_half = np.array([sum(map(operator.ne, half2_mate1, e)) + for e in min_tag_half2]) + + dist2 = dist_second_half.max() + max_index = np.where(dist_second_half == dist_second_half.max())[0] # get index of max HD + max_tag = min_tag_array2[max_index] + + # tags which have identical parts: + if min_value == 0 or dist2 == 0: + min_tags_list_zeros.append(tag) + chimera_tags.append(max_tag) + # chimeric = True + # else: + # chimeric = False + + # if mate_b is False: + # text = "pos {}: sample tag: {}; HD a = {}; HD b' = {}; similar tag(s): {}; chimeric = {}".format(pos, sample_tag, min_value, dist2, list(max_tag), chimeric) + # else: + # text = "pos {}: sample tag: {}; HD a' = {}; HD b = {}; similar tag(s): {}; chimeric = {}".format(pos, sample_tag, dist2, min_value, list(max_tag), chimeric) + i += 1 + chimera_tags = [x for x in chimera_tags if x != []] + chimera_tags_new = [] + for i in chimera_tags: + if len(i) > 1: + for t in i: + chimera_tags_new.append(t) + else: + chimera_tags_new.extend(i) + chimera_tags_new = np.asarray(chimera_tags_new) + chimera = ", ".join(chimera_tags_new) + else: + chimera_tags_new = [] + chimera = "" + + trimmed = False + contradictory = False + chimeric_dcs = False + + if chimera_correction and len(chimera_tags_new) > 0: # chimeras + alt1ff = 0 + alt4ff = 0 + alt2ff = 0 + alt3ff = 0 + chimeric_dcs = True + trimmed = False + contradictory = False + elif ((all(float(ij) >= 0.5 for ij in [alt1ff, alt4ff]) & # contradictory variant + all(float(ij) == 0. for ij in [alt2ff, alt3ff])) | + (all(float(ij) >= 0.5 for ij in [alt2ff, alt3ff]) & + all(float(ij) == 0. for ij in [alt1ff, alt4ff]))): + alt1ff = 0 + alt4ff = 0 + alt2ff = 0 + alt3ff = 0 + trimmed = False + contradictory = True + chimeric_dcs = False + else: + if ((read_pos1 >= 0) and ((read_pos1 <= trim) | (abs(read_len_median1 - read_pos1) <= trim))): + beg1 = total1new + total1new = 0 + alt1ff = 0 + alt1f = 0 + trimmed = True + + if ((read_pos4 >= 0) and ((read_pos4 <= trim) | (abs(read_len_median4 - read_pos4) <= trim))): + beg4 = total4new + total4new = 0 + alt4ff = 0 + alt4f = 0 + trimmed = True + + if ((read_pos2 >= 0) and ((read_pos2 <= trim) | (abs(read_len_median2 - read_pos2) <= trim))): + beg2 = total2new + total2new = 0 + alt2ff = 0 + alt2f = 0 + trimmed = True + + if ((read_pos3 >= 0) and ((read_pos3 <= trim) | (abs(read_len_median3 - read_pos3) <= trim))): + beg3 = total3new + total3new = 0 + alt3ff = 0 + alt3f = 0 + trimmed = True + details1 = (total1, total4, total1new, total4new, ref1, ref4, alt1, alt4, ref1f, ref4f, alt1f, alt4f, na1, na4, lowq1, lowq4, beg1, beg4) + details2 = (total2, total3, total2new, total3new, ref2, ref3, alt2, alt3, ref2f, ref3f, alt2f, alt3f, na2, na3, lowq2, lowq3, beg2, beg3) + + # assign tiers + if ((all(int(ij) >= 3 for ij in [total1new, total4new]) & + all(float(ij) >= 0.75 for ij in [alt1ff, alt4ff])) | + (all(int(ij) >= 3 for ij in [total2new, total3new]) & + all(float(ij) >= 0.75 for ij in [alt2ff, alt3ff]))): + tier = "1.1" + counter_tier11 += 1 + tier_dict[key1]["tier 1.1"] += 1 + + elif (all(int(ij) >= 1 for ij in [total1new, total2new, total3new, total4new]) & + any(int(ij) >= 3 for ij in [total1new, total4new]) & + any(int(ij) >= 3 for ij in [total2new, total3new]) & + all(float(ij) >= 0.75 for ij in [alt1ff, alt2ff, alt3ff, alt4ff])): + tier = "1.2" + counter_tier12 += 1 + tier_dict[key1]["tier 1.2"] += 1 + + elif ((all(int(ij) >= 1 for ij in [total1new, total4new]) & + any(int(ij) >= 3 for ij in [total1new, total4new]) & + all(float(ij) >= 0.75 for ij in [alt1ff, alt4ff])) | + (all(int(ij) >= 1 for ij in [total2new, total3new]) & + any(int(ij) >= 3 for ij in [total2new, total3new]) & + all(float(ij) >= 0.75 for ij in [alt2ff, alt3ff]))): + tier = "2.1" + counter_tier21 += 1 + tier_dict[key1]["tier 2.1"] += 1 + + elif (all(int(ij) >= 1 for ij in [total1new, total2new, total3new, total4new]) & + all(float(ij) >= 0.75 for ij in [alt1ff, alt2ff, alt3ff, alt4ff])): + tier = "2.2" + counter_tier22 += 1 + tier_dict[key1]["tier 2.2"] += 1 + + elif ((all(int(ij) >= 1 for ij in [total1new, total4new]) & + any(int(ij) >= 3 for ij in [total2new, total3new]) & + all(float(ij) >= 0.75 for ij in [alt1ff, alt4ff]) & + any(float(ij) >= 0.75 for ij in [alt2ff, alt3ff])) | + (all(int(ij) >= 1 for ij in [total2new, total3new]) & + any(int(ij) >= 3 for ij in [total1new, total4new]) & + all(float(ij) >= 0.75 for ij in [alt2ff, alt3ff]) & + any(float(ij) >= 0.75 for ij in [alt1ff, alt4ff]))): + tier = "2.3" + counter_tier23 += 1 + tier_dict[key1]["tier 2.3"] += 1 + + elif ((all(int(ij) >= 1 for ij in [total1new, total4new]) & + all(float(ij) >= 0.75 for ij in [alt1ff, alt4ff])) | + (all(int(ij) >= 1 for ij in [total2new, total3new]) & + all(float(ij) >= 0.75 for ij in [alt2ff, alt3ff]))): + tier = "2.4" + counter_tier24 += 1 + tier_dict[key1]["tier 2.4"] += 1 + + elif ((len(pure_tags_dict_short[key1]) > 1) & + (all(float(ij) >= 0.5 for ij in [alt1ff, alt4ff]) | + all(float(ij) >= 0.5 for ij in [alt2ff, alt3ff]))): + tier = "3.1" + counter_tier31 += 1 + tier_dict[key1]["tier 3.1"] += 1 + + elif ((all(int(ij) >= 1 for ij in [total1new, total4new]) & + all(float(ij) >= 0.5 for ij in [alt1ff, alt4ff])) | + (all(int(ij) >= 1 for ij in [total2new, total3new]) & + all(float(ij) >= 0.5 for ij in [alt2ff, alt3ff]))): + tier = "3.2" + counter_tier32 += 1 + tier_dict[key1]["tier 3.2"] += 1 + + elif (trimmed): + tier = "4.1" + counter_tier41 += 1 + tier_dict[key1]["tier 4.1"] += 1 + + elif (contradictory): + tier = "4.2" + counter_tier42 += 1 + tier_dict[key1]["tier 4.2"] += 1 + + elif chimera_correction and chimeric_dcs: + tier = "4.3" + counter_tier43 += 1 + tier_dict[key1]["tier 4.3"] += 1 + + else: + tier = "5" + counter_tier5 += 1 + tier_dict[key1]["tier 5"] += 1 + + if (read_pos1 == -1): + read_pos1 = read_len_median1 = None + if (read_pos4 == -1): + read_pos4 = read_len_median4 = None + if (read_pos2 == -1): + read_pos2 = read_len_median2 = None + if (read_pos3 == -1): + read_pos3 = read_len_median3 = None + line = (var_id, tier, key2[:-5], 'ab1.ba2', read_pos1, read_pos4, read_len_median1, read_len_median4, dcs_median) + details1 + (sscs_mut_ab, sscs_mut_ba, sscs_ref_ab, sscs_ref_ba, add_mut14, chimera) + ws1.write_row(row, 0, line) + line = ("", "", key2[:-5], 'ab2.ba1', read_pos2, read_pos3, read_len_median2, read_len_median3, dcs_median) + details2 + (sscs_mut_ab, sscs_mut_ba, sscs_ref_ab, sscs_ref_ba, add_mut23, chimera) + ws1.write_row(row + 1, 0, line) + + ws1.conditional_format('L{}:M{}'.format(row + 1, row + 2), + {'type': 'formula', + 'criteria': '=OR($B${}="1.1", $B${}="1.2")'.format(row + 1, row + 1), + 'format': format1, + 'multi_range': 'L{}:M{} T{}:U{} B{}'.format(row + 1, row + 2, row + 1, row + 2, row + 1, row + 2)}) + ws1.conditional_format('L{}:M{}'.format(row + 1, row + 2), + {'type': 'formula', + 'criteria': '=OR($B${}="2.1", $B${}="2.2", $B${}="2.3", $B${}="2.4")'.format(row + 1, row + 1, row + 1, row + 1), + 'format': format3, + 'multi_range': 'L{}:M{} T{}:U{} B{}'.format(row + 1, row + 2, row + 1, row + 2, row + 1, row + 2)}) + ws1.conditional_format('L{}:M{}'.format(row + 1, row + 2), + {'type': 'formula', + 'criteria': '=$B${}>="3"'.format(row + 1), + 'format': format2, + 'multi_range': 'L{}:M{} T{}:U{} B{}'.format(row + 1, row + 2, row + 1, row + 2, row + 1, row + 2)}) + + row += 3 + + if chimera_correction: + header_line2 = ('variant ID', 'cvrg', 'AC alt (all tiers)', 'AF (all tiers)', 'cvrg (tiers 1.1-2.4)', 'AC alt (tiers 1.1-2.4)', 'AF (tiers 1.1-2.4)', 'AC alt (orginal DCS)', 'AF (original DCS)', + 'tier 1.1', 'tier 1.2', 'tier 2.1', 'tier 2.2', 'tier 2.3', 'tier 2.4', + 'tier 3.1', 'tier 3.2', 'tier 4.1', 'tier 4.2', 'tier 4.3', 'tier 5', 'AF 1.1-1.2', 'AF 1.1-2.1', 'AF 1.1-2.2', + 'AF 1.1-2.3', 'AF 1.1-2.4', 'AF 1.1-3.1', 'AF 1.1-3.2', 'AF 1.1-4.1', 'AF 1.1-4.2', 'AF 1.1-4.3', 'AF 1.1-5') + else: + header_line2 = ('variant ID', 'cvrg', 'AC alt (all tiers)', 'AF (all tiers)', 'cvrg (tiers 1.1-2.4)', 'AC alt (tiers 1.1-2.4)', 'AF (tiers 1.1-2.4)', 'AC alt (orginal DCS)', 'AF (original DCS)', + 'tier 1.1', 'tier 1.2', 'tier 2.1', 'tier 2.2', 'tier 2.3', 'tier 2.4', + 'tier 3.1', 'tier 3.2', 'tier 4.1', 'tier 4.2', 'tier 5', 'AF 1.1-1.2', 'AF 1.1-2.1', 'AF 1.1-2.2', + 'AF 1.1-2.3', 'AF 1.1-2.4', 'AF 1.1-3.1', 'AF 1.1-3.2', 'AF 1.1-4.1', 'AF 1.1-4.2', 'AF 1.1-5') + + ws2.write_row(0, 0, header_line2) + row = 0 + + for key1, value1 in sorted(tier_dict.items()): + if key1 in pure_tags_dict_short.keys(): + i = np.where(np.array(['#'.join(str(i) for i in z) + for z in zip(mut_array[:, 1], mut_array[:, 2])]) == key1)[0][0] + ref = mut_array[i, 9] + alt = mut_array[i, 10] + chrom, pos = re.split(r'\#', key1) + ref_count = cvrg_dict[key1][0] + alt_count = cvrg_dict[key1][1] + cvrg = ref_count + alt_count + + var_id = '-'.join([chrom, pos, ref, alt]) + lst = [var_id, cvrg] + used_tiers = [] + cum_af = [] + for key2, value2 in sorted(value1.items()): + # calculate cummulative AF + used_tiers.append(value2) + if len(used_tiers) > 1: + cum = safe_div(sum(used_tiers), cvrg) + cum_af.append(cum) + lst.extend([sum(used_tiers), safe_div(sum(used_tiers), cvrg)]) + if chimera_correction: + lst.extend([(cvrg - sum(used_tiers[-6:])), sum(used_tiers[0:6]), safe_div(sum(used_tiers[0:6]), (cvrg - sum(used_tiers[-6:]))), alt_count, safe_div(alt_count, cvrg)]) + else: + lst.extend([(cvrg - sum(used_tiers[-5:])), sum(used_tiers[0:6]), safe_div(sum(used_tiers[0:6]), (cvrg - sum(used_tiers[-5:]))), alt_count, safe_div(alt_count, cvrg)]) + lst.extend(used_tiers) + lst.extend(cum_af) + lst = tuple(lst) + ws2.write_row(row + 1, 0, lst) + ws2.conditional_format('J{}:K{}'.format(row + 2, row + 2), {'type': 'formula', 'criteria': '=$J$1="tier 1.1"', 'format': format1, 'multi_range': 'J{}:K{} J1:K1'.format(row + 2, row + 2)}) + ws2.conditional_format('L{}:O{}'.format(row + 2, row + 2), {'type': 'formula', 'criteria': '=$L$1="tier 2.1"', 'format': format3, 'multi_range': 'L{}:O{} L1:O1'.format(row + 2, row + 2)}) + if chimera_correction: + ws2.conditional_format('P{}:U{}'.format(row + 2, row + 2), {'type': 'formula', 'criteria': '=$P$1="tier 3.1"', 'format': format2, 'multi_range': 'P{}:U{} P1:U1'.format(row + 2, row + 2)}) + else: + ws2.conditional_format('P{}:T{}'.format(row + 2, row + 2), {'type': 'formula', 'criteria': '=$P$1="tier 3.1"', 'format': format2, 'multi_range': 'P{}:T{} P1:T1'.format(row + 2, row + 2)}) + + row += 1 + + # sheet 3 + if chimera_correction: + sheet3 = [("tier 1.1", counter_tier11), ("tier 1.2", counter_tier12), ("tier 2.1", counter_tier21), + ("tier 2.2", counter_tier22), ("tier 2.3", counter_tier23), ("tier 2.4", counter_tier24), + ("tier 3.1", counter_tier31), ("tier 3.2", counter_tier32), ("tier 4.1", counter_tier41), + ("tier 4.2", counter_tier42), ("tier 4.3", counter_tier43), ("tier 5", counter_tier5)] + else: + sheet3 = [("tier 1.1", counter_tier11), ("tier 1.2", counter_tier12), ("tier 2.1", counter_tier21), + ("tier 2.2", counter_tier22), ("tier 2.3", counter_tier23), ("tier 2.4", counter_tier24), + ("tier 3.1", counter_tier31), ("tier 3.2", counter_tier32), ("tier 4.1", counter_tier41), + ("tier 4.2", counter_tier42), ("tier 5", counter_tier5)] + + + header = ("tier", "count") + ws3.write_row(0, 0, header) + + for i in range(len(sheet3)): + ws3.write_row(i + 1, 0, sheet3[i]) + ws3.conditional_format('A{}:B{}'.format(i + 2, i + 2), + {'type': 'formula', + 'criteria': '=OR($A${}="tier 1.1", $A${}="tier 1.2")'.format(i + 2, i + 2), + 'format': format1}) + ws3.conditional_format('A{}:B{}'.format(i + 2, i + 2), + {'type': 'formula', + 'criteria': '=OR($A${}="tier 2.1", $A${}="tier 2.2", $A${}="tier 2.3", $A${}="tier 2.4")'.format(i + 2, i + 2, i + 2, i + 2), + 'format': format3}) + ws3.conditional_format('A{}:B{}'.format(i + 2, i + 2), + {'type': 'formula', + 'criteria': '=$A${}>="3"'.format(i + 2), + 'format': format2}) + if chimera_correction: + description_tiers = [("Tier 1.1", "both ab and ba SSCS present (>75% of the sites with alternative base) and minimal FS>=3 for both SSCS in at least one mate"), ("", ""), ("Tier 1.2", "both ab and ba SSCS present (>75% of the sites with alt. base) and mate pair validation (min. FS=1) and minimal FS>=3 for at least one of the SSCS"), ("Tier 2.1", "both ab and ba SSCS present (>75% of the sites with alt. base) and minimal FS>=3 for at least one of the SSCS in at least one mate"), ("Tier 2.2", "both ab and ba SSCS present (>75% of the sites with alt. base) and mate pair validation (min. FS=1)"), ("Tier 2.3", "both ab and ba SSCS present (>75% of the sites with alt. base) and minimal FS=1 for both SSCS in one mate and minimal FS>=3 for at least one of the SSCS in the other mate"), ("Tier 2.4", "both ab and ba SSCS present (>75% of the sites with alt. base) and minimal FS=1 for both SSCS in at least one mate"), ("Tier 3.1", "both ab and ba SSCS present (>50% of the sites with alt. base) and recurring mutation on this position"), ("Tier 3.2", "both ab and ba SSCS present (>50% of the sites with alt. base) and minimal FS>=1 for both SSCS in at least one mate"), ("Tier 4.1", "variants at the start or end of the reads"), ("Tier 4.2", "mates with contradictory information"), ("Tier 4.3", "variants that are chimeric"), ("Tier 5", "remaining variants")] + else: + description_tiers = [("Tier 1.1", "both ab and ba SSCS present (>75% of the sites with alternative base) and minimal FS>=3 for both SSCS in at least one mate"), ("", ""), ("Tier 1.2", "both ab and ba SSCS present (>75% of the sites with alt. base) and mate pair validation (min. FS=1) and minimal FS>=3 for at least one of the SSCS"), ("Tier 2.1", "both ab and ba SSCS present (>75% of the sites with alt. base) and minimal FS>=3 for at least one of the SSCS in at least one mate"), ("Tier 2.2", "both ab and ba SSCS present (>75% of the sites with alt. base) and mate pair validation (min. FS=1)"), ("Tier 2.3", "both ab and ba SSCS present (>75% of the sites with alt. base) and minimal FS=1 for both SSCS in one mate and minimal FS>=3 for at least one of the SSCS in the other mate"), ("Tier 2.4", "both ab and ba SSCS present (>75% of the sites with alt. base) and minimal FS=1 for both SSCS in at least one mate"), ("Tier 3.1", "both ab and ba SSCS present (>50% of the sites with alt. base) and recurring mutation on this position"), ("Tier 3.2", "both ab and ba SSCS present (>50% of the sites with alt. base) and minimal FS>=1 for both SSCS in at least one mate"), ("Tier 4.1", "variants at the start or end of the reads"), ("Tier 4.2", "mates with contradictory information"), ("Tier 5", "remaining variants")] + + examples_tiers = [[("Chr5:5-20000-11068-C-G", "1.1", "AAAAAGATGCCGACTACCTT", "ab1.ba2", "254", "228", "287", "288", "289", + "3", "6", "3", "6", "0", "0", "3", "6", "0", "0", "1", "1", "0", "0", "0", "0", "0", "0", + "4081", "4098", "5", "10", "", ""), + ("", "", "AAAAAGATGCCGACTACCTT", "ab2.ba1", None, None, None, None, + "289", "0", "0", "0", "0", "0", "0", "0", "0", None, None, None, None, + "0", "0", "0", "0", "0", "0", "4081", "4098", "5", "10", "", "")], + [("Chr5:5-20000-11068-C-G", "1.1", "AAAAATGCGTAGAAATATGC", "ab1.ba2", "254", "228", "287", "288", "289", + "33", "43", "33", "43", "0", "0", "33", "43", "0", "0", "1", "1", "0", "0", "0", "0", "0", + "0", "4081", "4098", "5", "10", "", ""), + ("", "", "AAAAATGCGTAGAAATATGC", "ab2.ba1", "268", "268", "270", "288", "289", + "11", "34", "10", "27", "0", "0", "10", "27", "0", "0", "1", "1", "0", "0", "1", + "7", "0", "0", "4081", "4098", "5", "10", "", "")], + [("Chr5:5-20000-10776-G-T", "1.2", "CTATGACCCGTGAGCCCATG", "ab1.ba2", "132", "132", "287", "288", "290", + "4", "1", "4", "1", "0", "0", "4", "1", "0", "0", "1", "1", "0", "0", "0", "0", + "0", "0", "1", "6", "47170", "41149", "", ""), + ("", "", "CTATGACCCGTGAGCCCATG", "ab2.ba1", "77", "132", "233", "200", "290", + "4", "1", "4", "1", "0", "0", "4", "1", "0", "0", "1", "1", "0", "0", "0", "0", + "0", "0", "1", "6", "47170", "41149", "", "")], + [("Chr5:5-20000-11068-C-G", "2.1", "AAAAAAACATCATACACCCA", "ab1.ba2", "246", "244", "287", "288", "289", + "2", "8", "2", "8", "0", "0", "2", "8", "0", "0", "1", "1", "0", "0", "0", "0", "0", "0", + "4081", "4098", "5", "10", "", ""), + ("", "", "AAAAAAACATCATACACCCA", "ab2.ba1", None, None, None, None, + "289", "0", "0", "0", "0", "0", "0", "0", "0", None, None, None, None, "0", "0", + "0", "0", "0", "0", "4081", "4098", "5", "10", "", "")], + [("Chr5:5-20000-11068-C-G", "2.2", "ATCAGCCATGGCTATTATTG", "ab1.ba2", "72", "72", "217", "288", "289", + "1", "1", "1", "1", "0", "0", "1", "1", "0", "0", "1", "1", "0", "0", "0", "0", "0", "0", + "4081", "4098", "5", "10", "", ""), + ("", "", "ATCAGCCATGGCTATTATTG", "ab2.ba1", "153", "164", "217", "260", "289", + "1", "1", "1", "1", "0", "0", "1", "1", "0", "0", "1", "1", "0", "0", "0", "0", "0", "0", + "4081", "4098", "5", "10", "", "")], + [("Chr5:5-20000-11068-C-G", "2.3", "ATCAATATGGCCTCGCCACG", "ab1.ba2", None, None, None, None, + "289", "0", "5", "0", "5", "0", "0", "0", "5", None, None, None, "1", "0", + "0", "0", "0", "0", "0", "4081", "4098", "5", "10", "", ""), + ("", "", "ATCAATATGGCCTCGCCACG", "ab2.ba1", "202", "255", "277", "290", "289", + "1", "3", "1", "3", "0", "0", "1", "3", "0", "0", "1", "1", "0", "0", "0", "0", + "0", "0", "4081", "4098", "5", "10", "", "")], + [("Chr5:5-20000-11068-C-G", "2.4", "ATCAGCCATGGCTATTTTTT", "ab1.ba2", "72", "72", "217", "288", "289", + "1", "1", "1", "1", "0", "0", "1", "1", "0", "0", "1", "1", "0", "0", "0", "0", "0", "0", "4081", + "4098", "5", "10", "", ""), + ("", "", "ATCAGCCATGGCTATTTTTT", "ab2.ba1", "153", "164", "217", "260", "289", + "1", "1", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "1", "1", "0", "0", "0", "0", "4081", + "4098", "5", "10", "", "")], + [("Chr5:5-20000-10776-G-T", "3.1", "ATGCCTACCTCATTTGTCGT", "ab1.ba2", "46", "15", "287", "288", "290", + "3", "3", "3", "2", "3", "1", "0", "1", "1", "0.5", "0", "0.5", "0", "0", "0", "1", + "0", "0", "3", "3", "47170", "41149", "", ""), + ("", "", "ATGCCTACCTCATTTGTCGT", "ab2.ba1", None, "274", None, + "288", "290", "0", "3", "0", "2", "0", "1", "0", "1", None, "0.5", None, "0.5", + "0", "0", "0", "1", "0", "0", "3", "3", "47170", "41149", "", "")], + [("Chr5:5-20000-11315-C-T", "3.2", "ACAACATCACGTATTCAGGT", "ab1.ba2", "197", "197", "240", "255", "271", + "2", "3", "2", "3", "0", "1", "2", "2", "0", "0.333333333333333", "1", + "0.666666666666667", "0", "0", "0", "0", "0", "0", "1", "1", "6584", "6482", "", ""), + ("", "", "ACAACATCACGTATTCAGGT", "ab2.ba1", "35", "35", "240", "258", "271", + "2", "3", "2", "3", "0", "1", "2", "2", "0", "0.333333333333333", "1", + "0.666666666666667", "0", "0", "0", "0", "0", "0", "1", "1", "6584", "6482", "", "")], + [("Chr5:5-20000-13983-G-C", "4.1", "AAAAAAAGAATAACCCACAC", "ab1.ba2", "0", "100", "255", "276", "269", + "5", "6", "0", "6", "0", "0", "5", "6", "0", "0", "0", "1", "0", "0", "0", "0", "5", "0", "1", "1", "5348", "5350", "", ""), + ("", "", "AAAAAAAGAATAACCCACAC", "ab2.ba1", None, None, None, None, + "269", "0", "0", "0", "0", "0", "0", "0", "0", None, None, None, None, "0", + "0", "0", "0", "0", "0", "1", "1", "5348", "5350", "", "")], + [("Chr5:5-20000-13963-T-C", "4.2", "TTTTTAAGAATAACCCACAC", "ab1.ba2", "38", "38", "240", "283", "263", + "110", "54", "110", "54", "0", "0", "110", "54", "0", "0", "1", "1", "0", "0", "0", + "0", "0", "0", "1", "1", "5348", "5350", "", ""), + ("", "", "TTTTTAAGAATAACCCACAC", "ab2.ba1", "100", "112", "140", "145", "263", + "7", "12", "7", "12", "7", "12", "0", "0", "1", "1", "0", + "0", "0", "0", "0", "0", "0", "0", "1", "1", "5348", "5350", "", "")]] + + if chimera_correction: + examples_tiers.extend([[("Chr5:5-20000-13963-T-C", "4.3", "TTTTTAAGAAGCTATTTTTT", "ab1.ba2", "72", "72", "217", "288", "289", + "1", "1", "1", "1", "0", "0", "1", "1", "0", "0", "1", "1", "0", "0", "0", "0", "0", "0", "4081", + "4098", "5", "10", "", "TTTTTAAGAATAACCCACAC"), + ("", "", "TTTTTAAGAAGCTATTTTTT", "ab2.ba1", "153", "164", "217", "260", "289", + "1", "1", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "1", "1", "0", "0", "0", "0", "4081", + "4098", "5", "10", "", "TTTTTAAGAATAACCCACAC")], + [("Chr5:5-20000-13983-G-C", "5", "ATGTTGTGAATAACCCACAC", "ab1.ba2", None, "186", None, "276", "269", + "0", "6", "0", "6", "0", "0", "0", "6", "0", "0", "0", "1", "0", "0", "0", "0", "0", + "0", "1", "1", "5348", "5350", "", ""), + ("", "", "ATGTTGTGAATAACCCACAC", "ab2.ba1", None, None, None, None, + "269", "0", "0", "0", "0", "0", "0", "0", "0", None, None, None, None, "0", + "0", "0", "0", "0", "0", "1", "1", "5348", "5350", "", "")]]) + else: + examples_tiers.extend([ + [("Chr5:5-20000-13983-G-C", "5", "ATGTTGTGAATAACCCACAC", "ab1.ba2", None, "186", None, "276", "269", + "0", "6", "0", "6", "0", "0", "0", "6", "0", "0", "0", "1", "0", "0", "0", "0", "0", + "0", "1", "1", "5348", "5350", "", ""), + ("", "", "ATGTTGTGAATAACCCACAC", "ab2.ba1", None, None, None, None, + "269", "0", "0", "0", "0", "0", "0", "0", "0", None, None, None, None, "0", + "0", "0", "0", "0", "0", "1", "1", "5348", "5350", "", "")]]) + + start_row = 15 + ws3.write(start_row, 0, "Description of tiers with examples") + ws3.write_row(start_row + 1, 0, header_line) + row = 0 + for i in range(len(description_tiers)): + ws3.write_row(start_row + 2 + row + i + 1, 0, description_tiers[i]) + ex = examples_tiers[i] + for k in range(len(ex)): + ws3.write_row(start_row + 2 + row + i + k + 2, 0, ex[k]) + ws3.conditional_format('L{}:M{}'.format(start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 3), {'type': 'formula', 'criteria': '=OR($B${}="1.1", $B${}="1.2")'.format(start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 2), 'format': format1, 'multi_range': 'L{}:M{} T{}:U{} B{}'.format(start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 3, start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 3, start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 3)}) + ws3.conditional_format('L{}:M{}'.format(start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 3), + {'type': 'formula', 'criteria': '=OR($B${}="2.1",$B${}="2.2", $B${}="2.3", $B${}="2.4")'.format(start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 2), + 'format': format3, + 'multi_range': 'L{}:M{} T{}:U{} B{}'.format(start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 3, start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 3, start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 3)}) + ws3.conditional_format('L{}:M{}'.format(start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 3), + {'type': 'formula', + 'criteria': '=$B${}>="3"'.format(start_row + 2 + row + i + k + 2), + 'format': format2, + 'multi_range': 'L{}:M{} T{}:U{} B{}'.format(start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 3, start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 3, start_row + 2 + row + i + k + 2, start_row + 2 + row + i + k + 3)}) + row += 3 + workbook.close() + + +if __name__ == '__main__': + sys.exit(read2mut(sys.argv)) diff -r 000000000000 -r e5953c54cfb5 read2mut.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/read2mut.xml Sun Oct 04 17:19:39 2020 +0000 @@ -0,0 +1,84 @@ + + + Looks for reads with mutation at known positions and calculates frequencies and stats. + + va_macros.xml + + + xlsxwriter + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + `_. + +**Dataset 3:** JSON file generated by the **DCS mutations to tags/reads** tool +containing dictonaries of the tags of reads containing mutations +in the DCS. + +**Dataset 4:** JSON file generated by the **DCS mutations to SSCS stats** tool +stats of tags that carry a mutation in the SSCS at the same position a mutation +is called in the DCS. + +**Output** + +The output is an XLSX file containing frequencies stats for DCS mutations based +on information from the raw reads. In addition to that a tier based +classification is provided based on the amout of support for a true variant call. + + ]]> + + + diff -r 000000000000 -r e5953c54cfb5 test-data/Aligned_Families_test_data_VA.tabular --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/Aligned_Families_test_data_VA.tabular Sun Oct 04 17:19:39 2020 +0000 @@ -0,0 +1,41 @@ +GATAACCTTGCTTCGTGATTAATC ab 1 M01897:257:000000000-AYB6W:1:2112:28792:17250 2:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAATGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCGCCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGTCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC GGGGGGGGGGGGGGFGFFGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGAFGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGFGG8<:>>FD>FFC=4=:0<;DD>6461992<)892:=0306(4)-.42((44(667(449?0, +GATAACCTTGCTTCGTGATTAATC ba 2 M01897:257:000000000-AYB6W:1:1108:16316:3620 1:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCGCCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGTCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGGGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGG GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGFGGGGGGGGGGFGGGGGGGGGGGFGGGGGGFCFGGGGFGGGGEGGFGGGGGGGFFGGGGGGGGGGGGGGGGBB?CFGGGGGGCCGGFGGGGGGGCFGECC79CEECGDDD99CFGGGGF4>>GG>BDE@BBFG5=B?AF98::A423BBDFFF392?G)-96<2<:<:44<232:B3:F>6??F0(34A248:>?1,(.404-,((4( +GATAACCTTGCTTCGTGATTAATC ba 2 M01897:257:000000000-AYB6W:1:1118:5518:20674 1:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCGCCCGGGAGACCCCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGTCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCCAGAAGCGGGACGGCCGTAAGTCCCAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC FGEF9FGGCCG@FGGGGGGCFCC@EGGGCAF9E8CFCEFFFFF,?F=F8FFD=,DFFE+@CGGGFF7DFG***27)?::D)5557@>BFD@)/))).(().9<2((-29BF>F4(83,:12-)4)2,3??<<1:(7>((, +GATAACCTTGCTTCGTGATTAATC ab 2 M01897:257:000000000-AYB6W:1:2112:28792:17250 1:N:0:1 GATCCTGCCGTGTGGACTCTGTGCGGTGCCCGCAGGGCGGTGCTGGCGCTCGCCTATCGCTCTGCTCTCTCTTTGTAGACGGCGGGCGCTAACACCACCGACAAGGAGCTAGAGGTTCTCTCCTTGCACAACGTCACCTTTGAGGACGCCGGGGAGTACACCTGCCTGGCGGGCAATTCTATTGGGTTTTCTCATCACTCTGCGTGGCTGGTGGTGCTGCCAGGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCTCGGGACACG FGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGFGGGGGFFGGGGGGGGGGGGGGGGCFGGFFGGGCGGGGGGGGGCEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDFGGGGGGGGGGGGGGGGGGGGGGGGEGFGGFGFFGFGFGFFBEFGFFFF7F?FFB?DF>FDFFB:)9>FBFFF?F099E>;F0;BB1 +GATAACCTTGCTTCGTGATTAATC ba 1 M01897:257:000000000-AYB6W:1:1108:16316:3620 2:N:0:1 GATCCTGCCGTGTGGACTCTGTGCGGTGCCCGCAGGGCGGTGCTGGCGCTCGCCTATCGCTCTGCTCTCTCTTTGTAGACGGCGGGCGCTAACACCACCGACAAGGAGCTAGAGGTTCTCTCCTTGCACAACGTCACCTTTGAGGACGCCGGGGAGTACACCTGCCTGGCGGGCAATTCTATTGGGTTTTCTCATCACTCTGCGTGGCTGGTGGTGCTGCCAGGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCGCGGGACACG GGGGGGGGGGGGGGGGGFGGGGGGGGGGGFGGGGGGGGGGGGGGDGGGGGGGGGGGFGGGGGGGGFGGGGGGGGGGGGDGCEEFGGGGGCGGGGGGGGGGGGGFGGGGGGGGFGGGGGGGGGGFGGGCGFFGGGGGGGGGGFFFFGGGGCEG==CECFFCDGGGGGGGGGGGGG597*+*:=6FDFF4CFFF9B204>G?FE)5FAF?:7>FBBB69>;B))6<926(82 +GATAACCTTGCTTCGTGATTAATC ba 1 M01897:257:000000000-AYB6W:1:1118:22651:3876 2:N:0:1 GATCCTGCCGTGTGGACTCTGTGCGGTGCCCGCAGGGCGGTGCTGGCGCTCGCCTATCGCTCTGCTCTCTCTTTGTAGACGGCGGGCGCTAACACCACCGACAAGGAGCTAGAGGTTCTCTCCTTGCACAACGTCACCTTTGAGGACGCCGGGGAGTACACCTGCCTGGCGGGCAATTCTATTGGGTTTTCTCATCACTCTGCGTGGCTGGTGGTTCTGCCAGGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCCCGGGACACG FDCCF9FFDFGGGGGGGGGGGGFGGGGGGFGGGGGGGGGGGEGFECFGGGGGGGGGGGGG8EGFGEGFG,F@GGGGGEGGEG7FGGGGG@BFFGCCGGEGGGGGGGFGGGGGGFGGG@FFF9CFGGGGGGGGGGGGGGGGGGGFFGF5E*CECC>EGFGG7EGD==?E8:E7CCE3C+?:C?FFG@D3B5:>78)/C6=FFF<>B>>0:@EBF3))14>B?20>?A<2:>99>F<>>FAA9A:,403>BF?2;B(46(4((((( +GATTGGATAACGTTGTGGCAATTG ab 1 M01897:257:000000000-AYB6W:1:1101:14310:2734 2:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGGAAGTCACCGGAATCCGGGACGTCCTGGCAGCTAGGGCGGGCCCCGAGCCAGG GGGGCCFGGGGGGGGGFFGGGGGGGFGGGGEGGGGGGGFGGGGGGGGFGGGGGGGGGGGGGGGGGGGG<@FGGCCGGGGGGCFFGGGGGGGGC,@ECFBFDDGGG@FGGGGG9CFG@CCFF@DCDFC>=CEGGDEGCC@CDC*=CC*=5>FGCFEGGDFGG?:CGFD>5)/)47C@4),85:B:DF?(8)448:D:,5?7**430;>01661(-4((74,94:)(,-(-18(--(-(-(-(=01,442))(.(,(2((-(-(-((,(( +GATTGGATAACGTTGTGGCAATTG ab 1 M01897:257:000000000-AYB6W:1:1112:4840:17845 2:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCCGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCATACTGGGCACAGGGCCAGGCGTGAGGGCTCAAGAAGCGGGACCGCCGTCAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTCGGCGTGTCCCTAGCTAGC GGGC6<,CFDG8FG,CF6C,9=FFFG:@=BC7CCEFGFDGGGG788CEF66EFGGG7CF*:**=C5=FGG5AC=+:C*2:EFF*7*2/97DD>FC)7>@G@5(704(255005FFFB??FFB39((,--32()(./6>B<(())9))-38>0,43(-((((33<)-,((--(.4)).43)).( +GATTGGATAACGTTGTGGCAATTG ab 1 M01897:257:000000000-AYB6W:1:1118:8154:20084 2:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC GGGGGGGGGGGGGGFGGGGGGGGGGGGFGGGG@EDFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG=FGGGGGFFBFGGGGGGGFGGFGGGGFGGGFGGGGGEGGGGGGG:FGGGG5EGGGG:FEFGEGGGGGGGGGGGGGGGD:EG?FGFFGFCEGG>GFFFCGGFFFGEFE:>?(7.()44>B>G*=F<7:F9>D>9>F03;26:6)6>B<9(38<7A?FB2>>?FF(=:?(((.2:A:)-4((.63-,49>:?0 +GATTGGATAACGTTGTGGCAATTG ab 1 M01897:257:000000000-AYB6W:1:2110:10849:23965 2:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG@CGGGGGGGFGGGGECGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGEGGGCCCGGGGGGGGGGGGGGGGGDGGEGGGDGGGGDFGGGGGFDGDDFFFGGGFFFGFGFE@:?GFFFFFGFFFFFD2?BFFFF09>B9>F(7)2.9A2)6:44<@A7BF?>BF?>6>:((,(,5AF?F91(-:B<>,(3>00( +GATTGGATAACGTTGTGGCAATTG ba 2 M01897:257:000000000-AYB6W:1:1113:13084:11145 1:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGTACGGCCGTAAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGGCFGGGGFFFFBGFFFBEFFFFGFF?@AFFB?FFFFFFFFFFFFFFFB00:?FFFAFFFFFFFFF66>FFF<1 +GATTGGATAACGTTGTGGCAATTG ba 2 M01897:257:000000000-AYB6W:1:1115:14952:6061 1:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC GFGGGGGEGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGGGGGGGGGGGGGGFGGGGGGGGGGGGG,FGFFFFFGAFGGGGGGGGGGGGFGFGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGEGFGGGGGEGFGGGDCBEGGGGGGGGGGGGGGGGFEGGFGFGGGGGGFDGGGCDGD9DFFFGD4>FGG4FF@9DD>DD>>FFDBGFFDFFFFFFFFF=?:8?F>?F?FBFF?7>F>DB<>?B9>>?9>9?F1 +GATTGGATAACGTTGTGGCAATTG ba 2 M01897:257:000000000-AYB6W:1:2101:12835:23979 1:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGACGTAAGTCCAAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGGGTCCCGAGCCAGC GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG8FFGGGGFFGGGGGGGGGGFGGGGGGDFGGGGGFGGGGGGGGGGGGGEGG7FFGGGG@FGGGGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGFGGFGGGGC6CFFFGGCFFGGGGGG:C:47FFD3*1<677<6<;EGB@><)-3:>55-9))).:12<6)4430;>3>0(*4??F1(.:7?>(,((-.8B1999B?1 +GATTGGATAACGTTGTGGCAATTG ba 2 M01897:257:000000000-AYB6W:1:2102:25716:13556 1:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACAGGCTTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGCGCCAGC ,C@FCGGGGGE8C@FCFFGEECGCFGGDFFEFFGFFGGFGGGFFGFCF@FFAFGGGGGCGEG,EFGGG?=F@EGGGC,<=FF?DEGFFFGFFGFFDGCGG?FDDG>EFGFGA9?EFE@FGGGDFGFCFFGC+@EEE@:F:E7C1:FBCF@7<2CFFF**::8CFEEGE7C8BFF?CFGGC<9CFGCEG+CCC8:CFDCDC=:**202:65*CF5CGFD)6?5))).753>><:5>9@466-.((.9::0)4B>)8><>:0(80:2))501(--3:FF(,4>02( +GATTGGATAACGTTGTGGCAATTG ba 2 M01897:257:000000000-AYB6W:1:2107:10919:21008 1:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACCGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGCEGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG>EGEGGGFGGGGGGGGGGGGGGGGGFFGGGGGGGGGGGGGGCGGGGGGGGFGGFGGGFG@EGGBDGGFFFFGFFFFF)*4<:B@F?G6<>9BFFFFB?BFF?DF?BFF00:BFFFB;BFBB2 +GATTGGATAACGTTGTGGCAATTG ab 2 M01897:257:000000000-AYB6W:1:1101:14310:2734 1:N:0:1 GGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCGCTGTCTGGGGGACGCTGGCTCGGGACACGCCAAAGCTGCCAGGACGGACGGGAATCCTGTGACTTACGGCCGTCCCGCTTCTTGAGCCCTCACTCCTGGCCCTGTGCCCAGTGTGGGGACAAAGTTGGCCTGGCCCGGTCCTGGTCCCAGAGGGGCCCCCTCAGCCCCCTCGAGCCCACTTCCCATCTGGGTCCCCAAAGGCCTCTCCTGTGGCTCTGGCGTCTCCCGGGTGCCTGGTGGCGGTGTGGGGG GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGFFEGGGGGFFGGGGGGGGGGGGGGFFEE:=FED=FFFFFFEEEGGGGGG:FC:FFGGGGCFGGGGGGGFGFGGGCCG9FGGGGGGGGGFFEFBFGGEEEEEGDGFDC4>EDDFGGEBEFGE5>CGGFFF*)7:773(-766223:(( +GATTGGATAACGTTGTGGCAATTG ab 2 M01897:257:000000000-AYB6W:1:1112:4840:17845 1:N:0:1 GGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCTCGGGACACGCCAAAGCTGCCAGGACGGACGGGAATCCTGTGACTTACGGCCGTCCCGCTTCTTGAGCCCTCACTCCTGGCCCTGTGCCCAGTGTGGGGACAAAGTTGGCTTGGCCCGGTCCTGGCCCCAGAGGGGCCCCCTCAGCCCCCTCGAGCCCACTTCCCATCCGGGTCCCCAAAGGCCTCTCCTGTGGCTCTGGCGTCTCCCGGGTGCCTGGTGGCGGTGTGGGAC FCEFGG@@FGGGGEFGGFFGGDCF8CG86ECEGGGFGF,C,C@@FFFGGE:,CFFDG7CFFGECE=7:FCGGG:@FC6B,=EE7FFE>EGG9C?*=5CC7887*/=*:?C5E76C*::*//C>D*8D4377CF*;:?)055.547;FF?4*7F)<2@?>(41((4))8>7:0(--,-338(( +GATTGGATAACGTTGTGGCAATTG ab 2 M01897:257:000000000-AYB6W:1:1118:8154:20084 1:N:0:1 GGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCTCGGGACACGCCAAAGCTGCCAGGACGGACGGGAATCCTGTGACTTACGGCCGTCCCGCTTCTTGAGCCCTCACTCCTGGCCCTGTGCCCAGTGTGGGGACAAAGTTGGCCTGGCCCGGTCCTGGTCCCAGAGGGGCCCCCTCAGCCCCCTCGAGCCCACTTCCCATCTGGGTCCCCAAAGGCCTCTCCTGTGGCTCTGGTGTCTCCCGGGTGCCTGGTGGCGGTGTGGGAC GGGGGGGGGGGGGGGGGGGDDGGFFGGDFGGGGGGGGGGGGFEECFGGCFGEGGGGGGGGGGGGGGGFGGFCCFFFFGGGDGGGDGGGFGGGGGCFFFGGFGGGGGGGGFGC7EGG<=BFFGEGFGGGGGDGGGDGFFFGFGGGFGGFFF6@FFFFFFFF?B2>B<09B?AF:?:0(3139:FFGGFDBGFGFC?FFGABFFFB>G?AFFBFB@>>?BB>BEF?AB:??0:B;71??F?(. +GATTGGATAACGTTGTGGCAATTG ba 1 M01897:257:000000000-AYB6W:1:1113:13084:11145 2:N:0:1 GGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCTCGGGACACGCCAAAGCTGCCAGGACGGACGGGAATCCTGTGACTTACGGCCGTACCGCTTCTTGAGCCCTCACTCCTGGCCCTGTGCCCAGTGTGGGGACAAAGTTGGCCTGGCCCGGTCCTGGTCCCAGAGGGGCCCCCTCAGCCCCCTCGAGCCCACTTCCCATCTGGGTCCCCAAAGGCCTCTCCTGTGGCTCTGGTGTCTCCCGGGTGCCTGGTGGCGGTGTGGGAC GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGFGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGFCDFBGGGFFFFFFFGF:?D>FFBF>?F@FFFFFBFFFF??FFFF?62>:?FF>?FDFFFF?0:?F?ABFF>09BF?AFFF?:?>>9>?FF::0 +GATTGGATAACGTTGTGGCAATTG ba 1 M01897:257:000000000-AYB6W:1:1115:14952:6061 2:N:0:1 GGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCTCGGGACACGCCAAAGCTGCCAGGACGGACGGGAATCCTGTGACTTACGGCCGTCCCGCTTCTTGAGCCCTCACTCCTGGCCCTGTGCCCAGTGTGGGGACAAAGTTGGCCTGGCCCGGTCCTGGTCCCAGAGGGGCCCCCTCAGCCCCCTCGAGCCCACTTCCCATCTGGGTCCCCAAAGGCCTCTCCTGTGGCTCTGGTGTCTCCCGGGTGCCTGGTGGCGGTGTGGGAC GFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDGGGCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEDGF@FGGGGGGGGFGGGGGGCGGGFGGGGFGFFFCD@7DGF58FFFFFGF3:>D:6>>GBBFF474:79?28508?>>4<04>AA<09>0>F?:66C?3C.76)0319:*4)57?F*5?:7:93:B)1)21>99(489<1(((,751681(8-( +GATTGGATAACGTTGTGGCAATTG ba 1 M01897:257:000000000-AYB6W:1:2102:25716:13556 2:N:0:1 GGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTTGCTCGGGACACGCCAAAGCTGCCAGGACGGACGGGAATCCTGTGACTTACGGCCGTCCCGCCTCTTGAGCCCTCACTCCTGGCCCTGTGCCCAGTGTGGGGACAAAGTTGGCCTGGCCCGGTCCTGGTCCCAGAGGGGCCCCCTCAGCCCCCTCGCGCCCACTTCCCATCTGCGTCCCCACAGGCCTCTCCTGTTGCTCTGGTGTCTCCCGGGTGCCTGGTGGCGGTGTGGTCC ;,EACFFEGD6FFGGG8F<@BF6+,9BFA+4EDFFFEFGCC,BFCB:FFGD==C?CFG,,CEE9E7CGGE@FCCGD+8:CC<9DFGG,@:B9:F9BC@5DD5FFG;@DFBCECEE7EC7,,?CFDCDGGGFFFFC://=?=4CF4C+C558DDF5EDGBB5/*.0?@B2?04)4)1699<<:>247667<340;>B?A71((-49@(-((4()) +CCTCCCGGCAGTGCGAAAATGTCA ab 1 M01897:257:000000000-AYB6W:1:1108:17396:7377 2:N:0:1 CTAGACTCACTGGCGTTACTGACTGCGAGACCCTCCAGACAAGGCGCGTGCTGAGGTTCTGAGCCCCCTTCCGCTCCCAGTGGTGCCTGCGGCTCTGGGCCAGGGGCATCCATGGGAGCCCCGTGGGGGGGGGGGCCAGGCCGGGCCTCAACGCCCATGTCTTTGCAGCCGAGGAGGAGCTGGTGGTGGCTGCCCAGGCGGCCTGTTTTTTTGCAGGCTCCCTACGCTACGGGGTGGGCTTTTTCCGTTTCATCTTGGTGTTGCCGGCTGGGACGCCTTGCGCC GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG7FGGGGGGGGGGGCC5*:/C:*:<+2/>C:*:*+*<>?+*+0<5:/>E5<35***<6293*935=DC)))1707C5)(1*))())()*06)(((0,(*(,(,(-4(9),4D6(4((5)4*(,).2))-).5)5:228))-1(-(((((-((,()5(-( +CCTCCCGGCAGTGCGAAAATGTCA ba 2 M01897:257:000000000-AYB6W:1:1106:22053:22582 1:N:0:1 CTAGACTCACTGGCGTTACTGACTGCGAGACCCTCCAGACAAGGCGCGTGCTGAGGTTCTGAGCCCCCTTCCGCTCCCAGTGGTGCCTGCGGCTCTGGGCCAGGGGCATCCATGGGAGCCCCGTGGGGGGGGGGCCAGGCCAGGCCTCAACGCCCATGTCTTTGCAGCCGAGGAGGAGCTGGTGGAGGCTGACGACGCGGGCAGTGTGTATGCAGTCATCCTCAGCTACGGGCTGGGCTTCTTCCTGTTTATCCTGGTGGTGGCGGCTGTGTCGCTCTGCCGTC CFGGAFCFCFGGGFDGDDDGGDGGGG;F:BFGEGFGGGGFF8A:@;EFG8>:EEGE0>7+CGG>?F:?4*7?FG6).-))7)/DFGGFGGGGGGFFFFFFFFFF@FFFFCDFGF?FFAFFFDAAFFBFB9?FFD08<<6?BFFF;F?2B>9 +CCTCCCGGCAGTGCGAAAATGTCA ba 2 M01897:257:000000000-AYB6W:1:1111:28216:18792 1:N:0:1 CTAGACTCACTGGCGTTACTGACTGCGAGACCCTCCAGACAAGGCGCGTGCTGAGGTTCTGAGCCCCCTTCCGCTCCCAGTGGTGCCTGCGGCTCTGGGCCAGGGGCATCCATGGGAGCCCCGTGGGGGGGGGGCCAGGCCAGGCCTCAACGCCCATGTCTTTGCAGCCGAGGAGGAGCTGGTGGAGGCTGACGAGGCGGGCAGTGTGTATGCAGGCATCCTCAGCTACGGGGTGGGCTTCTTCCTGTTCATCCTGGTGGTGGCGGCTGTGACGCTCTGCCGCC GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG=8FGGEGDGGGGGGGGCFFGGGGGGGGGFFFFFFFFFGFFFFB58AABAFF<9<5FBF?):F:B?:2@FFFF1.54EFBFFFGFFFFF:BFFF?F?FFFFFF?FF7F?+A7FF+==CC@F?CCFFCFC@C,26,3224@C@C,,?CG+<+2CFC*:*:);C7E*21*9CE**>DDFC7+:0=/))5C)1)(*)00>*9:(.4(,577:*=47)721),,),(-(4(47()((43460(.)(0..).))).4(()(,(,)6)((((,4((((4(-(((((( +CCTCCCGGCAGTGCGAAAATGTCA ba 1 M01897:257:000000000-AYB6W:1:1109:21342:10199 2:N:0:1 CTAGGCTCTACATGGTGAGCAGAGACGAGGAGAGGGGAGCCCGCCTGGCTGCAGAGAGGGCTCACACAGCCCAGGACCAGCGTGGGCCGAGGTGGGGCTCCAGGAGGCCTGGCGGGCAGGCAGCTCAGAACCTGGTATCTACTTTCTGTTACCTGTCGCTTGAGCGGGAAGCGGGAGATCTTGTGCGCGGTGGGGGAGCCCAGGCCTTTCTTGGGGGGGCTGCGCAGGCGGCAGAGCGTCACAGCCGCCACCACCAGGATGAACAGGAAGCAGCCCAGCCCGT- GGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGFGGGGGGGGGFGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGGGGGGGGGGGGGFEGDFGFGGGGFGFGGFGGGEG?FGCDGGEGGGGGGGGG6>FEGFDFGGFFGGGEE3DFF@=@FFGF2?>FB9FFFFFBFFFBFFFFFF9>>F>F68?>>?:BABFFFFF6B??:?BF5<>BB<49?:?:?(4?:0:0(.3399 +CCTAGTCTTTGATTGGCCACTTTT ab 1 M01897:257:000000000-AYB6W:1:1106:12553:14962 2:N:0:1 TAGGCTCTACATGGTGAGCAGAGACGAGGAGAGGGGAGCCCGCCTGGCTGCAGAGAGGGCTCACACAGCCCAGGACCAGCGTGGGCCGAGGTGGGGCTCCAGGAGGCCTGGCGGGCAGGCAGCTCAGAACCTGGTATCTACTTTCTGTTACCTGTCGCTTGAGCGGGAAGCGGGAGATCTTGTGCACGGTGGGGGAGCCCAGGCCTTTCTTGGGGGGGCTGCGCAGGCGGCAGAGCGTCACAGCCGCCACCACCAGGATGAACAGGAAGAAGCCCACCCCGT-- GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEEEGGGGGGGGFGGGGGGFGGFFGGD@9DFFFGGFFFFF7AFFGFFFFFFFFFFFFFFBFF6>FFFFBBB>FBBADBFFF((428?F:?:DFA:1:?FF7??F10(3755>:DF>AFFGA=>FF>>09>09@6>>BFBF(3.116>F:0)-:1)37(:>)4:D>F6:AFBBFFFFFF0(8,1<6?7(6FGACCCCFGGE7EGGGGGGGFGCCFF;@A<;E;CFGFFFFG5:*=+588C8>57EEGGDEFGGGCFC?*9@FGGGB>)97)<=@F?/))3:BGFFFE:)5925?>326909>>6>54,7((-8)-8-71(--24641:B)47445270,(3124(.,(,:<>(. +CCTAGTCTTTGATTGGCCACTTTT ab 2 M01897:257:000000000-AYB6W:1:1106:12553:14962 1:N:0:1 CTAGACTCACTGGCGTTACTGACTGCGAGACCCTCCAGACAAGGCGCGTGCTGAGGTTCTGAGCCCCCTTCCGCTCCCAGTGGTGCCTGCGGCTCTGGGCCAGGGGCATCCATGGGAGCCCCGTGGGGGGGGGGG--CCAGGCAGGGCCCCAAGCCCCTTGTCT-TGCAGCCGGGGGGGGGGCGGTGGGAGCCTAACAAGCGGGGCGGGGGGTTGGAGGCCTCCCCAAGTTCGGGGGTGGCTTCTTCCTGTTCATCCTTGGTGTGGGGGCTGTGACGCCTTTGCGGC GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG* ?8*;*;**:*;***2***2A***00+< C++0++;***:**:*****://:**;**0++*++2*/:E/*1**)))/)1)+*1**))9))**)/**)03>))8D)(8().5<*)-7))1)67)6/.8118((4(-,.()-(()(-)).(,- +CCTAGTCTTTGATTGGCCACTTTT ab 2 M01897:257:000000000-AYB6W:1:1106:15615:18803 1:N:0:1 CTAGACTCACTGGCGTTACTGACTGCGAGACCCTCCAGACAAGGCGCGTGCTGAGGTTCTGAGCCCCCTTCCGCTCCCAGTGGTGCCTGCGGCTCTGGGCCAGGGGCATCCATGGGAGCCCCGTGGGGGGGGGGA--CCAGGCCAGGCCTCAACGCCCATGTCTTTGCAGCCGAGGAGGAGCTGGTGGAGGCTGACGAGGCGGGCAGTGTGTATGCAGGCATCCTCAGCTACGGGGTGGGCTTCTTCCTGTTCATCCTGGTGGTGGCGGCTGTGACGCTCTGCC-GC GGGGGFGGGGFGGGGGGGGGFGGGGGGGGGGEGGGGGGGGGDEFGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGCFGGGGGGGCFGEGFGEEGGGGGGGDGGDGGCFGGGGGGGDGGGGGGFGGGGGGGGGGG GDGGGGGGGFGGGGFGGGGGGGGGFGG9FFGGGGGGGGGCEGECGGCFCEEGGGGGGGGGGGCGGCEC3*:C>DG=FCB0?B><:D243847 10 +CCTAGTCTTTGATTGGCCACTTTT ab 2 M01897:257:000000000-AYB6W:1:1110:11692:17499 1:N:0:1 CTAGACTCACTGGCGTTACTGACTGCGAGACCCTCCAGACAAGGCGCGTGCTGAGGTTCTGAGCCCCCTTCCGCTCCCAGTGGTGCCTGCGGCTCTGGGCCAGGGGCATCCATGGGAGCCCCGTGGGGGGGGGGGGGCCAGGCCAGGCCCCAACGCCCATGTCTTTGCAGCCGAGGGGGAGCTGGTTGGGGCTGACGAGGCGGGCAGTGGGTATGCAGGCATCCTCAGCTACGGGGTGGGCTTCTTCCTGTTCATCCTGGGGGTGGCGGCTTGTACCCTCTTCC--- GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG8>5*;:8C8E@;?*:;88*2CE*8*++3/4;<31<2)9:4=).0))/69?<((213(7:960(,1.-))))(()-) +GATAAGCCAACTGCCATCTAGAAT ab 1 M01897:257:000000000-AYB6W:1:1105:25798:19415 2:N:0:1 CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCGCCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGTCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGTACGTCACAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTTTCCCGAGCCAGC GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGGGGGFGGGFGGGGGGGGGGGGGGEFCFFGGDEGG8EGGGGFGGGGEGFEGFFGFFCDBF)7:>7FF:EF?B01>BF;FFF*4(,2:24?FBBF>?F?FFBF0;B2:0(:??FF7:BF?03:2>CEFDGB7B5/<<:9<>>>>>:>:D(47:6<26)402346>2<>(-49??0 +GATAAGCCAACTGCCATCTAGAAT ab 2 M01897:257:000000000-AYB6W:1:1105:25798:19415 1:N:0:1 GATCCTGCCGTGTGGACTCTGTGCGGTGCCCGCAGGGCGGTGCTGGCGCTCGCCTATCGCTCTGCTCTCTCTTTGTAGACGGCGGGCGCTAACACCACCGACAAGGAGCTAGAGGTTCTCTCCTTGCACAACGTCACCTTTGAGGACGCCGGGGAGTACACCTGCCTGGCGGGCAATTCTATTGGGTTTTCTCATCACTCTGCGTGGCTGGTGGTGCTGCCAGGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCTCGGGACACG GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFFFFGGFFGDB9FFGFFFF0F?:?FFFFFFFFFFFFBF@FFBA?B9;9B9>BB>FFF>FF>><:>>FD>FFC=4=:0<;DD>6461992<)892:=0306(4)-.42((44(667(449?0, +@GATAACCTTGCTTCGTGATTAATC.ba.2 +CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCGCCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGTCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGGGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGG ++ +GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGFGGGGGGGGGGFGGGGGGGGGGGFGGGGGGFCFGGGGFGGGGEGGFGGGGGGGFFGGGGGGGGGGGGGGGGBB?CFGGGGGGCCGGFGGGGGGGCFGECC79CEECGDDD99CFGGGGF4>>GG>BDE@BBFG5=B?AF98::A423BBDFFF392?G)-96<2<:<:44<232:B3:F>6??F0(34A248:>?1,(.404-,((4( +@GATAACCTTGCTTCGTGATTAATC.ba.2 +CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCGCCCGGGAGACCCCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGTCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCCAGAAGCGGGACGGCCGTAAGTCCCAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC ++ +FGEF9FGGCCG@FGGGGGGCFCC@EGGGCAF9E8CFCEFFFFF,?F=F8FFD=,DFFE+@CGGGFF7DFG***27)?::D)5557@>BFD@)/))).(().9<2((-29BF>F4(83,:12-)4)2,3??<<1:(7>((, +@GATAACCTTGCTTCGTGATTAATC.ab.2 +GATCCTGCCGTGTGGACTCTGTGCGGTGCCCGCAGGGCGGTGCTGGCGCTCGCCTATCGCTCTGCTCTCTCTTTGTAGACGGCGGGCGCTAACACCACCGACAAGGAGCTAGAGGTTCTCTCCTTGCACAACGTCACCTTTGAGGACGCCGGGGAGTACACCTGCCTGGCGGGCAATTCTATTGGGTTTTCTCATCACTCTGCGTGGCTGGTGGTGCTGCCAGGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCTCGGGACACG ++ +FGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGFGGGGGFFGGGGGGGGGGGGGGGGCFGGFFGGGCGGGGGGGGGCEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDFGGGGGGGGGGGGGGGGGGGGGGGGEGFGGFGFFGFGFGFFBEFGFFFF7F?FFB?DF>FDFFB:)9>FBFFF?F099E>;F0;BB1 +@GATAACCTTGCTTCGTGATTAATC.ba.1 +GATCCTGCCGTGTGGACTCTGTGCGGTGCCCGCAGGGCGGTGCTGGCGCTCGCCTATCGCTCTGCTCTCTCTTTGTAGACGGCGGGCGCTAACACCACCGACAAGGAGCTAGAGGTTCTCTCCTTGCACAACGTCACCTTTGAGGACGCCGGGGAGTACACCTGCCTGGCGGGCAATTCTATTGGGTTTTCTCATCACTCTGCGTGGCTGGTGGTGCTGCCAGGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCGCGGGACACG ++ +GGGGGGGGGGGGGGGGGFGGGGGGGGGGGFGGGGGGGGGGGGGGDGGGGGGGGGGGFGGGGGGGGFGGGGGGGGGGGGDGCEEFGGGGGCGGGGGGGGGGGGGFGGGGGGGGFGGGGGGGGGGFGGGCGFFGGGGGGGGGGFFFFGGGGCEG==CECFFCDGGGGGGGGGGGGG597*+*:=6FDFF4CFFF9B204>G?FE)5FAF?:7>FBBB69>;B))6<926(82 +@GATAACCTTGCTTCGTGATTAATC.ba.1 +GATCCTGCCGTGTGGACTCTGTGCGGTGCCCGCAGGGCGGTGCTGGCGCTCGCCTATCGCTCTGCTCTCTCTTTGTAGACGGCGGGCGCTAACACCACCGACAAGGAGCTAGAGGTTCTCTCCTTGCACAACGTCACCTTTGAGGACGCCGGGGAGTACACCTGCCTGGCGGGCAATTCTATTGGGTTTTCTCATCACTCTGCGTGGCTGGTGGTTCTGCCAGGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCCCGGGACACG ++ +FDCCF9FFDFGGGGGGGGGGGGFGGGGGGFGGGGGGGGGGGEGFECFGGGGGGGGGGGGG8EGFGEGFG,F@GGGGGEGGEG7FGGGGG@BFFGCCGGEGGGGGGGFGGGGGGFGGG@FFF9CFGGGGGGGGGGGGGGGGGGGFFGF5E*CECC>EGFGG7EGD==?E8:E7CCE3C+?:C?FFG@D3B5:>78)/C6=FFF<>B>>0:@EBF3))14>B?20>?A<2:>99>F<>>FAA9A:,403>BF?2;B(46(4((((( +@GATTGGATAACGTTGTGGCAATTG.ab.1 +CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGGAAGTCACCGGAATCCGGGACGTCCTGGCAGCTAGGGCGGGCCCCGAGCCAGG ++ +GGGGCCFGGGGGGGGGFFGGGGGGGFGGGGEGGGGGGGFGGGGGGGGFGGGGGGGGGGGGGGGGGGGG<@FGGCCGGGGGGCFFGGGGGGGGC,@ECFBFDDGGG@FGGGGG9CFG@CCFF@DCDFC>=CEGGDEGCC@CDC*=CC*=5>FGCFEGGDFGG?:CGFD>5)/)47C@4),85:B:DF?(8)448:D:,5?7**430;>01661(-4((74,94:)(,-(-18(--(-(-(-(=01,442))(.(,(2((-(-(-((,(( +@GATTGGATAACGTTGTGGCAATTG.ab.1 +CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCCGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCATACTGGGCACAGGGCCAGGCGTGAGGGCTCAAGAAGCGGGACCGCCGTCAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTCGGCGTGTCCCTAGCTAGC ++ +GGGC6<,CFDG8FG,CF6C,9=FFFG:@=BC7CCEFGFDGGGG788CEF66EFGGG7CF*:**=C5=FGG5AC=+:C*2:EFF*7*2/97DD>FC)7>@G@5(704(255005FFFB??FFB39((,--32()(./6>B<(())9))-38>0,43(-((((33<)-,((--(.4)).43)).( +@GATTGGATAACGTTGTGGCAATTG.ab.1 +CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC ++ +GGGGGGGGGGGGGGFGGGGGGGGGGGGFGGGG@EDFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG=FGGGGGFFBFGGGGGGGFGGFGGGGFGGGFGGGGGEGGGGGGG:FGGGG5EGGGG:FEFGEGGGGGGGGGGGGGGGD:EG?FGFFGFCEGG>GFFFCGGFFFGEFE:>?(7.()44>B>G*=F<7:F9>D>9>F03;26:6)6>B<9(38<7A?FB2>>?FF(=:?(((.2:A:)-4((.63-,49>:?0 +@GATTGGATAACGTTGTGGCAATTG.ab.1 +CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC ++ +GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG@CGGGGGGGFGGGGECGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGEGGGCCCGGGGGGGGGGGGGGGGGDGGEGGGDGGGGDFGGGGGFDGDDFFFGGGFFFGFGFE@:?GFFFFFGFFFFFD2?BFFFF09>B9>F(7)2.9A2)6:44<@A7BF?>BF?>6>:((,(,5AF?F91(-:B<>,(3>00( +@GATTGGATAACGTTGTGGCAATTG.ba.2 +CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGTACGGCCGTAAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC ++ +GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGGCFGGGGFFFFBGFFFBEFFFFGFF?@AFFB?FFFFFFFFFFFFFFFB00:?FFFAFFFFFFFFF66>FFF<1 +@GATTGGATAACGTTGTGGCAATTG.ba.2 +CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC ++ +GFGGGGGEGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGGGGGGGGGGGGGGFGGGGGGGGGGGGG,FGFFFFFGAFGGGGGGGGGGGGFGFGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGEGFGGGGGEGFGGGDCBEGGGGGGGGGGGGGGGGFEGGFGFGGGGGGFDGGGCDGD9DFFFGD4>FGG4FF@9DD>DD>>FFDBGFFDFFFFFFFFF=?:8?F>?F?FBFF?7>F>DB<>?B9>>?9>9?F1 +@GATTGGATAACGTTGTGGCAATTG.ba.2 +CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGACGTAAGTCCAAGGATTCCCGTCCGTCCTGGCAGCTTTGGCGGGTCCCGAGCCAGC ++ +GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG8FFGGGGFFGGGGGGGGGGFGGGGGGDFGGGGGFGGGGGGGGGGGGGEGG7FFGGGG@FGGGGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGFGGFGGGGC6CFFFGGCFFGGGGGG:C:47FFD3*1<677<6<;EGB@><)-3:>55-9))).:12<6)4430;>3>0(*4??F1(.:7?>(,((-.8B1999B?1 +@GATTGGATAACGTTGTGGCAATTG.ba.2 +CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACAGGCTTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGCGCCAGC ++ +,C@FCGGGGGE8C@FCFFGEECGCFGGDFFEFFGFFGGFGGGFFGFCF@FFAFGGGGGCGEG,EFGGG?=F@EGGGC,<=FF?DEGFFFGFFGFFDGCGG?FDDG>EFGFGA9?EFE@FGGGDFGFCFFGC+@EEE@:F:E7C1:FBCF@7<2CFFF**::8CFEEGE7C8BFF?CFGGC<9CFGCEG+CCC8:CFDCDC=:**202:65*CF5CGFD)6?5))).753>><:5>9@466-.((.9::0)4B>)8><>:0(80:2))501(--3:FF(,4>02( +@GATTGGATAACGTTGTGGCAATTG.ba.2 +CTAGAGGGCCAGACCCTGGAGAGAAGGAGCCCAGCAGAGCCAGCCAGTCCCACACCGCCACCAGGCACCCGGGAGACACCAGAGCCACAGGAGAGGCCTTTGGGGACCCAGATGGGAAGTGGGCTCGAGGGGGCTGAGGGGGCCCCTCTGGGACCAGGACCGGGCCAGGCCAACTTTGTCCCCACACTGGGCACAGGGCCAGGAGTGAGGGCTCAAGAAGCGGGACGGCCGTAAGTCACCGGATTCCCGTCCGTCCTGGCAGCTTTGGCGTGTCCCGAGCCAGC ++ +GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGCEGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG>EGEGGGFGGGGGGGGGGGGGGGGGFFGGGGGGGGGGGGGGCGGGGGGGGFGGFGGGFG@EGGBDGGFFFFGFFFFF)*4<:B@F?G6<>9BFFFFB?BFF?DF?BFF00:BFFFB;BFBB2 +@GATTGGATAACGTTGTGGCAATTG.ab.2 +GGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCGCTGTCTGGGGGACGCTGGCTCGGGACACGCCAAAGCTGCCAGGACGGACGGGAATCCTGTGACTTACGGCCGTCCCGCTTCTTGAGCCCTCACTCCTGGCCCTGTGCCCAGTGTGGGGACAAAGTTGGCCTGGCCCGGTCCTGGTCCCAGAGGGGCCCCCTCAGCCCCCTCGAGCCCACTTCCCATCTGGGTCCCCAAAGGCCTCTCCTGTGGCTCTGGCGTCTCCCGGGTGCCTGGTGGCGGTGTGGGGG ++ +GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGFFEGGGGGFFGGGGGGGGGGGGGGFFEE:=FED=FFFFFFEEEGGGGGG:FC:FFGGGGCFGGGGGGGFGFGGGCCG9FGGGGGGGGGFFEFBFGGEEEEEGDGFDC4>EDDFGGEBEFGE5>CGGFFF*)7:773(-766223:(( +@GATTGGATAACGTTGTGGCAATTG.ab.2 +GGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCTCGGGACACGCCAAAGCTGCCAGGACGGACGGGAATCCTGTGACTTACGGCCGTCCCGCTTCTTGAGCCCTCACTCCTGGCCCTGTGCCCAGTGTGGGGACAAAGTTGGCTTGGCCCGGTCCTGGCCCCAGAGGGGCCCCCTCAGCCCCCTCGAGCCCACTTCCCATCCGGGTCCCCAAAGGCCTCTCCTGTGGCTCTGGCGTCTCCCGGGTGCCTGGTGGCGGTGTGGGAC ++ +FCEFGG@@FGGGGEFGGFFGGDCF8CG86ECEGGGFGF,C,C@@FFFGGE:,CFFDG7CFFGECE=7:FCGGG:@FC6B,=EE7FFE>EGG9C?*=5CC7887*/=*:?C5E76C*::*//C>D*8D4377CF*;:?)055.547;FF?4*7F)<2@?>(41((4))8>7:0(--,-338(( +@GATTGGATAACGTTGTGGCAATTG.ab.2 +GGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCTCGGGACACGCCAAAGCTGCCAGGACGGACGGGAATCCTGTGACTTACGGCCGTCCCGCTTCTTGAGCCCTCACTCCTGGCCCTGTGCCCAGTGTGGGGACAAAGTTGGCCTGGCCCGGTCCTGGTCCCAGAGGGGCCCCCTCAGCCCCCTCGAGCCCACTTCCCATCTGGGTCCCCAAAGGCCTCTCCTGTGGCTCTGGTGTCTCCCGGGTGCCTGGTGGCGGTGTGGGAC ++ +GGGGGGGGGGGGGGGGGGGDDGGFFGGDFGGGGGGGGGGGGFEECFGGCFGEGGGGGGGGGGGGGGGFGGFCCFFFFGGGDGGGDGGGFGGGGGCFFFGGFGGGGGGGGFGC7EGG<=BFFGEGFGGGGGDGGGDGFFFGFGGGFGGFFF6@FFFFFFFF?B2>B<09B?AF:?:0(3139:FFGGFDBGFGFC?FFGABFFFB>G?AFFBFB@>>?BB>BEF?AB:??0:B;71??F?(. +@GATTGGATAACGTTGTGGCAATTG.ba.1 +GGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCTCGGGACACGCCAAAGCTGCCAGGACGGACGGGAATCCTGTGACTTACGGCCGTACCGCTTCTTGAGCCCTCACTCCTGGCCCTGTGCCCAGTGTGGGGACAAAGTTGGCCTGGCCCGGTCCTGGTCCCAGAGGGGCCCCCTCAGCCCCCTCGAGCCCACTTCCCATCTGGGTCCCCAAAGGCCTCTCCTGTGGCTCTGGTGTCTCCCGGGTGCCTGGTGGCGGTGTGGGAC ++ +GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGFGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGFCDFBGGGFFFFFFFGF:?D>FFBF>?F@FFFFFBFFFF??FFFF?62>:?FF>?FDFFFF?0:?F?ABFF>09BF?AFFF?:?>>9>?FF::0 +@GATTGGATAACGTTGTGGCAATTG.ba.1 +GGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTGGCTCGGGACACGCCAAAGCTGCCAGGACGGACGGGAATCCTGTGACTTACGGCCGTCCCGCTTCTTGAGCCCTCACTCCTGGCCCTGTGCCCAGTGTGGGGACAAAGTTGGCCTGGCCCGGTCCTGGTCCCAGAGGGGCCCCCTCAGCCCCCTCGAGCCCACTTCCCATCTGGGTCCCCAAAGGCCTCTCCTGTGGCTCTGGTGTCTCCCGGGTGCCTGGTGGCGGTGTGGGAC ++ +GFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDGGGCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEDGF@FGGGGGGGGFGGGGGGCGGGFGGGGFGFFFCD@7DGF58FFFFFGF3:>D:6>>GBBFF474:79?28508?>>4<04>AA<09>0>F?:66C?3C.76)0319:*4)57?F*5?:7:93:B)1)21>99(489<1(((,751681(8-( +@GATTGGATAACGTTGTGGCAATTG.ba.1 +GGTACCGGCTTCTGCTGCTGCTGCTGCTCCGCACTGTCTGGGGGACGCTTGCTCGGGACACGCCAAAGCTGCCAGGACGGACGGGAATCCTGTGACTTACGGCCGTCCCGCCTCTTGAGCCCTCACTCCTGGCCCTGTGCCCAGTGTGGGGACAAAGTTGGCCTGGCCCGGTCCTGGTCCCAGAGGGGCCCCCTCAGCCCCCTCGCGCCCACTTCCCATCTGCGTCCCCACAGGCCTCTCCTGTTGCTCTGGTGTCTCCCGGGTGCCTGGTGGCGGTGTGGTCC ++ +;,EACFFEGD6FFGGG8F<@BF6+,9BFA+4EDFFFEFGCC,BFCB:FFGD==C?CFG,,CEE9E7CGGE@FCCGD+8:CC<9DFGG,@:B9:F9BC@5DD5FFG;@DFBCECEE7EC7,,?CFDCDGGGFFFFC://=?=4CF4C+C558DDF5EDGBB5/*.0?@B2?04)4)1699<<:>247667<340;>B?A71((-49@(-((4()) +@CCTCCCGGCAGTGCGAAAATGTCA.ab.1 +CTAGACTCACTGGCGTTACTGACTGCGAGACCCTCCAGACAAGGCGCGTGCTGAGGTTCTGAGCCCCCTTCCGCTCCCAGTGGTGCCTGCGGCTCTGGGCCAGGGGCATCCATGGGAGCCCCGTGGGGGGGGGGGCCAGGCCGGGCCTCAACGCCCATGTCTTTGCAGCCGAGGAGGAGCTGGTGGTGGCTGCCCAGGCGGCCTGTTTTTTTGCAGGCTCCCTACGCTACGGGGTGGGCTTTTTCCGTTTCATCTTGGTGTTGCCGGCTGGGACGCCTTGCGCC ++ +GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG7FGGGGGGGGGGGCC5*:/C:*:<+2/>C:*:*+*<>?+*+0<5:/>E5<35***<6293*935=DC)))1707C5)(1*))())()*06)(((0,(*(,(,(-4(9),4D6(4((5)4*(,).2))-).5)5:228))-1(-(((((-((,()5(-( +@CCTCCCGGCAGTGCGAAAATGTCA.ba.2 +CTAGACTCACTGGCGTTACTGACTGCGAGACCCTCCAGACAAGGCGCGTGCTGAGGTTCTGAGCCCCCTTCCGCTCCCAGTGGTGCCTGCGGCTCTGGGCCAGGGGCATCCATGGGAGCCCCGTGGGGGGGGGGCCAGGCCAGGCCTCAACGCCCATGTCTTTGCAGCCGAGGAGGAGCTGGTGGAGGCTGACGACGCGGGCAGTGTGTATGCAGTCATCCTCAGCTACGGGCTGGGCTTCTTCCTGTTTATCCTGGTGGTGGCGGCTGTGTCGCTCTGCCGTC ++ +CFGGAFCFCFGGGFDGDDDGGDGGGG;F:BFGEGFGGGGFF8A:@;EFG8>:EEGE0>7+CGG>?F:?4*7?FG6).-))7)/DFGGFGGGGGGFFFFFFFFFF@FFFFCDFGF?FFAFFFDAAFFBFB9?FFD08<<6?BFFF;F?2B>9 +@CCTCCCGGCAGTGCGAAAATGTCA.ba.2 +CTAGACTCACTGGCGTTACTGACTGCGAGACCCTCCAGACAAGGCGCGTGCTGAGGTTCTGAGCCCCCTTCCGCTCCCAGTGGTGCCTGCGGCTCTGGGCCAGGGGCATCCATGGGAGCCCCGTGGGGGGGGGGCCAGGCCAGGCCTCAACGCCCATGTCTTTGCAGCCGAGGAGGAGCTGGTGGAGGCTGACGAGGCGGGCAGTGTGTATGCAGGCATCCTCAGCTACGGGGTGGGCTTCTTCCTGTTCATCCTGGTGGTGGCGGCTGTGACGCTCTGCCGCC ++ +GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG=8FGGEGDGGGGGGGGCFFGGGGGGGGGFFFFFFFFFGFFFFB58AABAFF<9<5FBF?):F:B?:2@FFFF1.54EFBFFFGFFFFF:BFFF?F?FFFFFF?FF7F?+A7FF+==CC@F?CCFFCFC@C,26,3224@C@C,,?CG+<+2CFC*:*:);C7E*21*9CE**>DDFC7+:0=/))5C)1)(*)00>*9:(.4(,577:*=47)721),,),(-(4(47()((43460(.)(0..).))).4(()(,(,)6)((((,4((((4(-(((((( +@CCTCCCGGCAGTGCGAAAATGTCA.ba.1 +CTAGGCTCTACATGGTGAGCAGAGACGAGGAGAGGGGAGCCCGCCTGGCTGCAGAGAGGGCTCACACAGCCCAGGACCAGCGTGGGCCGAGGTGGGGCTCCAGGAGGCCTGGCGGGCAGGCAGCTCAGAACCTGGTATCTACTTTCTGTTACCTGTCGCTTGAGCGGGAAGCGGGAGATCTTGTGCGCGGTGGGGGAGCCCAGGCCTTTCTTGGGGGGGCTGCGCAGGCGGCAGAGCGTCACAGCCGCCACCACCAGGATGAACAGGAAGCAGCCCAGCCCGT ++ +GGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGFGGGGGGGGGFGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGGGGGGGGGGGGGFEGDFGFGGGGFGFGGFGGGEG?FGCDGGEGGGGGGGGG6>FEGFDFGGFFGGGEE3DFF@=@FFGF2?>FB9FFFFFBFFFBFFFFFF9>>F>F68?>>?:BABFFFFF6B??:?BF5<>BB<49?:?:?(4?:0:0(.3399 diff -r 000000000000 -r e5953c54cfb5 test-data/Interesting_Reads_test_data_VA.trim.bam Binary file test-data/Interesting_Reads_test_data_VA.trim.bam has changed diff -r 000000000000 -r e5953c54cfb5 test-data/Interesting_Reads_test_data_VA.trim.bam.bai Binary file test-data/Interesting_Reads_test_data_VA.trim.bam.bai has changed diff -r 000000000000 -r e5953c54cfb5 test-data/SSCS_counts_test_data_VA.json --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/SSCS_counts_test_data_VA.json Sun Oct 04 17:19:39 2020 +0000 @@ -0,0 +1,1 @@ +[{"ACH_TDII_5regions#505": {"ab": 2, "ba": 1}, "ACH_TDII_5regions#571": {"ab": 1, "ba": 1}, "ACH_TDII_5regions#958": {"ab": 1, "ba": 1}}, {"ACH_TDII_5regions#505": {"ab": 1, "ba": 1}, "ACH_TDII_5regions#571": {"ab": 2, "ba": 1}, "ACH_TDII_5regions#958": {"ab": 1}}] \ No newline at end of file diff -r 000000000000 -r e5953c54cfb5 test-data/SSCS_test_data_VA.bam Binary file test-data/SSCS_test_data_VA.bam has changed diff -r 000000000000 -r e5953c54cfb5 test-data/SSCS_test_data_VA.bam.bai Binary file test-data/SSCS_test_data_VA.bam.bai has changed diff -r 000000000000 -r e5953c54cfb5 test-data/mutant_reads_summary_short_trim_test_data_VA.xlsx Binary file test-data/mutant_reads_summary_short_trim_test_data_VA.xlsx has changed diff -r 000000000000 -r e5953c54cfb5 test-data/tag_count_dict_test_data_VA.json --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/tag_count_dict_test_data_VA.json Sun Oct 04 17:19:39 2020 +0000 @@ -0,0 +1,1 @@ +[{"GATAACCTTGCTTCGTGATTAATC": {"ACH_TDII_5regions#505": "A"}, "GATTGGATAACGTTGTGGCAATTG": {"ACH_TDII_5regions#571": "T"}, "CCTCCCGGCAGTGCGAAAATGTCA": {"ACH_TDII_5regions#958": "C"}}, {"ACH_TDII_5regions#505": [1, 1, 173.0], "ACH_TDII_5regions#571": [1, 1, 143.0], "ACH_TDII_5regions#958": [0, 1, 195.0]}] \ No newline at end of file diff -r 000000000000 -r e5953c54cfb5 va_macros.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/va_macros.xml Sun Oct 04 17:19:39 2020 +0000 @@ -0,0 +1,20 @@ + + + + +@misc{duplex, + author = {Povysil, Gundula and Heinzl, Monika and Salazar, Renato and Stoler, Nicholas and Nekrutenko, Anton and Tiemann-Boege, Irene}, + year = {2019}, + title = {{Variant Analyzer: a quality control for variant calling in duplex sequencing data (manuscript)}} + } + + + + + + matplotlib + pysam + + + +