Mercurial > repos > mheinzl > variant_analyzer2
changeset 78:fdfe9a919ff7 draft
planemo upload for repository https://github.com/Single-Molecule-Genetics/VariantAnalyzerGalaxy/tree/master/tools/variant_analyzer commit ee4a8e6cf290e6c8a4d55f9cd2839d60ab3b11c8-dirty
| author | mheinzl |
|---|---|
| date | Fri, 22 Jul 2022 09:19:44 +0000 |
| parents | 1797e461d674 |
| children | d7aea14291e8 |
| files | mut2read.py mut2read.xml mut2sscs.py mut2sscs.xml read2mut.py read2mut.xml test-data/Interesting_Reads_test.trim.bai test-data/Interesting_Reads_test.trim.bam test-data/Variant_Analyzer_allele_frequencies_test.xlsx test-data/Variant_Analyzer_summary_test.csv test-data/Variant_Analyzer_summary_test.xlsx test-data/Variant_Analyzer_tiers_test.xlsx test-data/tag_count_dict_test.json |
| diffstat | 13 files changed, 658 insertions(+), 333 deletions(-) [+] |
line wrap: on
line diff
--- a/mut2read.py Mon Mar 29 09:22:57 2021 +0000 +++ b/mut2read.py Fri Jul 22 09:19:44 2022 +0000 @@ -20,6 +20,7 @@ import argparse import json import os +import re import sys import numpy as np @@ -68,6 +69,7 @@ # get tags tag_dict = {} cvrg_dict = {} + tag_dict_ref = {} for variant in VCF(file1): chrom = variant.CHROM @@ -77,54 +79,108 @@ continue else: alt = variant.ALT[0] + alt = alt.upper() + ref = ref.upper() + if "N" in alt: # skip indels with N in alt allele --> it is not an indel but just a mismatch at the position where the N is (checked this in IGV) + continue chrom_stop_pos = str(chrom) + "#" + str(stop_pos) + "#" + ref + "#" + alt dcs_len = [] - if len(ref) == len(alt): - for pileupcolumn in bam.pileup(chrom, stop_pos - 1, stop_pos + 1, max_depth=100000000): - if pileupcolumn.reference_pos == stop_pos: - 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 + for pileupcolumn in bam.pileup(chrom, stop_pos - 1, stop_pos + 1, max_depth=100000000): + if pileupcolumn.reference_pos == stop_pos: + count_alt = 0 + count_ref = 0 + count_indel = 0 + count_n = 0 + count_other = 0 + count_lowq = 0 + for pileupread in pileupcolumn.pileups: + if not pileupread.is_refskip: + if pileupread.is_del: + p = pileupread.query_position_or_next + e = p + len(alt) - 1 + else: + p = pileupread.query_position + e = p + len(alt) + s = p + split_cigar = re.split('(\d+)', pileupread.alignment.cigarstring) + if len(ref) < len(alt): + if "I" in split_cigar: + all_insertions = [inser_i for inser_i, ins in enumerate(split_cigar) if ins == "I"] + for ai in all_insertions: # if multiple insertions in DCS + ins_index = [int(ci) for ci in split_cigar[:ai - 1] if ci.isdigit()] + ins_count = split_cigar[ai - 1] # nr of insertions should match with alt allele + if "I" in split_cigar and sum(ins_index) == p + 1 and len(alt) - 1 == int(ins_count): # if position in read matches and length of insertion + nuc = pileupread.alignment.query_sequence[s:e] + break + else: + nuc = pileupread.alignment.query_sequence[s] 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) + nuc = pileupread.alignment.query_sequence[s] + elif len(ref) > len(alt): + ref_positions = pileupread.alignment.get_reference_positions(full_length=True)[s:p + len(ref)] + if "D" in split_cigar: + all_deletions = [del_i for del_i, dele in enumerate(split_cigar) if dele == "D"] + for di, ai in enumerate(all_deletions): # if multiple insertions in DCS + if di > 0: # more than 1 deletion, don't count previous deletion to position + all_deletions_mod = split_cigar[:ai - 1] + prev_del_idx = [all_deletions_mod.index("D") - 1, all_deletions_mod.index("D")] + split_cigar_no_prev = [ad for i, ad in enumerate(all_deletions_mod) if i not in prev_del_idx] + del_index = [int(ci) for ci in split_cigar_no_prev[:ai - 1] if ci.isdigit()] + else: # first deletion in read, sum all previous (mis)matches and insertions to position + del_index = [int(ci) for ci in split_cigar[:ai - 1] if ci.isdigit()] + del_count = split_cigar[ai - 1] # nr of deletions should match with ref allele + if "D" in split_cigar and sum(del_index) == p + 1 and len(ref) - 1 == int(del_count): + nuc = pileupread.alignment.query_sequence[s:e] + if nuc == "": + nuc = str(alt) + break + else: + nuc = pileupread.alignment.query_sequence[s:s + len(ref)] + elif len(ref_positions) < len(ref): # DCS has reference but the position is at the very end of the DCS and therefore not the full reference positions are there + nuc = pileupread.alignment.get_reference_sequence()[s:s + len(ref)] + if nuc.upper() == ref[:len(nuc)]: + nuc = str(ref) + else: + nuc = pileupread.alignment.query_sequence[s:s + len(ref)] + else: # SNV: query position is None if is_del or is_refskip is set. + nuc = pileupread.alignment.query_sequence[s] - 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)) - else: - print("indels are currently not evaluated") + nuc = nuc.upper() + tag = pileupread.alignment.query_name + if "_" in tag: + tag = re.split('_', tag)[0] + + if nuc == alt: + count_alt += 1 + 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 + if tag in tag_dict_ref: + tag_dict_ref[tag][chrom_stop_pos] = ref + else: + tag_dict_ref[tag] = {} + tag_dict_ref[tag][chrom_stop_pos] = ref + elif nuc == "N": + count_n += 1 + elif nuc == "lowQ": + count_lowq += 1 + else: + count_other += 1 + dcs_len.append(len(pileupread.alignment.query_sequence)) + + 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) + json.dump((tag_dict, cvrg_dict, tag_dict_ref), f) # create fastq from aligned reads with open(outfile, 'w') as out: @@ -134,12 +190,11 @@ splits = line.split('\t') tag = splits[0] - if tag in tag_dict: + if tag in tag_dict or tag in tag_dict_ref: 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")
--- a/mut2read.xml Mon Mar 29 09:22:57 2021 +0000 +++ b/mut2read.xml Fri Jul 22 09:19:44 2022 +0000 @@ -1,5 +1,5 @@ <?xml version="1.0" encoding="UTF-8"?> -<tool id="mut2read" name="DCS mutations to tags/reads:" version="2.1.1" profile="19.01"> +<tool id="mut2read" name="DCS mutations to tags/reads:" version="3.0.0" profile="19.01"> <description>Extracts all tags that carry a mutation in the duplex consensus sequence (DCS)</description> <macros> <import>va_macros.xml</import> @@ -17,7 +17,7 @@ ]]> </command> <inputs> - <param name="file1" type="data" format="vcf" label="DCS Mutation File" optional="false" help="VCF file with DCS mutations. See Help section below for a detailed explanation."/> + <param name="file1" type="data" format="vcf" label="DCS Mutation File" optional="false" help="VCF file with DCS mutations. See the Help section below for a detailed explanation."/> <param name="file2" type="data" format="bam" label="DCS BAM File" optional="false" help="BAM file with aligned DCS reads."/> <param name="file3" type="data" format="tabular" label="Aligned Families File" optional="false" help="TABULAR file with aligned families."/> </inputs> @@ -27,11 +27,11 @@ </outputs> <tests> <test> - <param name="file1" value="FreeBayes_test.vcf" lines_diff="2"/> + <param name="file1" value="FreeBayes_test.vcf"/> <param name="file2" value="DCS_test.bam"/> <param name="file3" value="Aligned_Families_test.tabular"/> - <output name="output_fastq" file="Interesting_Reads_test.fastq" lines_diff="136"/> - <output name="output_json" file="tag_count_dict_test.json" lines_diff="2"/> + <output name="output_fastq" file="Interesting_Reads_test.fastq"/> + <output name="output_json" file="tag_count_dict_test.json"/> </test> </tests> <help> <![CDATA[ @@ -39,12 +39,12 @@ Takes a VCF file with mutations, a BAM file of aligned DCS reads, and a tabular file with aligned families as input and prints all tags of reads that -carry a mutation to a user specified output file and creates a fastq file of -reads of tags with a mutation. +carry a mutation or have the reference allele to a user-specified output file and creates a fastq file of +reads of tags with a mutation and the reference allele. **Input** -**Dataset 1:** VCF file with duplex consesus sequence (DCS) mutations. E.g. +**Dataset 1:** VCF file with duplex consensus sequence (DCS) mutations. E.g. generated by the `FreeBayes <https://arxiv.org/abs/1207.3907>`_ or `LoFreq <https://academic.oup.com/nar/article/40/22/11189/1152727>`_ variant caller. @@ -57,7 +57,7 @@ **Output** -The output is a json file containing dictonaries of the tags of reads containing mutations +The output is a json file containing dictionaries of the tags of reads containing mutations in the DCS and a fastq file of all reads of these tags. ]]>
--- a/mut2sscs.py Mon Mar 29 09:22:57 2021 +0000 +++ b/mut2sscs.py Fri Jul 22 09:19:44 2022 +0000 @@ -23,6 +23,7 @@ import argparse import json import os +import re import sys import numpy as np @@ -71,60 +72,111 @@ continue else: alt = variant.ALT[0] + alt = alt.upper() + ref = ref.upper() + if "N" in alt: # skip indels with N in alt allele --> it is not an indel but just a mismatch at the position where the N is (checked this in IGV) + continue chrom_stop_pos = str(chrom) + "#" + str(stop_pos) + "#" + ref + "#" + alt - if len(ref) == len(alt): - for pileupcolumn in bam.pileup(chrom, stop_pos - 1, stop_pos + 1, max_depth=1000000000): - if pileupcolumn.reference_pos == stop_pos: - 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 + for pileupcolumn in bam.pileup(chrom, stop_pos - 1, stop_pos + 1, max_depth=1000000000): + if pileupcolumn.reference_pos == stop_pos: + count_alt = 0 + count_ref = 0 + count_indel = 0 + for pileupread in pileupcolumn.pileups: + if not pileupread.is_refskip: + tag = pileupread.alignment.query_name + abba = tag[-2:] + if pileupread.is_del: + p = pileupread.query_position_or_next + e = p + len(alt) - 1 + else: + p = pileupread.query_position + e = p + len(alt) + tag = pileupread.alignment.query_name + if "_" in tag: + tag = re.split('_', tag)[0] + s = p + split_cigar = re.split('(\d+)', pileupread.alignment.cigarstring) + if len(ref) < len(alt): # insertion + if "I" in split_cigar: + all_insertions = [inser_i for inser_i, ins in enumerate(split_cigar) if ins == "I"] + for ai in all_insertions: # if multiple insertions in DCS + ins_index = [int(ci) for ci in split_cigar[:ai - 1] if ci.isdigit()] + ins_count = split_cigar[ai - 1] # nr of insertions should match with alt allele + if "I" in split_cigar and sum(ins_index) == p + 1 and len(alt) - 1 == int(ins_count): + nuc = pileupread.alignment.query_sequence[s:e] + break else: - mut_pos_dict[chrom_stop_pos][abba] = 1 - else: - mut_pos_dict[chrom_stop_pos] = {} - mut_pos_dict[chrom_stop_pos][abba] = 1 - if chrom_stop_pos not in ref_pos_dict: - ref_pos_dict[chrom_stop_pos] = {} - ref_pos_dict[chrom_stop_pos][abba] = 0 + nuc = pileupread.alignment.query_sequence[s] + else: + nuc = pileupread.alignment.query_sequence[s] + elif len(ref) > len(alt): # deletion + ref_positions = pileupread.alignment.get_reference_positions(full_length=True)[s:p + len(ref)] + if "D" in split_cigar: + all_deletions = [del_i for del_i, dele in enumerate(split_cigar) if dele == "D"] + for di, ai in enumerate(all_deletions): # if multiple insertions in DCS + if di > 0: # more than 1 deletion, don't count previous deletion to position + all_deletions_mod = split_cigar[:ai - 1] + prev_del_idx = [all_deletions_mod.index("D") - 1, all_deletions_mod.index("D")] + split_cigar_no_prev = [ad for i, ad in enumerate(all_deletions_mod) if i not in prev_del_idx] + del_index = [int(ci) for ci in split_cigar_no_prev[:ai - 1] if ci.isdigit()] + else: # first deletion in read, sum all previous (mis)matches and insertions to position + del_index = [int(ci) for ci in split_cigar[:ai - 1] if ci.isdigit()] + del_count = split_cigar[ai - 1] # nr of deletions should match with ref allele + if "D" in split_cigar and sum(del_index) == p + 1 and len(ref) - 1 == int(del_count): + nuc = pileupread.alignment.query_sequence[s:e] + if nuc == "": + nuc = str(alt) + break + else: + nuc = pileupread.alignment.query_sequence[s:s + len(ref)] + elif len(ref_positions) < len(ref): # DCS has reference but the position is at the very end of the DCS and therefore not the full reference positions are there + nuc = pileupread.alignment.get_reference_sequence()[s:s + len(ref)] + if nuc.upper() == ref[:len(nuc)]: + nuc = str(ref) + else: + nuc = pileupread.alignment.query_sequence[s:s + len(ref)] + else: # SNV: query position is None if is_del or is_refskip is set. + nuc = pileupread.alignment.query_sequence[s] - 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 + nuc = nuc.upper() + + if nuc == 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: - ref_pos_dict[chrom_stop_pos] = {} + mut_pos_dict[chrom_stop_pos][abba] = 1 + else: + mut_pos_dict[chrom_stop_pos] = {} + mut_pos_dict[chrom_stop_pos][abba] = 1 + if chrom_stop_pos not in ref_pos_dict: + ref_pos_dict[chrom_stop_pos] = {} + ref_pos_dict[chrom_stop_pos][abba] = 0 + elif nuc == 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: - 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)) + 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, other = %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 - else: - print("indels are currently not evaluated") + # 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
--- a/mut2sscs.xml Mon Mar 29 09:22:57 2021 +0000 +++ b/mut2sscs.xml Fri Jul 22 09:19:44 2022 +0000 @@ -1,6 +1,6 @@ <?xml version="1.0" encoding="UTF-8"?> -<tool id="mut2sscs" name="DCS mutations to SSCS stats:" version="2.1.1" profile="19.01"> - <description>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</description> +<tool id="mut2sscs" name="DCS mutations to SSCS stats:" version="3.0.0" profile="19.01"> + <description>Extracts all tags from the single-stranded consensus sequence (SSCS) bam file that carries a mutation at the same position a mutation is called in the duplex consensus sequence (DCS) and calculates their frequencies</description> <macros> <import>va_macros.xml</import> </macros> @@ -15,7 +15,7 @@ ]]> </command> <inputs> - <param name="file1" type="data" format="vcf" label="DCS Mutation File" optional="false" help="VCF file with DCS mutations. See Help section below for a detailed explanation."/> + <param name="file1" type="data" format="vcf" label="DCS Mutation File" optional="false" help="VCF file with DCS mutations. See the Help section below for a detailed explanation."/> <param name="file2" type="data" format="bam" label="SSCS BAM File" optional="false" help="BAM file with aligned SSCS reads."/> </inputs> <outputs> @@ -25,29 +25,29 @@ <test> <param name="file1" value="FreeBayes_test.vcf"/> <param name="file2" value="SSCS_test.bam"/> - <output name="output_json" file="SSCS_counts_test.json" lines_diff="2"/> + <output name="output_json" file="SSCS_counts_test.json"/> </test> </tests> <help> <![CDATA[ **What it does** Takes a VCF file with DCS mutations and a BAM file of aligned SSCS reads -as input and writes statistics about tags of reads that carry a mutation in the -SSCS at the same position a mutation is called in the DCS to a user specified output file.. +as input and writes statistics about tags of reads that carry a mutation and the reference allele in the +SSCS at the same position a mutation is called in the DCS to a user-specified output file. **Input** -**Dataset 1:** VCF file with duplex consesus sequence (DCS) mutations. E.g. +**Dataset 1:** VCF file with duplex consensus sequence (DCS) mutations. E.g. generated by the `FreeBayes <https://arxiv.org/abs/1207.3907>`_ or `LoFreq <https://academic.oup.com/nar/article/40/22/11189/1152727>`_ variant caller. -**Dataset 2:** BAM file of aligned single stranded consensus sequence (SSCS) +**Dataset 2:** BAM file of the aligned single stranded consensus sequence (SSCS) reads. This file can be obtained by the tool `Map with BWA-MEM <https://arxiv.org/abs/1303.3997>`_. **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. +The output is a json file containing dictionaries with stats of tags that carry a mutation in the SSCS +at the same position, a mutation is called in the DCS. ]]> </help>
--- a/read2mut.py Mon Mar 29 09:22:57 2021 +0000 +++ b/read2mut.py Fri Jul 22 09:19:44 2022 +0000 @@ -87,6 +87,7 @@ outfile2 = args.outputFile2 outfile3 = args.outputFile3 outputFile_csv = args.outputFile_csv + thresh = args.thresh phred_score = args.phred trim = args.trim @@ -111,7 +112,7 @@ # load dicts with open(json_file, "r") as f: - (tag_dict, cvrg_dict) = json.load(f) + (tag_dict, cvrg_dict, tag_dict_ref) = json.load(f) with open(sscs_json, "r") as f: (mut_pos_dict, ref_pos_dict) = json.load(f) @@ -138,106 +139,207 @@ continue else: alt = variant.ALT[0] - chrom_stop_pos = str(chrom) + "#" + str(stop_pos) + "#" + ref + "#" + alt + + alt = alt.upper() + ref = ref.upper() + if "N" in alt: # skip indels with N in alt allele --> it is not an indel but just a mismatch at the position where the N is (checked this in IGV) + continue - if len(ref) == len(alt): - mut_array.append([chrom, stop_pos, ref, alt]) - i += 1 - mut_dict[chrom_stop_pos] = {} - mut_read_pos_dict[chrom_stop_pos] = {} - reads_dict[chrom_stop_pos] = {} - mut_read_cigar_dict[chrom_stop_pos] = {} - real_start_end[chrom_stop_pos] = {} + chrom_stop_pos = str(chrom) + "#" + str(stop_pos) + "#" + ref + "#" + alt + mut_array.append([chrom, stop_pos, ref, alt]) + i += 1 + mut_dict[chrom_stop_pos] = {} + mut_read_pos_dict[chrom_stop_pos] = {} + reads_dict[chrom_stop_pos] = {} + mut_read_cigar_dict[chrom_stop_pos] = {} + real_start_end[chrom_stop_pos] = {} + for pileupcolumn in bam.pileup(chrom, stop_pos - 1, stop_pos + 1, max_depth=1000000000): + if pileupcolumn.reference_pos == stop_pos: + count_alt = 0 + count_ref = 0 + count_indel = 0 + count_n = 0 + count_other = 0 + count_lowq = 0 + n = 0 + for pileupread in pileupcolumn.pileups: + n += 1 + if not pileupread.is_refskip: + if pileupread.is_del: + p = pileupread.query_position_or_next + e = p + len(alt) - 1 + else: + p = pileupread.query_position + e = p + len(alt) + s = p + tag = pileupread.alignment.query_name + split_cigar = re.split('(\d+)', pileupread.alignment.cigarstring) - for pileupcolumn in bam.pileup(chrom, stop_pos - 1, stop_pos + 1, max_depth=100000000): - if pileupcolumn.reference_pos == stop_pos: - count_alt = 0 - count_ref = 0 - count_indel = 0 - count_n = 0 - count_other = 0 - count_lowq = 0 - n = 0 - 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 read is softclipped, store real position in reference - if "S" in pileupread.alignment.cigarstring: - # spftclipped at start - if re.search(r"^[0-9]+S", pileupread.alignment.cigarstring): - start = pileupread.alignment.reference_start - int(pileupread.alignment.cigarstring.split("S")[0]) - end = pileupread.alignment.reference_end - # softclipped at end - elif re.search(r"S$", pileupread.alignment.cigarstring): - end = pileupread.alignment.reference_end + int(re.split("[A-Z]", str(pileupread.alignment.cigarstring))[-2]) - start = pileupread.alignment.reference_start - else: + if len(ref) < len(alt): + if "I" in split_cigar: + all_insertions = [inser_i for inser_i, ins in enumerate(split_cigar) if ins == "I"] + for ai in all_insertions: # if multiple insertions in DCS + ins_index = [int(ci) for ci in split_cigar[:ai - 1] if ci.isdigit()] + ins_count = split_cigar[ai - 1] # nr of insertions should match with alt allele + if "I" in split_cigar and sum(ins_index) == p + 1 and int(ins_count) >= len(alt) - 1: # if pe read matches exatcly to insertion + nuc = pileupread.alignment.query_sequence[s:e] + phred = ord(pileupread.alignment.qual[s]) - 33 + break + elif "I" in split_cigar and sum(ins_index) == p + 1 and int(ins_count) < len(alt) - 1: # if pe read has shorter insertion -- not alt + nuc = pileupread.alignment.query_sequence[s:s+int(ins_count)+1] + phred = ord(pileupread.alignment.qual[s]) - 33 + break + else: # insertion in pe reads but not at the desired position + nuc = pileupread.alignment.query_sequence[s] + phred = ord(pileupread.alignment.qual[s]) - 33 + elif "D" in split_cigar: # if deletion in pe read, don't count + all_deletions = [del_i for del_i, dele in enumerate(split_cigar) if dele == "D"] + for di, ai in enumerate(all_deletions): # if multiple insertions in DCS + if di > 0: # more than 1 deletion, don't count previous deletion to position + all_deletions_mod = split_cigar[:ai - 1] + prev_del_idx = [all_deletions_mod.index("D") - 1, all_deletions_mod.index("D")] + split_cigar_no_prev = [ad for i, ad in enumerate(all_deletions_mod) if i not in prev_del_idx] + del_index = [int(ci) for ci in split_cigar_no_prev[:ai - 1] if ci.isdigit()] + else: # first deletion in read, sum all previous (mis)matches and insertions to position + del_index = [int(ci) for ci in split_cigar[:ai - 1] if ci.isdigit()] + if "D" in split_cigar and sum(del_index) == p + 1: # if deletion at that position + nuc = "D" + phred = ord(pileupread.alignment.qual[s]) - 33 + break + else: + nuc = pileupread.alignment.query_sequence[s] + phred = ord(pileupread.alignment.qual[s]) - 33 + else: # insertion in pe reads but not at the desired position + nuc = pileupread.alignment.query_sequence[s] + phred = ord(pileupread.alignment.qual[s]) - 33 + + elif len(ref) > len(alt): + ref_positions = pileupread.alignment.get_reference_positions(full_length=True)[s:p + len(ref)] + if "D" in split_cigar: + all_deletions = [del_i for del_i, dele in enumerate(split_cigar) if dele == "D"] + for di, ai in enumerate(all_deletions): # if multiple insertions in DCS + if di > 0: # more than 1 deletion, don't count previous deletion to position + all_deletions_mod = split_cigar[:ai - 1] + prev_del_idx = [all_deletions_mod.index("D") - 1, all_deletions_mod.index("D")] + split_cigar_no_prev = [ad for i, ad in enumerate(all_deletions_mod) if i not in prev_del_idx] + del_index = [int(ci) for ci in split_cigar_no_prev[:ai - 1] if ci.isdigit()] + else: # first deletion in read, sum all previous (mis)matches and insertions to position + del_index = [int(ci) for ci in split_cigar[:ai - 1] if ci.isdigit()] + del_count = split_cigar[ai - 1] # deletion on that position but does not necesserily match in the length + if "D" in split_cigar and sum(del_index) == p + 1 and int(del_count) >= len(ref) - 1: # if pe read matches exatcly to deletion + nuc = pileupread.alignment.query_sequence[s:e] + phred = ord(pileupread.alignment.qual[s]) - 33 + if nuc == "": + nuc = str(alt) + break + elif "D" in split_cigar and sum(del_index) == p + 1 and int(del_count) < len(ref) - 1: # if pe read has shorter deletion --> not alt + nuc = str(ref)[:int(del_count)+1] + phred = ord(pileupread.alignment.qual[s]) - 33 + break + else: # deletion in pe reads but not at the desired position + nuc = pileupread.alignment.query_sequence[s:s + len(ref)] + phred = ord(pileupread.alignment.qual[s]) - 33 + elif "I" in split_cigar: # if pe read has insertion --> not count + all_insertions = [inser_i for inser_i, ins in enumerate(split_cigar) if ins == "I"] + for ai in all_insertions: # if multiple insertions in DCS + ins_index = [int(ci) for ci in split_cigar[:ai - 1] if ci.isdigit()] + ins_count = split_cigar[ai - 1] # nr of insertions should match with alt allele + if "I" in split_cigar and sum(ins_index) == p + 1: + nuc = "I" + phred = ord(pileupread.alignment.qual[s]) - 33 + break + else: + nuc = pileupread.alignment.query_sequence[s] + phred = ord(pileupread.alignment.qual[s]) - 33 + elif len(ref_positions) < len(ref): # DCS has reference but the position is at the very end of the DCS and therefore not the full reference positions are there + nuc = pileupread.alignment.get_reference_sequence()[s:s + len(ref)] + phred = ord(pileupread.alignment.qual[s]) - 33 + if nuc.upper() == ref[:len(nuc)]: + nuc = str(ref) + else: # deletion in pe reads but not at the desired position + nuc = pileupread.alignment.query_sequence[s:s + len(ref)] + phred = ord(pileupread.alignment.qual[s]) - 33 + else: # SNV: query position is None if is_del or is_refskip is set. + nuc = pileupread.alignment.query_sequence[s] + phred = ord(pileupread.alignment.qual[s]) - 33 + + nuc = nuc.upper() + + # if read is softclipped, store real position in reference + if "S" in pileupread.alignment.cigarstring: + # spftclipped at start + if re.search(r"^[0-9]+S", pileupread.alignment.cigarstring): + start = pileupread.alignment.reference_start - int(pileupread.alignment.cigarstring.split("S")[0]) end = pileupread.alignment.reference_end + # softclipped at end + elif re.search(r"S$", pileupread.alignment.cigarstring): + end = pileupread.alignment.reference_end + int(re.split("[A-Z]", str(pileupread.alignment.cigarstring))[-2]) start = pileupread.alignment.reference_start - - 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] = [pileupread.query_position + 1] - reads_dict[chrom_stop_pos][tag] = [len(pileupread.alignment.query_sequence)] - mut_read_cigar_dict[chrom_stop_pos][tag] = [pileupread.alignment.cigarstring] - real_start_end[chrom_stop_pos][tag] = [(start, end)] + else: + end = pileupread.alignment.reference_end + start = pileupread.alignment.reference_start + 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] = [s + 1] + reads_dict[chrom_stop_pos][tag] = [len(pileupread.alignment.query_sequence)] + mut_read_cigar_dict[chrom_stop_pos][tag] = [pileupread.alignment.cigarstring] + real_start_end[chrom_stop_pos][tag] = [(start, end)] + else: + mut_read_pos_dict[chrom_stop_pos][tag].append(s + 1) + reads_dict[chrom_stop_pos][tag].append(len(pileupread.alignment.query_sequence)) + mut_read_cigar_dict[chrom_stop_pos][tag].append(pileupread.alignment.cigarstring) + real_start_end[chrom_stop_pos][tag].append((start, end)) + 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, ref) else: - mut_read_pos_dict[chrom_stop_pos][tag].append(pileupread.query_position + 1) - reads_dict[chrom_stop_pos][tag].append(len(pileupread.alignment.query_sequence)) - mut_read_cigar_dict[chrom_stop_pos][tag].append(pileupread.alignment.cigarstring) - real_start_end[chrom_stop_pos][tag].append((start, end)) - 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, ref) - else: - mut_read_dict[tag][chrom_stop_pos] = (alt, ref) - elif nuc == ref: - count_ref += 1 - elif nuc == "N": - count_n += 1 - elif nuc == "lowQ": - count_lowq += 1 - else: - count_other += 1 + mut_read_dict[tag][chrom_stop_pos] = (alt, ref) + elif nuc == ref: + count_ref += 1 + elif nuc == "N": + count_n += 1 + elif nuc == "lowQ": + count_lowq += 1 else: - count_indel += 1 + count_other += 1 + else: + count_indel += 1 mut_array = np.array(mut_array) 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 + if pure_tag in tag_dict.keys(): # stored all ref and alt reads --> get only alt reads + 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() - # create pure_tags_dict pure_tags_dict = {} + pure_tags_dict_ref = {} 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[:, 0], mut_array[:, 1], mut_array[:, 2], mut_array[:, 3])]) == key1)[0][0] ref = mut_array[i, 2] alt = mut_array[i, 3] pure_tags_dict[key1] = {} + pure_tags_dict_ref[key1] = {} for key2, value2 in sorted(value1.items()): for key3, value3 in value2.items(): pure_tag = key2[:-5] @@ -247,6 +349,12 @@ else: pure_tags_dict[key1][pure_tag] = 1 + if key3 == ref: + if pure_tag in pure_tags_dict_ref[key1]: + pure_tags_dict_ref[key1][pure_tag] += 1 + else: + pure_tags_dict_ref[key1][pure_tag] = 1 + # create pure_tags_dict_short with thresh if thresh > 0: pure_tags_dict_short = {} @@ -279,7 +387,7 @@ format23 = workbook3.add_format({'bg_color': '#FFC7CE'}) # red format33 = workbook3.add_format({'bg_color': '#FACC2E'}) # yellow - header_line = ('variant ID', 'tier', 'tag', 'mate', 'read pos.ab', 'read pos.ba', 'read median length.ab', + header_line = ('variant ID', 'tier', 'allele', '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', @@ -288,6 +396,7 @@ 'in phase', 'chimeric tag') ws1.write_row(0, 0, header_line) csv_writer.writerow(header_line) + counter_tier11 = 0 counter_tier12 = 0 counter_tier21 = 0 @@ -308,12 +417,14 @@ row = 1 tier_dict = {} + tier_dict_ref = {} chimera_dict = {} for key1, value1 in sorted(mut_dict.items()): counts_mut = 0 chimeric_tag_list = [] chimeric_tag = {} - if key1 in pure_tags_dict_short.keys(): + if (key1 in pure_tags_dict_short.keys()) or (key1 in pure_tags_dict_ref.keys()): # ref or alt + change_tier_after_print = [] i = np.where(np.array(['#'.join(str(i) for i in z) for z in zip(mut_array[:, 0], mut_array[:, 1], mut_array[:, 2], mut_array[:, 3])]) == key1)[0][0] @@ -323,11 +434,13 @@ whole_array = list(pure_tags_dict_short[key1].keys()) tier_dict[key1] = {} + tier_dict_ref[key1] = {} 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 2.5", 0), ("tier 3.1", 0), ("tier 3.2", 0), ("tier 4", 0), ("tier 5.1", 0), ("tier 5.2", 0), ("tier 5.3", 0), ("tier 5.4", 0), ("tier 5.5", 0), ("tier 6", 0), ("tier 7", 0)] for k, v in values_tier_dict: tier_dict[key1][k] = v + tier_dict_ref[key1][k] = v used_keys = [] if 'ab' in mut_pos_dict[key1].keys(): @@ -349,7 +462,16 @@ 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] not in tag_dict.keys() and key2[:-5] not in tag_dict_ref.keys(): # skip reads that have not alt or ref + continue + + if ((key2[:-5] in tag_dict.keys() and key2[:-5] in pure_tags_dict_short[key1].keys() and key1 in tag_dict[key2[:-5]].keys()) or (key2[:-5] in tag_dict_ref.keys() and key2[:-5] in pure_tags_dict_ref[key1].keys() and key1 in tag_dict_ref[key2[:-5]].keys())) and (key2[:-5] not in used_keys): + + if key2[:-5] in pure_tags_dict_short[key1].keys(): + variant_type = "alt" + elif key2[:-5] in pure_tags_dict_ref[key1].keys(): + variant_type = "ref" + 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(): @@ -550,36 +672,59 @@ used_keys.append(key2[:-5]) counts_mut += 1 - if (alt1f + alt2f + alt3f + alt4f) > 0.5: + if (variant_type == "alt" and ((alt1f + alt2f + alt3f + alt4f) > 0.5)) or (variant_type == "ref" and ((ref1f + ref2f + ref3f + ref4f) > 0.5)): + if variant_type == "alt": + tier1ff, tier2ff, tier3ff, tier4ff = alt1f, alt2f, alt3f, alt4f + tier1ff_trim, tier2ff_trim, tier3ff_trim, tier4ff_trim = alt1f, alt2f, alt3f, alt4f + elif variant_type == "ref": + tier1ff, tier2ff, tier3ff, tier4ff = ref1f, ref2f, ref3f, ref4f + tier1ff_trim, tier2ff_trim, tier3ff_trim, tier4ff_trim = ref1f, ref2f, ref3f, ref4f + total1new_trim, total2new_trim, total3new_trim, total4new_trim = total1new, total2new, total3new, total4new if total1new == 0: ref1f = alt1f = None alt1ff = -1 alt1ff_trim = -1 + tier1ff = -1 + tier1ff_trim = -1 else: alt1ff = alt1f alt1ff_trim = alt1f + tier1ff = tier1ff + tier1ff_trim = tier1ff_trim if total2new == 0: ref2f = alt2f = None alt2ff = -1 alt2ff_trim = -1 + tier2ff = -1 + tier2ff_trim = -1 else: alt2ff = alt2f alt2ff_trim = alt2f + tier2ff = tier2ff + tier2ff_trim = tier2ff_trim if total3new == 0: ref3f = alt3f = None alt3ff = -1 alt3ff_trim = -1 + tier3ff = -1 + tier3ff_trim = -1 else: alt3ff = alt3f alt3ff_trim = alt3f + tier3ff = tier3ff + tier3ff_trim = tier3ff_trim if total4new == 0: ref4f = alt4f = None alt4ff = -1 alt4ff_trim = -1 + tier4ff = -1 + tier4ff_trim = -1 else: alt4ff = alt4f alt4ff_trim = alt4f + tier4ff = tier4ff + tier4ff_trim = tier4ff_trim beg1 = beg4 = beg2 = beg3 = 0 @@ -605,7 +750,7 @@ # mate 1 - SSCS ab softclipped_idx1 = [True if re.search(r"^[0-9]+S", string) or re.search(r"S$", string) else False for string in cigars_dcs1] safe_div_result = safe_div(sum(softclipped_idx1), float(len(softclipped_idx1))) - if (safe_div_result == None): + if (safe_div_result is None): ratio1 = False else: ratio1 = safe_div_result >= threshold_reads @@ -629,7 +774,7 @@ # mate 1 - SSCS ba softclipped_idx4 = [True if re.search(r"^[0-9]+S", string) or re.search(r"S$", string) else False for string in cigars_dcs4] safe_div_result = safe_div(sum(softclipped_idx4), float(len(softclipped_idx4))) - if (safe_div_result == None): + if (safe_div_result is None): ratio4 = False else: ratio4 = safe_div_result >= threshold_reads @@ -653,7 +798,7 @@ # mate 2 - SSCS ab softclipped_idx2 = [True if re.search(r"^[0-9]+S", string) or re.search(r"S$", string) else False for string in cigars_dcs2] safe_div_result = safe_div(sum(softclipped_idx2), float(len(softclipped_idx2))) - if (safe_div_result == None): + if (safe_div_result is None): ratio2 = False else: ratio2 = safe_div_result >= threshold_reads @@ -677,7 +822,7 @@ # mate 2 - SSCS ba softclipped_idx3 = [True if re.search(r"^[0-9]+S", string) or re.search(r"S$", string) else False for string in cigars_dcs3] safe_div_result = safe_div(sum(softclipped_idx3), float(len(softclipped_idx3))) - if (safe_div_result == None): + if (safe_div_result is None): ratio3 = False else: ratio3 = safe_div_result >= threshold_reads @@ -698,21 +843,21 @@ ratio_dist_start3 = safe_div(sum([True if x <= thr else False for x in dist_start_read3]), float(sum(softclipped_idx3))) >= threshold_reads ratio_dist_end3 = safe_div(sum([True if x <= thr else False for x in dist_end_read3]), float(sum(softclipped_idx3))) >= threshold_reads - if ((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 + if ((all(float(ij) >= 0.5 for ij in [tier1ff, tier4ff]) & # contradictory variant + all(float(ij) == 0. for ij in [tier2ff, tier3ff])) | + (all(float(ij) >= 0.5 for ij in [tier2ff, tier3ff]) & + all(float(ij) == 0. for ij in [tier1ff, tier4ff]))): + tier1ff = 0 + tier4ff = 0 + tier2ff = 0 + tier3ff = 0 trimmed = False contradictory = True # softclipping tiers # information of both mates available --> all reads for both mates and SSCS are softclipped elif (ratio1 & ratio4 & ratio2 & ratio3 & (ratio_dist_start1 | ratio_dist_end1) & (ratio_dist_start4 | ratio_dist_end4) & (ratio_dist_start2 | ratio_dist_end2) & (ratio_dist_start3 | ratio_dist_end3) & - all(float(ij) > 0. for ij in [alt1ff, alt2ff, alt3ff, alt4ff])): # all mates available + all(float(ij) > 0. for ij in [tier1ff, tier2ff, tier3ff, tier4ff])): # all mates available # if distance between softclipping and mutation is at start or end of the read smaller than threshold softclipped_mutation_allMates = True softclipped_mutation_oneOfTwoMates = False @@ -720,16 +865,16 @@ softclipped_mutation_oneOfTwoSSCS_diffMates = False softclipped_mutation_oneMate = False softclipped_mutation_oneMateOneSSCS = False - alt1ff = 0 - alt4ff = 0 - alt2ff = 0 - alt3ff = 0 + tier1ff = 0 + tier4ff = 0 + tier2ff = 0 + tier3ff = 0 trimmed = False contradictory = False # information of both mates available --> only one mate softclipped elif (((ratio1 & ratio4 & (ratio_dist_start1 | ratio_dist_end1) & (ratio_dist_start4 | ratio_dist_end4)) | (ratio2 & ratio3 & (ratio_dist_start2 | ratio_dist_end2) & (ratio_dist_start3 | ratio_dist_end3))) & - all(float(ij) > 0. for ij in [alt1ff, alt2ff, alt3ff, alt4ff])): # all mates available + all(float(ij) > 0. for ij in [tier1ff, tier2ff, tier3ff, tier4ff])): # all mates available # if distance between softclipping and mutation is at start or end of the read smaller than threshold min_start1 = min(min([ij[0] for ij in ref_positions1]), min([ij[0] for ij in ref_positions4])) # red min_start2 = min(min([ij[0] for ij in ref_positions2]), min([ij[0] for ij in ref_positions3])) # blue @@ -763,10 +908,10 @@ read_len_median3 = read_len_median3 - n_spacer_barcode else: softclipped_mutation_oneOfTwoMates = True - alt1ff = 0 - alt4ff = 0 - alt2ff = 0 - alt3ff = 0 + tier1ff = 0 + tier4ff = 0 + tier2ff = 0 + tier3ff = 0 trimmed = False contradictory = False softclipped_mutation_allMates = False @@ -780,6 +925,7 @@ total1new = 0 alt1ff = 0 alt1f = 0 + tier1ff = 0 trimmed = True if ((read_pos4 >= 0) and ((read_pos4 <= trim) | (abs(read_len_median4 - read_pos4) <= trim))): @@ -787,6 +933,7 @@ total4new = 0 alt4ff = 0 alt4f = 0 + tier4ff = 0 trimmed = True if ((read_pos2 >= 0) and ((read_pos2 <= trim) | (abs(read_len_median2 - read_pos2) <= trim))): @@ -794,6 +941,7 @@ total2new = 0 alt2ff = 0 alt2f = 0 + tier2ff = 0 trimmed = True if ((read_pos3 >= 0) and ((read_pos3 <= trim) | (abs(read_len_median3 - read_pos3) <= trim))): @@ -801,6 +949,7 @@ total3new = 0 alt3ff = 0 alt3f = 0 + tier3ff = 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) @@ -808,7 +957,7 @@ # information of both mates available --> only one mate softclipped elif (((ratio1 & (ratio_dist_start1 | ratio_dist_end1)) | (ratio4 & (ratio_dist_start4 | ratio_dist_end4))) & ((ratio2 & (ratio_dist_start2 | ratio_dist_end2)) | (ratio3 & (ratio_dist_start3 | ratio_dist_end3))) & - all(float(ij) > 0. for ij in [alt1ff, alt2ff, alt3ff, alt4ff])): # all mates available + all(float(ij) > 0. for ij in [tier1ff, tier2ff, tier3ff, tier4ff])): # all mates available # if distance between softclipping and mutation is at start or end of the read smaller than threshold softclipped_mutation_allMates = False softclipped_mutation_oneOfTwoMates = False @@ -816,18 +965,18 @@ softclipped_mutation_oneOfTwoSSCS_diffMates = False softclipped_mutation_oneMate = False softclipped_mutation_oneMateOneSSCS = False - alt1ff = 0 - alt4ff = 0 - alt2ff = 0 - alt3ff = 0 + tier1ff = 0 + tier4ff = 0 + tier2ff = 0 + tier3ff = 0 trimmed = False contradictory = False # information of one mate available --> all reads of one mate are softclipped elif ((ratio1 & ratio4 & (ratio_dist_start1 | ratio_dist_end1) & (ratio_dist_start4 | ratio_dist_end4) & - all(float(ij) < 0. for ij in [alt2ff, alt3ff]) & all(float(ij) > 0. for ij in [alt1ff, alt4ff])) | + all(float(ij) < 0. for ij in [tier2ff, tier3ff]) & all(float(ij) > 0. for ij in [tier1ff, tier4ff])) | (ratio2 & ratio3 & (ratio_dist_start2 | ratio_dist_end2) & (ratio_dist_start3 | ratio_dist_end3) & - all(float(ij) < 0. for ij in [alt1ff, alt4ff]) & all(float(ij) > 0. for ij in [alt2ff, alt3ff]))): # all mates available + all(float(ij) < 0. for ij in [tier1ff, tier4ff]) & all(float(ij) > 0. for ij in [tier2ff, tier3ff]))): # all mates available # if distance between softclipping and mutation is at start or end of the read smaller than threshold softclipped_mutation_allMates = False softclipped_mutation_oneOfTwoMates = False @@ -835,17 +984,17 @@ softclipped_mutation_oneOfTwoSSCS_diffMates = False softclipped_mutation_oneMate = True softclipped_mutation_oneMateOneSSCS = False - alt1ff = 0 - alt4ff = 0 - alt2ff = 0 - alt3ff = 0 + tier1ff = 0 + tier4ff = 0 + tier2ff = 0 + tier3ff = 0 trimmed = False contradictory = False # information of one mate available --> only one SSCS is softclipped elif ((((ratio1 & (ratio_dist_start1 | ratio_dist_end1)) | (ratio4 & (ratio_dist_start4 | ratio_dist_end4))) & - (all(float(ij) < 0. for ij in [alt2ff, alt3ff]) & all(float(ij) > 0. for ij in [alt1ff, alt4ff]))) | + (all(float(ij) < 0. for ij in [tier2ff, tier3ff]) & all(float(ij) > 0. for ij in [tier1ff, tier4ff]))) | (((ratio2 & (ratio_dist_start2 | ratio_dist_end2)) | (ratio3 & (ratio_dist_start3 | ratio_dist_end3))) & - (all(float(ij) < 0. for ij in [alt1ff, alt4ff]) & all(float(ij) < 0. for ij in [alt2ff, alt3ff])))): # all mates available + (all(float(ij) < 0. for ij in [tier1ff, tier4ff]) & all(float(ij) < 0. for ij in [tier2ff, tier3ff])))): # all mates available # if distance between softclipping and mutation is at start or end of the read smaller than threshold softclipped_mutation_allMates = False softclipped_mutation_oneOfTwoMates = False @@ -853,10 +1002,10 @@ softclipped_mutation_oneOfTwoSSCS_diffMates = False softclipped_mutation_oneMate = False softclipped_mutation_oneMateOneSSCS = True - alt1ff = 0 - alt4ff = 0 - alt2ff = 0 - alt3ff = 0 + tier1ff = 0 + tier4ff = 0 + tier2ff = 0 + tier3ff = 0 trimmed = False contradictory = False @@ -866,6 +1015,7 @@ total1new = 0 alt1ff = 0 alt1f = 0 + tier1ff = 0 trimmed = True if ((read_pos4 >= 0) and ((read_pos4 <= trim) | (abs(read_len_median4 - read_pos4) <= trim))): @@ -873,6 +1023,7 @@ total4new = 0 alt4ff = 0 alt4f = 0 + tier4ff = 0 trimmed = True if ((read_pos2 >= 0) and ((read_pos2 <= trim) | (abs(read_len_median2 - read_pos2) <= trim))): @@ -880,6 +1031,7 @@ total2new = 0 alt2ff = 0 alt2f = 0 + tier2ff = 0 trimmed = True if ((read_pos3 >= 0) and ((read_pos3 <= trim) | (abs(read_len_median3 - read_pos3) <= trim))): @@ -887,117 +1039,145 @@ total3new = 0 alt3ff = 0 alt3f = 0 + tier3ff = 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(float(ij) >= 0.75 for ij in [tier1ff, tier4ff])) | (all(int(ij) >= 3 for ij in [total2new, total3new]) & - all(float(ij) >= 0.75 for ij in [alt2ff, alt3ff]))): + all(float(ij) >= 0.75 for ij in [tier2ff, tier3ff]))): tier = "1.1" counter_tier11 += 1 - tier_dict[key1]["tier 1.1"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 1.1"] += 1 + elif variant_type == "ref": + tier_dict_ref[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])): + all(float(ij) >= 0.75 for ij in [tier1ff, tier2ff, tier3ff, tier4ff])): tier = "1.2" counter_tier12 += 1 - tier_dict[key1]["tier 1.2"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 1.2"] += 1 + elif variant_type == "ref": + tier_dict_ref[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(float(ij) >= 0.75 for ij in [tier1ff, tier4ff])) | (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]))): + all(float(ij) >= 0.75 for ij in [tier2ff, tier3ff]))): tier = "2.1" counter_tier21 += 1 - tier_dict[key1]["tier 2.1"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 2.1"] += 1 + elif variant_type == "ref": + tier_dict_ref[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])): + all(float(ij) >= 0.75 for ij in [tier1ff, tier2ff, tier3ff, tier4ff])): tier = "2.2" counter_tier22 += 1 - tier_dict[key1]["tier 2.2"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 2.2"] += 1 + elif variant_type == "ref": + tier_dict_ref[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(float(ij) >= 0.75 for ij in [tier1ff, tier4ff]) & + any(float(ij) >= 0.75 for ij in [tier2ff, tier3ff])) | (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]))): + all(float(ij) >= 0.75 for ij in [tier2ff, tier3ff]) & + any(float(ij) >= 0.75 for ij in [tier1ff, tier4ff]))): tier = "2.3" counter_tier23 += 1 - tier_dict[key1]["tier 2.3"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 2.3"] += 1 + elif variant_type == "ref": + tier_dict_ref[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(float(ij) >= 0.75 for ij in [tier1ff, tier4ff])) | (all(int(ij) >= 1 for ij in [total2new, total3new]) & - all(float(ij) >= 0.75 for ij in [alt2ff, alt3ff]))): + all(float(ij) >= 0.75 for ij in [tier2ff, tier3ff]))): tier = "2.4" counter_tier24 += 1 - tier_dict[key1]["tier 2.4"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 2.4"] += 1 + elif variant_type == "ref": + tier_dict_ref[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]))): + (all(float(ij) >= 0.5 for ij in [tier1ff, tier4ff]) | + all(float(ij) >= 0.5 for ij in [tier2ff, tier3ff]))): tier = "3.1" counter_tier31 += 1 - tier_dict[key1]["tier 3.1"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 3.1"] += 1 + elif variant_type == "ref": + tier_dict_ref[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(float(ij) >= 0.5 for ij in [tier1ff, tier4ff])) | (all(int(ij) >= 1 for ij in [total2new, total3new]) & - all(float(ij) >= 0.5 for ij in [alt2ff, alt3ff]))): + all(float(ij) >= 0.5 for ij in [tier2ff, tier3ff]))): tier = "3.2" counter_tier32 += 1 - tier_dict[key1]["tier 3.2"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 3.2"] += 1 + elif variant_type == "ref": + tier_dict_ref[key1]["tier 3.2"] += 1 elif (trimmed): tier = "4" counter_tier4 += 1 - tier_dict[key1]["tier 4"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 4"] += 1 + elif variant_type == "ref": + tier_dict_ref[key1]["tier 4"] += 1 # assign tiers if ((all(int(ij) >= 3 for ij in [total1new_trim, total4new_trim]) & - all(float(ij) >= 0.75 for ij in [alt1ff_trim, alt4ff_trim])) | + all(float(ij) >= 0.75 for ij in [tier1ff_trim, tier4ff_trim])) | (all(int(ij) >= 3 for ij in [total2new_trim, total3new_trim]) & - all(float(ij) >= 0.75 for ij in [alt2ff_trim, alt3ff_trim]))): + all(float(ij) >= 0.75 for ij in [tier2ff_trim, tier3ff_trim]))): trimmed_actual_high_tier = True elif (all(int(ij) >= 1 for ij in [total1new_trim, total2new_trim, total3new_trim, total4new_trim]) & any(int(ij) >= 3 for ij in [total1new_trim, total4new_trim]) & any(int(ij) >= 3 for ij in [total2new_trim, total3new_trim]) & - all(float(ij) >= 0.75 for ij in [alt1ff_trim, alt2ff_trim, alt3ff_trim, alt4ff_trim])): + all(float(ij) >= 0.75 for ij in [tier1ff_trim, tier2ff_trim, tier3ff_trim, tier4ff_trim])): trimmed_actual_high_tier = True elif ((all(int(ij) >= 1 for ij in [total1new_trim, total4new_trim]) & any(int(ij) >= 3 for ij in [total1new_trim, total4new_trim]) & - all(float(ij) >= 0.75 for ij in [alt1ff_trim, alt4ff_trim])) | + all(float(ij) >= 0.75 for ij in [tier1ff_trim, tier4ff_trim])) | (all(int(ij) >= 1 for ij in [total2new_trim, total3new_trim]) & any(int(ij) >= 3 for ij in [total2new_trim, total3new_trim]) & - all(float(ij) >= 0.75 for ij in [alt2ff_trim, alt3ff_trim]))): + all(float(ij) >= 0.75 for ij in [tier2ff_trim, tier3ff_trim]))): trimmed_actual_high_tier = True elif (all(int(ij) >= 1 for ij in [total1new_trim, total2new_trim, total3new_trim, total4new_trim]) & - all(float(ij) >= 0.75 for ij in [alt1ff_trim, alt2ff_trim, alt3ff_trim, alt4ff_trim])): + all(float(ij) >= 0.75 for ij in [tier1ff_trim, tier2ff_trim, tier3ff_trim, tier4ff_trim])): trimmed_actual_high_tier = True elif ((all(int(ij) >= 1 for ij in [total1new_trim, total4new_trim]) & any(int(ij) >= 3 for ij in [total2new_trim, total3new_trim]) & - all(float(ij) >= 0.75 for ij in [alt1ff_trim, alt4ff_trim]) & - any(float(ij) >= 0.75 for ij in [alt2ff_trim, alt3ff_trim])) | + all(float(ij) >= 0.75 for ij in [tier1ff_trim, tier4ff_trim]) & + any(float(ij) >= 0.75 for ij in [tier2ff_trim, tier3ff_trim])) | (all(int(ij) >= 1 for ij in [total2new_trim, total3new_trim]) & any(int(ij) >= 3 for ij in [total1new_trim, total4new_trim]) & - all(float(ij) >= 0.75 for ij in [alt2ff_trim, alt3ff_trim]) & - any(float(ij) >= 0.75 for ij in [alt1ff_trim, alt4ff_trim]))): + all(float(ij) >= 0.75 for ij in [tier2ff_trim, tier3ff_trim]) & + any(float(ij) >= 0.75 for ij in [tier1ff_trim, tier4ff_trim]))): trimmed_actual_high_tier = True elif ((all(int(ij) >= 1 for ij in [total1new_trim, total4new_trim]) & - all(float(ij) >= 0.75 for ij in [alt1ff_trim, alt4ff_trim])) | + all(float(ij) >= 0.75 for ij in [tier1ff_trim, tier4ff_trim])) | (all(int(ij) >= 1 for ij in [total2new_trim, total3new_trim]) & - all(float(ij) >= 0.75 for ij in [alt2ff_trim, alt3ff_trim]))): + all(float(ij) >= 0.75 for ij in [tier2ff_trim, tier3ff_trim]))): trimmed_actual_high_tier = True else: trimmed_actual_high_tier = False @@ -1005,37 +1185,58 @@ elif softclipped_mutation_allMates: tier = "5.1" counter_tier51 += 1 - tier_dict[key1]["tier 5.1"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 5.1"] += 1 + elif variant_type == "ref": + tier_dict_ref[key1]["tier 5.1"] += 1 elif softclipped_mutation_oneOfTwoMates: tier = "5.2" counter_tier52 += 1 - tier_dict[key1]["tier 5.2"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 5.2"] += 1 + elif variant_type == "ref": + tier_dict_ref[key1]["tier 5.2"] += 1 elif softclipped_mutation_oneOfTwoSSCS: tier = "5.3" counter_tier53 += 1 - tier_dict[key1]["tier 5.3"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 5.3"] += 1 + elif variant_type == "ref": + tier_dict_ref[key1]["tier 5.3"] += 1 elif softclipped_mutation_oneMate: tier = "5.4" counter_tier54 += 1 - tier_dict[key1]["tier 5.4"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 5.4"] += 1 + elif variant_type == "ref": + tier_dict_ref[key1]["tier 5.4"] += 1 elif softclipped_mutation_oneMateOneSSCS: tier = "5.5" counter_tier55 += 1 - tier_dict[key1]["tier 5.5"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 5.5"] += 1 + elif variant_type == "ref": + tier_dict_ref[key1]["tier 5.5"] += 1 elif (contradictory): tier = "6" counter_tier6 += 1 - tier_dict[key1]["tier 6"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 6"] += 1 + elif variant_type == "ref": + tier_dict_ref[key1]["tier 6"] += 1 else: tier = "7" counter_tier7 += 1 - tier_dict[key1]["tier 7"] += 1 + if variant_type == "alt": + tier_dict[key1]["tier 7"] += 1 + elif variant_type == "ref": + tier_dict_ref[key1]["tier 7"] += 1 chrom, pos, ref_a, alt_a = re.split(r'\#', key1) var_id = '-'.join([chrom, str(int(pos)+1), ref, alt]) @@ -1114,13 +1315,12 @@ 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) + line = (var_id, tier, variant_type, 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) # csv_writer.writerow(line) - line2 = ("", "", 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) + line2 = ("", "", "", 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, line2) # csv_writer.writerow(line2) - # 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), @@ -1155,12 +1355,19 @@ # write to file # move tier 4 counts to tier 2.5 if there other mutations with tier <= 2.4 sum_highTiers = sum([tier_dict[key1][ij] for ij in list(sorted(tier_dict[key1].keys()))[:6]]) + sum_highTiers_ref = sum([tier_dict_ref[key1][ij] for ij in list(sorted(tier_dict_ref[key1].keys()))[:6]]) correct_tier = False + correct_tier_ref = False if tier_dict[key1]["tier 4"] > 0 and sum_highTiers > 0: tier_dict[key1]["tier 2.5"] = tier_dict[key1]["tier 4"] tier_dict[key1]["tier 4"] = 0 correct_tier = True + if tier_dict_ref[key1]["tier 4"] > 0 and sum_highTiers_ref > 0: + tier_dict_ref[key1]["tier 2.5"] = tier_dict_ref[key1]["tier 4"] + tier_dict_ref[key1]["tier 4"] = 0 + correct_tier_ref = True + for sample in change_tier_after_print: row_number = sample[0] line1 = sample[1] @@ -1168,44 +1375,47 @@ actual_high_tier = sample[3] current_tier = list(line1)[1] - if correct_tier and (current_tier == "4") and actual_high_tier: + if line1[2] == "alt" and correct_tier and (current_tier == "4") and actual_high_tier: + line1 = list(line1) + line1[1] = "2.5" + line1 = tuple(line1) + counter_tier25 += 1 + counter_tier4 -= 1 + if line1[2] == "ref" and correct_tier_ref and (current_tier == "4") and actual_high_tier: line1 = list(line1) line1[1] = "2.5" line1 = tuple(line1) counter_tier25 += 1 counter_tier4 -= 1 + ws1.write_row(row_number, 0, line1) csv_writer.writerow(line1) ws1.write_row(row_number + 1, 0, line2) csv_writer.writerow(line2) + if line1[2] == "alt": + ws1.conditional_format('M{}:N{}'.format(row_number + 1, row_number + 2), {'type': 'formula', 'criteria': '=OR($B${}="1.1", $B${}="1.2")'.format(row_number + 1, row_number + 1), 'format': format1, 'multi_range': 'M{}:N{} U{}:V{} B{}'.format(row_number + 1, row_number + 2, row_number + 1, row_number + 2, row_number + 1, row_number + 2)}) + ws1.conditional_format('M{}:N{}'.format(row_number + 1, row_number + 2), {'type': 'formula', 'criteria': '=OR($B${}="2.1", $B${}="2.2", $B${}="2.3", $B${}="2.4", $B${}="2.5")'.format(row_number + 1, row_number + 1, row_number + 1, row_number + 1, row_number + 1), 'format': format3, 'multi_range': 'M{}:N{} U{}:V{} B{}'.format(row_number + 1, row_number + 2, row_number + 1, row_number + 2, row_number + 1, row_number + 2)}) + ws1.conditional_format('M{}:N{}'.format(row_number + 1, row_number + 2), {'type': 'formula', 'criteria': '=$B${}>="3"'.format(row_number + 1), 'format': format2, 'multi_range': 'M{}:N{} U{}:V{} B{}'.format(row_number + 1, row_number + 2, row_number + 1, row_number + 2, row_number + 1, row_number + 2)}) + elif line1[2] == "ref": + ws1.conditional_format('M{}:N{}'.format(row_number + 1, row_number + 2), {'type': 'formula', 'criteria': '=OR($B${}="1.1", $B${}="1.2")'.format(row_number + 1, row_number + 1), 'format': format1, 'multi_range': 'M{}:N{} S{}:T{} B{}'.format(row_number + 1, row_number + 2, row_number + 1, row_number + 2, row_number + 1, row_number + 2)}) + ws1.conditional_format('M{}:N{}'.format(row_number + 1, row_number + 2), {'type': 'formula', 'criteria': '=OR($B${}="2.1", $B${}="2.2", $B${}="2.3", $B${}="2.4", $B${}="2.5")'.format(row_number + 1, row_number + 1, row_number + 1, row_number + 1, row_number + 1), 'format': format3, 'multi_range': 'M{}:N{} S{}:T{} B{}'.format(row_number + 1, row_number + 2, row_number + 1, row_number + 2, row_number + 1, row_number + 2)}) + ws1.conditional_format('M{}:N{}'.format(row_number + 1, row_number + 2), {'type': 'formula', 'criteria': '=$B${}>="3"'.format(row_number + 1), 'format': format2, 'multi_range': 'M{}:N{} S{}:T{} B{}'.format(row_number + 1, row_number + 2, row_number + 1, row_number + 2, row_number + 1, row_number + 2)}) - ws1.conditional_format('L{}:M{}'.format(row_number + 1, row_number + 2), - {'type': 'formula', - 'criteria': '=OR($B${}="1.1", $B${}="1.2")'.format(row_number + 1, row_number + 1), - 'format': format1, - 'multi_range': 'L{}:M{} T{}:U{} B{}'.format(row_number + 1, row_number + 2, row_number + 1, row_number + 2, row_number + 1, row_number + 2)}) - ws1.conditional_format('L{}:M{}'.format(row_number + 1, row_number + 2), - {'type': 'formula', - 'criteria': '=OR($B${}="2.1", $B${}="2.2", $B${}="2.3", $B${}="2.4", $B${}="2.5")'.format(row_number + 1, row_number + 1, row_number + 1, row_number + 1, row_number + 1), - 'format': format3, - 'multi_range': 'L{}:M{} T{}:U{} B{}'.format(row_number + 1, row_number + 2, row_number + 1, row_number + 2, row_number + 1, row_number + 2)}) - ws1.conditional_format('L{}:M{}'.format(row_number + 1, row_number + 2), - {'type': 'formula', - 'criteria': '=$B${}>="3"'.format(row_number + 1), - 'format': format2, - 'multi_range': 'L{}:M{} T{}:U{} B{}'.format(row_number + 1, row_number + 2, row_number + 1, row_number + 2, row_number + 1, row_number + 2)}) # sheet 2 if chimera_correction: - header_line2 = ('variant ID', 'cvrg', 'AC alt (all tiers)', 'AF (all tiers)', 'cvrg (tiers 1.1-2.5)', 'AC alt (tiers 1.1-2.5)', 'AF (tiers 1.1-2.5)', 'chimera-corrected cvrg (tiers 1.1-2.5)', 'chimeras in AC alt (tiers 1.1-2.5)', 'chimera-corrected AF (tiers 1.1-2.5)', '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 2.5', - 'tier 3.1', 'tier 3.2', 'tier 4', 'tier 5.1', 'tier 5.2', 'tier 5.3', 'tier 5.4', 'tier 5.5', 'tier 6', 'tier 7', '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-2.5', 'AF 1.1-3.1', 'AF 1.1-3.2', 'AF 1.1-4', 'AF 1.1-5.1', 'AF 1.1-5.2', 'AF 1.1-5.3', 'AF 1.1-5.4', 'AF 1.1-5.5', 'AF 1.1-6', 'AF 1.1-7') + header_line2 = ('variant ID', 'cvrg', 'AC alt (all tiers)', 'AF (all tiers)', 'cvrg (tiers 1.1-2.5)', 'AC ref (tiers 1.1-2.5)', 'AC alt (tiers 1.1-2.5)', 'AF (tiers 1.1-2.5)', 'chimera-corrected cvrg (tiers 1.1-2.5)', 'chimeras in AC alt (tiers 1.1-2.5)', 'chimera-corrected AF (tiers 1.1-2.5)', 'AC alt (orginal DCS)', 'AF (original DCS)', + 'tier 1.1 (alt)', 'tier 1.2 (alt)', 'tier 2.1 (alt)', 'tier 2.2 (alt)', 'tier 2.3 (alt)', 'tier 2.4 (alt)', 'tier 2.5 (alt)', + 'tier 3.1 (alt)', 'tier 3.2 (alt)', 'tier 4 (alt)', 'tier 5.1 (alt)', 'tier 5.2 (alt)', 'tier 5.3 (alt)', 'tier 5.4 (alt)', 'tier 5.5 (alt)', 'tier 6 (alt)', 'tier 7 (alt)', + 'tier 1.1 (ref)', 'tier 1.2 (ref)', 'tier 2.1 (ref)', 'tier 2.2 (ref)', 'tier 2.3 (ref)', 'tier 2.4 (ref)', 'tier 2.5 (ref)', + 'tier 3.1 (ref)', 'tier 3.2 (ref)', 'tier 4 (ref)', 'tier 5.1 (ref)', 'tier 5.2 (ref)', 'tier 5.3 (ref)', 'tier 5.4 (ref)', 'tier 5.5 (ref)', 'tier 6 (ref)', 'tier 7 (ref)' + ) else: - header_line2 = ('variant ID', 'cvrg', 'AC alt (all tiers)', 'AF (all tiers)', 'cvrg (tiers 1.1-2.5)', 'AC alt (tiers 1.1-2.5)', 'AF (tiers 1.1-2.5)', '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 2.5', - 'tier 3.1', 'tier 3.2', 'tier 4', 'tier 5.1', 'tier 5.2', 'tier 5.3', 'tier 5.4', 'tier 5.5', 'tier 6', 'tier 7', '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-2.5', 'AF 1.1-3.1', 'AF 1.1-3.2', 'AF 1.1-4', 'AF 1.1-5.1', 'AF 1.1-5.2', 'AF 1.1-5.3', 'AF 1.1-5.4', 'AF 1.1-5.5', 'AF 1.1-6', 'AF 1.1-7') - + header_line2 = ('variant ID', 'cvrg', 'AC alt (all tiers)', 'AF (all tiers)', 'cvrg (tiers 1.1-2.5)', 'AC ref (tiers 1.1-2.5)', 'AC alt (tiers 1.1-2.5)', 'AF (tiers 1.1-2.5)', 'AC alt (orginal DCS)', 'AF (original DCS)', + 'tier 1.1 (alt)', 'tier 1.2 (alt)', 'tier 2.1 (alt)', 'tier 2.2 (alt)', 'tier 2.3 (alt)', 'tier 2.4 (alt)', 'tier 2.5 (alt)', + 'tier 3.1 (alt)', 'tier 3.2 (alt)', 'tier 4 (alt)', 'tier 5.1 (alt)', 'tier 5.2 (alt)', 'tier 5.3 (alt)', 'tier 5.4 (alt)', 'tier 5.5 (alt)', 'tier 6 (alt)', 'tier 7 (alt)', + 'tier 1.1 (ref)', 'tier 1.2 (ref)', 'tier 2.1 (ref)', 'tier 2.2 (ref)', 'tier 2.3 (ref)', 'tier 2.4 (ref)', 'tier 2.5 (ref)', + 'tier 3.1 (ref)', 'tier 3.2 (ref)', 'tier 4 (ref)', 'tier 5.1 (ref)', 'tier 5.2 (ref)', 'tier 5.3 (ref)', 'tier 5.4 (ref)', 'tier 5.5 (ref)', 'tier 6 (ref)', 'tier 7 (ref)' + ) ws2.write_row(0, 0, header_line2) row = 0 @@ -1220,9 +1430,12 @@ alt_count = cvrg_dict[key1][1] cvrg = ref_count + alt_count + ref_tiers = tier_dict_ref[key1] + var_id = '-'.join([chrom, str(int(pos)+1), ref, alt]) lst = [var_id, cvrg] used_tiers = [] + used_tiers_ref = [t for k, t in sorted(ref_tiers.items())] cum_af = [] for key2, value2 in sorted(value1.items()): # calculate cummulative AF @@ -1233,24 +1446,25 @@ if sum(used_tiers) == 0: # skip mutations that are filtered by the VA in the first place continue lst.extend([sum(used_tiers), safe_div(sum(used_tiers), cvrg)]) - lst.extend([(cvrg - sum(used_tiers[-10:])), sum(used_tiers[0:7]), safe_div(sum(used_tiers[0:7]), (cvrg - sum(used_tiers[-10:])))]) + lst.extend([(sum(used_tiers_ref[0:7]) + sum(used_tiers[0:7])), sum(used_tiers_ref[0:7]), sum(used_tiers[0:7]), safe_div(sum(used_tiers[0:7]), (sum(used_tiers_ref[0:7]) + sum(used_tiers[0:7])))]) if chimera_correction: chimeras_all = chimera_dict[key1][1] new_alt = sum(used_tiers[0:7]) - chimeras_all fraction_chimeras = safe_div(chimeras_all, float(sum(used_tiers[0:7]))) if fraction_chimeras is None: fraction_chimeras = 0. - new_cvrg = (cvrg - sum(used_tiers[-10:])) * (1. - fraction_chimeras) + new_cvrg = (sum(used_tiers_ref[0:7]) + sum(used_tiers[0:7])) * (1. - fraction_chimeras) lst.extend([new_cvrg, chimeras_all, safe_div(new_alt, new_cvrg)]) lst.extend([alt_count, safe_div(alt_count, cvrg)]) lst.extend(used_tiers) - lst.extend(cum_af) + lst.extend(used_tiers_ref) + # lst.extend(cum_af) lst = tuple(lst) ws2.write_row(row + 1, 0, lst) if chimera_correction: - ws2.conditional_format('M{}:N{}'.format(row + 2, row + 2), {'type': 'formula', 'criteria': '=$M$1="tier 1.1"', 'format': format12, 'multi_range': 'M{}:N{} M1:N1'.format(row + 2, row + 2)}) - ws2.conditional_format('O{}:S{}'.format(row + 2, row + 2), {'type': 'formula', 'criteria': '=$O$1="tier 2.1"', 'format': format32, 'multi_range': 'O{}:S{} O1:S1'.format(row + 2, row + 2)}) - ws2.conditional_format('T{}:AC{}'.format(row + 2, row + 2), {'type': 'formula', 'criteria': '=$T$1="tier 3.1"', 'format': format22, 'multi_range': 'T{}:AC{} T1:AC1'.format(row + 2, row + 2)}) + ws2.conditional_format('N{}:O{}'.format(row + 2, row + 2), {'type': 'formula', 'criteria': '=$N$1="tier 1.1 (alt)"', 'format': format12, 'multi_range': 'N{}:O{} N1:O1 AE{}:AF{} AE1:AF1'.format(row + 2, row + 2, row + 2, row + 2)}) + ws2.conditional_format('P{}:T{}'.format(row + 2, row + 2), {'type': 'formula', 'criteria': '=$P$1="tier 2.1 (alt)"', 'format': format32, 'multi_range': 'P{}:T{} P1:T1 AG{}:AK{} AG1:AK1'.format(row + 2, row + 2, row + 2, row + 2)}) + ws2.conditional_format('U{}:AD{}'.format(row + 2, row + 2), {'type': 'formula', 'criteria': '=$U$1="tier 3.1 (alt)"', 'format': format22, 'multi_range': 'U{}:AD{} U1:AD1 AL{}:AU{} AL1:AU1'.format(row + 2, row + 2, row + 2, row + 2)}) else: ws2.conditional_format('J{}:K{}'.format(row + 2, row + 2), {'type': 'formula', 'criteria': '=$J$1="tier 1.1"', 'format': format12, 'multi_range': 'J{}:K{} J1:K1'.format(row + 2, row + 2)}) ws2.conditional_format('L{}:P{}'.format(row + 2, row + 2), {'type': 'formula', 'criteria': '=$L$1="tier 2.1"', 'format': format32, 'multi_range': 'L{}:P{} L1:P1'.format(row + 2, row + 2)})
--- a/read2mut.xml Mon Mar 29 09:22:57 2021 +0000 +++ b/read2mut.xml Fri Jul 22 09:19:44 2022 +0000 @@ -1,6 +1,6 @@ <?xml version="1.0" encoding="UTF-8"?> -<tool id="read2mut" name="Call specific mutations in reads:" version="2.1.4" profile="19.01"> - <description>Looks for reads with mutation at known positions and calculates frequencies and stats.</description> +<tool id="read2mut" name="Call specific mutations in reads:" version="3.0.0" profile="19.01"> + <description>Looks for reads with a mutation at known positions and calculates frequencies and stats.</description> <macros> <import>va_macros.xml</import> </macros> @@ -28,12 +28,12 @@ ]]> </command> <inputs> - <param name="file1" type="data" format="vcf" label="DCS Mutation File" optional="false" help="VCF file with DCS mutations. See Help section below for a detailed explanation."/> + <param name="file1" type="data" format="vcf" label="DCS Mutation File" optional="false" help="VCF file with DCS mutations. See the Help section below for a detailed explanation."/> <param name="file2" type="data" format="bam" label="BAM File of raw reads" optional="false" help="BAM file with aligned raw reads of selected tags."/> <param name="file3" type="data" format="json" label="JSON File with DCS tag stats" optional="false" help="JSON file generated by DCS mutations to tags/reads"/> <param name="file4" type="data" format="json" label="JSON File with SSCS tag stats" optional="false" help="JSON file generated by DCS mutations to SSCS stats."/> - <param name="thresh" type="integer" label="Tag count threshold" value="0" help="Integer threshold for displaying mutations. Only mutations occuring in DCS of less than thresh tags are displayed. Default of 0 displays all."/> - <param name="phred" type="integer" label="Phred quality score threshold" min="0" max="41" value="20" help="Integer threshold for Phred quality score. Only reads higher than this threshold are considered. Default = 20."/> + <param name="thresh" type="integer" label="Tag count threshold" value="0" help="Integer threshold for displaying mutations. Only mutations occurring in DCS of less than thresh tags are displayed. Default of 0 displays all."/> + <param name="phred" type="integer" label="Phred quality score threshold" min="0" max="41" value="20" help="Integer threshold for Phred quality score. Only reads higher than this threshold is considered. Default = 20."/> <param name="trim" type="integer" label="Trimming threshold" value="10" help="Integer threshold for assigning mutations at start and end of reads to lower tier. Default 10."/> <param name="chimera_correction" type="boolean" label="Apply chimera correction?" truevalue="--chimera_correction" falsevalue="" checked="False" help="Count chimeric variants and correct the variant frequencies."/> <param name="softclipping_dist" type="integer" label="Distance between artifact and softclipping of the reads" min="1" value="15" help="Count mutation as an artifact if mutation lies within this parameter away from the softclipping part of the reads. Default = 20"/> @@ -57,10 +57,10 @@ <param name="chimera_correction"/> <param name="softclipping_dist" value="15"/> <param name="reads_threshold" value="1.0"/> - <output name="output_xlsx" file="Variant_Analyzer_summary_test.xlsx" decompress="true" lines_diff="10"/> - <output name="outputFile_csv" file="Variant_Analyzer_summary_test.csv" decompress="true" lines_diff="10"/> - <output name="output_xlsx2" file="Variant_Analyzer_allele_frequencies_test.xlsx" decompress="true" lines_diff="10"/> - <output name="output_xlsx3" file="Variant_Analyzer_tiers_test.xlsx" decompress="true" lines_diff="10"/> + <output name="output_xlsx" file="Variant_Analyzer_summary_test.xlsx" decompress="true"/> + <output name="outputFile_csv" file="Variant_Analyzer_summary_test.csv" decompress="true"/> + <output name="output_xlsx2" file="Variant_Analyzer_allele_frequencies_test.xlsx" decompress="true"/> + <output name="output_xlsx3" file="Variant_Analyzer_tiers_test.xlsx" decompress="true"/> </test> </tests> <help> <![CDATA[ @@ -73,25 +73,25 @@ **Input** -**Dataset 1:** VCF file with duplex consesus sequence (DCS) mutations. E.g. +**Dataset 1:** VCF file with duplex consensus sequence (DCS) mutations. E.g. generated by the `FreeBayes <https://arxiv.org/abs/1207.3907>`_ or `LoFreq <https://academic.oup.com/nar/article/40/22/11189/1152727>`_ variant caller. **Dataset 2:** BAM file of aligned raw reads. This file can be obtained by the tool `Map with BWA-MEM <https://arxiv.org/abs/1303.3997>`_. **Dataset 3:** JSON file generated by the **DCS mutations to tags/reads** tool -containing dictonaries of the tags of reads containing mutations +containing dictionaries 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 +stats of tags that carry a mutation and the reference allele in the SSCS at the same position a mutation is called in the DCS. **Output** -The output are three XLSX files containing frequencies stats for DCS mutations based -on information from the raw reads and a CSV file containing the summary information without color-coding. In addition to that a tier based -classification is provided based on the amout of support for a true variant call. +The output is three XLSX files containing frequencies stats for DCS mutations based +on information from the raw reads and a CSV file containing the summary information without color-coding. In addition to that, a tier-based +classification is provided based on the amount of support for a true variant call. ]]>
--- a/test-data/Variant_Analyzer_summary_test.csv Mon Mar 29 09:22:57 2021 +0000 +++ b/test-data/Variant_Analyzer_summary_test.csv Fri Jul 22 09:19:44 2022 +0000 @@ -1,7 +1,11 @@ -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,in phase,chimeric tag -ACH_TDII_5regions-505-C-A,2.1,GATAACCTTGCTTCGTGATTAATC,ab1.ba2,132.0,131.0,264.0,263.0,173.0,1,3,1,3,0,0,1,3,0,0,1.0,1.0,0,0,0,0,0,0,2,1,1,1,, -,,GATAACCTTGCTTCGTGATTAATC,ab2.ba1,,,,,173.0,0,0,0,0,0,0,0,0,,,,,0,0,0,0,0,0,2,1,1,1,, -ACH_TDII_5regions-571-C-T,1.1,GATTGGATAACGTTGTGGCAATTG,ab1.ba2,129.0,218.0,195.0,284.0,143.0,4,5,4,5,0,0,4,5,0,0,1.0,1.0,0,0,0,0,0,0,1,1,2,1,, -,,GATTGGATAACGTTGTGGCAATTG,ab2.ba1,264.0,264.0,278.5,283.5,143.0,2,2,2,2,0,0,2,2,0,0,1.0,1.0,0,0,0,0,0,0,1,1,2,1,, -ACH_TDII_5regions-958-T-C,2.4,CCTCCCGGCAGTGCGAAAATGTCA,ab1.ba2,,,,,195.0,0,0,0,0,0,0,0,0,,,,,0,0,0,0,0,0,1,1,1,0,, -,,CCTCCCGGCAGTGCGAAAATGTCA,ab2.ba1,97.0,91.0,283.0,277.0,195.0,1,1,1,1,0,0,1,1,0,0,1.0,1.0,0,0,0,0,0,0,1,1,1,0,, +variant ID,tier,allele,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,in phase,chimeric tag +ACH_TDII_5regions-505-C-A,2.1,alt,GATAACCTTGCTTCGTGATTAATC,ab1.ba2,132.0,131.0,264.0,263.0,173.0,1,3,1,3,0,0,1,3,0,0,1.0,1.0,0,0,0,0,0,0,2,1,1,1,, +,,,GATAACCTTGCTTCGTGATTAATC,ab2.ba1,,,,,173.0,0,0,0,0,0,0,0,0,,,,,0,0,0,0,0,0,2,1,1,1,, +ACH_TDII_5regions-505-C-A,1.1,ref,GATTGGATAACGTTGTGGCAATTG,ab1.ba2,73.0,152.0,205.0,284.0,173.0,3,5,3,5,3,5,0,0,1.0,1.0,0,0,0,0,0,0,0,0,2,1,1,1,, +,,,GATTGGATAACGTTGTGGCAATTG,ab2.ba1,198.0,198.0,263.5,283.0,173.0,4,3,4,3,4,3,0,0,1.0,1.0,0,0,0,0,0,0,0,0,2,1,1,1,, +ACH_TDII_5regions-571-C-T,2.1,ref,GATAACCTTGCTTCGTGATTAATC,ab1.ba2,198.0,197.0,264.0,263.0,143.0,1,3,1,3,1,3,0,0,1.0,1.0,0,0,0,0,0,0,0,0,1,1,2,1,, +,,,GATAACCTTGCTTCGTGATTAATC,ab2.ba1,,,,,143.0,0,0,0,0,0,0,0,0,,,,,0,0,0,0,0,0,1,1,2,1,, +ACH_TDII_5regions-571-C-T,1.1,alt,GATTGGATAACGTTGTGGCAATTG,ab1.ba2,129.0,218.0,195.0,284.0,143.0,4,5,4,5,0,0,4,5,0,0,1.0,1.0,0,0,0,0,0,0,1,1,2,1,, +,,,GATTGGATAACGTTGTGGCAATTG,ab2.ba1,264.0,264.0,278.5,283.5,143.0,2,2,2,2,0,0,2,2,0,0,1.0,1.0,0,0,0,0,0,0,1,1,2,1,, +ACH_TDII_5regions-958-T-C,2.4,alt,CCTCCCGGCAGTGCGAAAATGTCA,ab1.ba2,,,,,195.0,0,0,0,0,0,0,0,0,,,,,0,0,0,0,0,0,1,1,1,0,, +,,,CCTCCCGGCAGTGCGAAAATGTCA,ab2.ba1,97.0,91.0,283.0,277.0,195.0,1,1,1,1,0,0,1,1,0,0,1.0,1.0,0,0,0,0,0,0,1,1,1,0,,
--- a/test-data/tag_count_dict_test.json Mon Mar 29 09:22:57 2021 +0000 +++ b/test-data/tag_count_dict_test.json Fri Jul 22 09:19:44 2022 +0000 @@ -1,1 +1,1 @@ -[{"GATAACCTTGCTTCGTGATTAATC": {"ACH_TDII_5regions#504#C#A": "A"}, "GATTGGATAACGTTGTGGCAATTG": {"ACH_TDII_5regions#570#C#T": "T"}, "CCTCCCGGCAGTGCGAAAATGTCA": {"ACH_TDII_5regions#957#T#C": "C"}}, {"ACH_TDII_5regions#957#T#C": [0, 1, 195.0], "ACH_TDII_5regions#504#C#A": [1, 1, 173.0], "ACH_TDII_5regions#570#C#T": [1, 1, 143.0]}] \ No newline at end of file +[{"GATAACCTTGCTTCGTGATTAATC": {"ACH_TDII_5regions#504#C#A": "A"}, "GATTGGATAACGTTGTGGCAATTG": {"ACH_TDII_5regions#570#C#T": "T"}, "CCTCCCGGCAGTGCGAAAATGTCA": {"ACH_TDII_5regions#957#T#C": "C"}}, {"ACH_TDII_5regions#957#T#C": [0, 1, 195.0], "ACH_TDII_5regions#504#C#A": [1, 1, 173.0], "ACH_TDII_5regions#570#C#T": [1, 1, 143.0]}, {"GATAACCTTGCTTCGTGATTAATC": {"ACH_TDII_5regions#570#C#T": "C"}, "GATTGGATAACGTTGTGGCAATTG": {"ACH_TDII_5regions#504#C#A": "C"}}] \ No newline at end of file