diff --git a/CHANGELOG.md b/CHANGELOG.md index a54814f2cb7affcf8472c6f50ae213bb5eef1202..e975fc257b08f53ab1e5081b01d3f397126d999d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## Unreleased +### Added +- Option in merge_vcf_files.py to merge VCF files of different samples without genotyping based on read depth computed by duphold +- Option in vcf_to_table.py to calculate number of deletions that are supported by a change in read depth compared to 1000 bp flanking regions ## [0.3.0] - 2019-12-16 ### Added diff --git a/scripts/.idea/workspace.xml b/scripts/.idea/workspace.xml index 129471360ed376b28b3bbfdbdc2f4cec2c6defdd..39ed2b8058a20e2de6d60ae898e781e7b7085cac 100644 --- a/scripts/.idea/workspace.xml +++ b/scripts/.idea/workspace.xml @@ -3,6 +3,8 @@ <component name="ChangeListManager"> <list default="true" id="38800356-6f91-47c6-8c6a-4d9b650fe843" name="Default" comment=""> <change beforePath="$PROJECT_DIR$/.idea/workspace.xml" afterPath="$PROJECT_DIR$/.idea/workspace.xml" /> + <change beforePath="$PROJECT_DIR$/convert/vcf_to_table.py" afterPath="$PROJECT_DIR$/convert/vcf_to_table.py" /> + <change beforePath="$PROJECT_DIR$/genotype/merge_vcf_files.py" afterPath="$PROJECT_DIR$/genotype/merge_vcf_files.py" /> </list> <option name="EXCLUDED_CONVERTED_TO_IGNORED" value="true" /> <option name="TRACKING_ENABLED" value="true" /> @@ -11,62 +13,26 @@ <option name="HIGHLIGHT_NON_ACTIVE_CHANGELIST" value="false" /> <option name="LAST_RESOLUTION" value="IGNORE" /> </component> + <component name="DatabaseView"> + <option name="SHOW_INTERMEDIATE" value="true" /> + <option name="GROUP_DATA_SOURCES" value="true" /> + <option name="GROUP_SCHEMA" value="true" /> + <option name="GROUP_CONTENTS" value="false" /> + <option name="SORT_POSITIONED" value="false" /> + <option name="SHOW_TABLE_DETAILS" value="true" /> + <option name="SHOW_EMPTY_GROUPS" value="false" /> + <option name="AUTO_SCROLL_FROM_SOURCE" value="false" /> + <expand /> + <select /> + </component> <component name="FileEditorManager"> <leaf SIDE_TABS_SIZE_LIMIT_KEY="300"> - <file leaf-file-name="merge_vcf_files.py" pinned="false" current-in-tab="false"> - <entry file="file://$PROJECT_DIR$/genotype/merge_vcf_files.py"> - <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="7272"> - <caret line="414" column="31" lean-forward="false" selection-start-line="414" selection-start-column="4" selection-end-line="414" selection-end-column="31" /> - <folding> - <element signature="e#132#147#0" expanded="true" /> - </folding> - </state> - </provider> - </entry> - </file> - <file leaf-file-name="predict_cnvs_model_insertions_probabilities.py" pinned="false" current-in-tab="false"> - <entry file="file://$PROJECT_DIR$/predict/predict_cnvs_model_insertions_probabilities.py"> - <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="1026"> - <caret line="60" column="8" lean-forward="false" selection-start-line="60" selection-start-column="8" selection-end-line="60" selection-end-column="8" /> - <folding> - <element signature="e#112#131#0" expanded="true" /> - </folding> - </state> - </provider> - </entry> - </file> - <file leaf-file-name="process_simple_cnvs.py" pinned="false" current-in-tab="false"> - <entry file="file://$PROJECT_DIR$/process/process_simple_cnvs.py"> - <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="936"> - <caret line="59" column="4" lean-forward="false" selection-start-line="59" selection-start-column="4" selection-end-line="59" selection-end-column="4" /> - <folding> - <element signature="e#120#135#0" expanded="true" /> - </folding> - </state> - </provider> - </entry> - </file> - <file leaf-file-name="internally_collapse_bedpe_intervals.py" pinned="false" current-in-tab="false"> - <entry file="file://$PROJECT_DIR$/collapse/internally_collapse_bedpe_intervals.py"> - <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="3636"> - <caret line="210" column="40" lean-forward="false" selection-start-line="210" selection-start-column="40" selection-end-line="210" selection-end-column="40" /> - <folding /> - </state> - </provider> - </entry> - </file> - <file leaf-file-name="filter_ref_sites_vcf.py" pinned="false" current-in-tab="true"> + <file leaf-file-name="filter_ref_sites_vcf.py" pinned="false" current-in-tab="false"> <entry file="file://$PROJECT_DIR$/filter/filter_ref_sites_vcf.py"> <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="18"> + <state relative-caret-position="198"> <caret line="11" column="0" lean-forward="false" selection-start-line="11" selection-start-column="0" selection-end-line="11" selection-end-column="0" /> - <folding> - <element signature="e#154#169#0" expanded="true" /> - </folding> + <folding /> </state> </provider> </entry> @@ -93,11 +59,23 @@ </provider> </entry> </file> + <file leaf-file-name="merge_vcf_files.py" pinned="false" current-in-tab="true"> + <entry file="file://$PROJECT_DIR$/genotype/merge_vcf_files.py"> + <provider selected="true" editor-type-id="text-editor"> + <state relative-caret-position="521"> + <caret line="102" column="34" lean-forward="true" selection-start-line="102" selection-start-column="34" selection-end-line="102" selection-end-column="34" /> + <folding> + <element signature="e#132#147#0" expanded="true" /> + </folding> + </state> + </provider> + </entry> + </file> <file leaf-file-name="vcf_to_table.py" pinned="false" current-in-tab="false"> <entry file="file://$PROJECT_DIR$/convert/vcf_to_table.py"> <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="2376"> - <caret line="135" column="0" lean-forward="false" selection-start-line="135" selection-start-column="0" selection-end-line="145" selection-end-column="0" /> + <state relative-caret-position="1602"> + <caret line="152" column="78" lean-forward="false" selection-start-line="152" selection-start-column="78" selection-end-line="152" selection-end-column="78" /> <folding> <element signature="e#108#123#0" expanded="true" /> </folding> @@ -144,6 +122,7 @@ <find>structural_variant_from</find> <find>keys</find> <find>.keys</find> + <find>genotyping</find> </findStrings> </component> <component name="Git.Settings"> @@ -178,12 +157,12 @@ <option value="$PROJECT_DIR$/convert/intersect_bedpe_to_feature_bedpe.py" /> <option value="$PROJECT_DIR$/filter/concat_duphold_to_calls.py" /> <option value="$PROJECT_DIR$/convert/bedpe_to_vcf.py" /> - <option value="$PROJECT_DIR$/convert/vcf_to_table.py" /> - <option value="$PROJECT_DIR$/genotype/merge_vcf_files.py" /> <option value="$PROJECT_DIR$/collapse/internally_collapse_bedpe_intervals.py" /> <option value="$PROJECT_DIR$/intersect/intersecting_bedpe_intervals.py" /> <option value="$PROJECT_DIR$/filter/filter_ref_sites_vcf.py" /> <option value="$PROJECT_DIR$/convert/vcf_to_site_only_vcf.py" /> + <option value="$PROJECT_DIR$/convert/vcf_to_table.py" /> + <option value="$PROJECT_DIR$/genotype/merge_vcf_files.py" /> </list> </option> </component> @@ -215,7 +194,6 @@ <foldersAlwaysOnTop value="true" /> </navigator> <panes> - <pane id="Scratches" /> <pane id="Scope" /> <pane id="ProjectPane"> <subPane> @@ -278,12 +256,13 @@ <select /> </subPane> </pane> + <pane id="Scratches" /> </panes> </component> <component name="PropertiesComponent"> <property name="nodejs_interpreter_path.stuck_in_default_project" value="undefined stuck path" /> <property name="WebServerToolWindowFactoryState" value="false" /> - <property name="last_opened_file_path" value="$USER_HOME$/Code/PycharmProjects/Polyhetero" /> + <property name="last_opened_file_path" value="$PROJECT_DIR$" /> </component> <component name="RunDashboard"> <option name="ruleStates"> @@ -323,10 +302,10 @@ <window_info id="Python Console" active="false" anchor="bottom" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="false" show_stripe_button="true" weight="0.32914045" sideWeight="0.5" order="7" side_tool="false" content_ui="tabs" /> <window_info id="Run" active="false" anchor="bottom" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="false" show_stripe_button="true" weight="0.33" sideWeight="0.5" order="2" side_tool="false" content_ui="tabs" /> <window_info id="Terminal" active="false" anchor="bottom" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="false" show_stripe_button="true" weight="0.32914045" sideWeight="0.5" order="7" side_tool="false" content_ui="tabs" /> - <window_info id="Project" active="false" anchor="left" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="true" show_stripe_button="true" weight="0.42999446" sideWeight="0.5" order="0" side_tool="false" content_ui="combo" /> + <window_info id="Project" active="false" anchor="left" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="true" show_stripe_button="true" weight="0.44659656" sideWeight="0.5" order="0" side_tool="false" content_ui="combo" /> <window_info id="Docker" active="false" anchor="bottom" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="false" show_stripe_button="false" weight="0.33" sideWeight="0.5" order="7" side_tool="false" content_ui="tabs" /> - <window_info id="Database" active="false" anchor="right" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="false" show_stripe_button="true" weight="0.33" sideWeight="0.5" order="3" side_tool="false" content_ui="tabs" /> - <window_info id="SciView" active="false" anchor="right" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="false" show_stripe_button="true" weight="0.33" sideWeight="0.5" order="3" side_tool="false" content_ui="tabs" /> + <window_info id="Database" active="false" anchor="right" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="false" show_stripe_button="true" weight="0.32982844" sideWeight="0.5125786" order="4" side_tool="true" content_ui="tabs" /> + <window_info id="SciView" active="false" anchor="right" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="false" show_stripe_button="true" weight="0.32982844" sideWeight="0.4874214" order="3" side_tool="false" content_ui="tabs" /> <window_info id="Structure" active="false" anchor="left" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="false" show_stripe_button="true" weight="0.25" sideWeight="0.5" order="1" side_tool="false" content_ui="tabs" /> <window_info id="Favorites" active="false" anchor="left" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="false" show_stripe_button="true" weight="0.33" sideWeight="0.5" order="2" side_tool="true" content_ui="tabs" /> <window_info id="Debug" active="false" anchor="bottom" auto_hide="false" internal_type="DOCKED" type="DOCKED" visible="false" show_stripe_button="true" weight="0.4" sideWeight="0.5" order="3" side_tool="false" content_ui="tabs" /> @@ -598,40 +577,31 @@ </state> </provider> </entry> - <entry file="file://$PROJECT_DIR$/genotype/merge_vcf_files.py"> + <entry file="file://$PROJECT_DIR$/collapse/internally_collapse_bedpe_intervals.py"> <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="7272"> - <caret line="414" column="31" lean-forward="false" selection-start-line="414" selection-start-column="4" selection-end-line="414" selection-end-column="31" /> - <folding> - <element signature="e#132#147#0" expanded="true" /> - </folding> + <state relative-caret-position="866"> + <caret line="210" column="40" lean-forward="false" selection-start-line="210" selection-start-column="40" selection-end-line="210" selection-end-column="40" /> </state> </provider> </entry> - <entry file="file://$PROJECT_DIR$/predict/predict_cnvs_model_insertions_probabilities.py"> + <entry file="file://$PROJECT_DIR$/process/process_simple_cnvs.py"> <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="1026"> - <caret line="60" column="8" lean-forward="false" selection-start-line="60" selection-start-column="8" selection-end-line="60" selection-end-column="8" /> - <folding> - <element signature="e#112#131#0" expanded="true" /> - </folding> + <state relative-caret-position="866"> + <caret line="59" column="4" lean-forward="false" selection-start-line="59" selection-start-column="4" selection-end-line="59" selection-end-column="4" /> </state> </provider> </entry> - <entry file="file://$PROJECT_DIR$/process/process_simple_cnvs.py"> + <entry file="file://$PROJECT_DIR$/predict/predict_cnvs_model_insertions_probabilities.py"> <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="936"> - <caret line="59" column="4" lean-forward="false" selection-start-line="59" selection-start-column="4" selection-end-line="59" selection-end-column="4" /> - <folding> - <element signature="e#120#135#0" expanded="true" /> - </folding> + <state relative-caret-position="866"> + <caret line="60" column="8" lean-forward="false" selection-start-line="60" selection-start-column="8" selection-end-line="60" selection-end-column="8" /> </state> </provider> </entry> - <entry file="file://$PROJECT_DIR$/collapse/internally_collapse_bedpe_intervals.py"> + <entry file="file://$PROJECT_DIR$/filter/filter_ref_sites_vcf.py"> <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="3636"> - <caret line="210" column="40" lean-forward="false" selection-start-line="210" selection-start-column="40" selection-end-line="210" selection-end-column="40" /> + <state relative-caret-position="198"> + <caret line="11" column="0" lean-forward="false" selection-start-line="11" selection-start-column="0" selection-end-line="11" selection-end-column="0" /> <folding /> </state> </provider> @@ -646,6 +616,16 @@ </state> </provider> </entry> + <entry file="file://$PROJECT_DIR$/gridss/fix_gridss_header.py"> + <provider selected="true" editor-type-id="text-editor"> + <state relative-caret-position="648"> + <caret line="38" column="34" lean-forward="false" selection-start-line="38" selection-start-column="34" selection-end-line="38" selection-end-column="34" /> + <folding> + <element signature="e#68#83#0" expanded="true" /> + </folding> + </state> + </provider> + </entry> <entry file="file://$PROJECT_DIR$/convert/vcf_to_bedpe.py"> <provider selected="true" editor-type-id="text-editor"> <state relative-caret-position="108"> @@ -656,30 +636,20 @@ </entry> <entry file="file://$PROJECT_DIR$/convert/vcf_to_table.py"> <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="2376"> - <caret line="135" column="0" lean-forward="false" selection-start-line="135" selection-start-column="0" selection-end-line="145" selection-end-column="0" /> + <state relative-caret-position="1602"> + <caret line="152" column="78" lean-forward="false" selection-start-line="152" selection-start-column="78" selection-end-line="152" selection-end-column="78" /> <folding> <element signature="e#108#123#0" expanded="true" /> </folding> </state> </provider> </entry> - <entry file="file://$PROJECT_DIR$/gridss/fix_gridss_header.py"> - <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="648"> - <caret line="38" column="34" lean-forward="false" selection-start-line="38" selection-start-column="34" selection-end-line="38" selection-end-column="34" /> - <folding> - <element signature="e#68#83#0" expanded="true" /> - </folding> - </state> - </provider> - </entry> - <entry file="file://$PROJECT_DIR$/filter/filter_ref_sites_vcf.py"> + <entry file="file://$PROJECT_DIR$/genotype/merge_vcf_files.py"> <provider selected="true" editor-type-id="text-editor"> - <state relative-caret-position="18"> - <caret line="11" column="0" lean-forward="false" selection-start-line="11" selection-start-column="0" selection-end-line="11" selection-end-column="0" /> + <state relative-caret-position="521"> + <caret line="102" column="34" lean-forward="true" selection-start-line="102" selection-start-column="34" selection-end-line="102" selection-end-column="34" /> <folding> - <element signature="e#154#169#0" expanded="true" /> + <element signature="e#132#147#0" expanded="true" /> </folding> </state> </provider> diff --git a/scripts/convert/vcf_to_table.py b/scripts/convert/vcf_to_table.py index 471237003b756491935e6477288a1a1e85b960e5..1bd2bda6a94eb3e043636672e69c10e34c4c9475 100755 --- a/scripts/convert/vcf_to_table.py +++ b/scripts/convert/vcf_to_table.py @@ -47,13 +47,15 @@ def vcf_to_table(input_fn, output_fn, fields=None, genotype_fields=None): samples = list(vcf.header.samples) samples.sort() header_fields = [] - standard_fields = ["CHROM", "POS", "REF", "ALT", "QUAL", "FILTER", "END", "HOM-VAR", "VAR"] + standard_fields = ["CHROM", "POS", "REF", "ALT", "QUAL", "FILTER", "END", "HOM-VAR", "VAR", "DEL_SUPPORTED"] standard_fields.extend(vcf.header.info) if fields: for field in fields: if field not in standard_fields: raise ValueError( "{} not defined as standard or INFO field".format(field)) + if field == "DEL_SUPPORTED": + header_fields.extend(["DEL_SUPPORTED", "DEL_UNSUPPORTED", "NON_DEL_SUPPORTED", "NON_DEL_UNSUPPORTED"]) else: header_fields.append(field) if genotype_fields: @@ -116,6 +118,39 @@ def vcf_to_table(input_fn, output_fn, fields=None, genotype_fields=None): if genotype not in non_variants: var_calls += 1 record_fields.append(str(var_calls)) + elif field == "DEL_SUPPORTED": + del_supported_calls = 0 + del_unsupported_calls = 0 + nondel_supported_calls = 0 + nondel_unsupported_calls = 0 + + if record.info["SVTYPE"] != "DEL": + record_fields.extend(["NA", "NA", "NA", "NA"]) + else: + for sample in samples: + genotype = record.samples[sample]["GT"] + dhffc = record.samples[sample]["DHFFC"] + non_variants = [(0, 0), (None, None), + (None, 0)] + if genotype in non_variants: + if dhffc >= 0.75: + nondel_supported_calls += 1 + else: + nondel_unsupported_calls += 1 + elif genotype == (0, 1): + if dhffc >= 0.25 and dhffc < 0.75: + del_supported_calls += 1 + else: + del_unsupported_calls += 1 + elif genotype == (1, 1): + if dhffc >= 0 and dhffc < 0.25: + del_supported_calls += 1 + else: + del_unsupported_calls += 1 + record_fields.append(str(del_supported_calls)) + record_fields.append(str(del_unsupported_calls)) + record_fields.append(str(nondel_supported_calls)) + record_fields.append(str(nondel_unsupported_calls)) elif field in record.info: record_fields.append(str(record.info[field])) else: diff --git a/scripts/genotype/merge_vcf_files.py b/scripts/genotype/merge_vcf_files.py index 512384745aa6c4fa607a4bb183c5b0b08837e6f8..fccaf9cd87b7ff00c36b73ee479d50a2515faf4f 100755 --- a/scripts/genotype/merge_vcf_files.py +++ b/scripts/genotype/merge_vcf_files.py @@ -38,21 +38,23 @@ def parse_cl_args(in_args): parser.add_argument("-r", "--reciprocal", type=float, default=0.5, help="Minimum fraction of reciprocal overlap needed " "to collapse calls") - parser.add_argument("-he", "--heterozygous", type=bool, default=True, - help="Boolean indicating whether calls will be " - "classified as heterozygozous or homozygous in " - "different samples. " - "All calls will be considered homozygous if false.") + parser.add_argument("-he", "--heterozygous", default=False, action='store_true', + help="If set, all heterozygous calls will be considered homozygous.") + parser.add_argument("-nge", "--no_genotyping", default=True, action='store_false', + help="If set, calls will not be genotyped based on read depth " + "computations made by duphold.") args = parser.parse_args(in_args) return args -def obtain_sites_and_genotypes(input_fns, heterozygous=True): +def obtain_sites_and_genotypes(input_fns, heterozygous=True, genotyping=True): """ Obtain all CNV sites from a list of VCF files :param input_fns: List of paths to VCF files :param heterozygous: Boolean indicating whether calls will be classified into heterozygous or homozygous. All calls will be homozygous if false. + :param genotyping: Boolean indicating whether calls will be genotyped + based on read depth computed by duphold. :return: Set containing all sites. Sites are tuples (chrom, pos, alt, end, sv_type, inschrom, inspos). Dict containing sites as keys and genotype lines of samples as values. @@ -96,27 +98,44 @@ def obtain_sites_and_genotypes(input_fns, heterozygous=True): dhffc = record.samples[0]["DHFFC"] # set genotypes based on format gt = [".", "."] - non_calls = [(None, None), (0, 0)] + non_calls = [(".", "."), (0, 0)] called_gt = record.samples[0]["GT"] - if not heterozygous: - gt[0] = "1" - gt[1] = "1" + if not genotyping: + if called_gt[0] == None: + gt[0] = "." + else: + gt[0] = str(called_gt[0]) + if called_gt[1] == None: + gt[1] = "." + else: + gt[1] = str(called_gt[1]) else: if sv_type == "INS": # do not call INS if single sample shows reference only or no call if called_gt in non_calls: - gt[0] = called_gt[0] - gt[1] = called_gt[1] + if called_gt[0] == None: + gt[0] = "." + else: + gt[0] = str(called_gt[0]) + if called_gt[1] == None: + gt[1] = "." + else: + gt[1] = str(called_gt[1]) else: # call can be only be ./1, based on available evidence gt[1] = "1" elif sv_type == "DEL": + print("Wrong") if dhfc >= 0 and dhfc < 0.25: gt[0] = "1" gt[1] = "1" elif dhfc >= 0.25 and dhfc < 0.75: - gt[0] = "0" - gt[1] = "1" + if not heterozygous: + gt[0] = "1" + gt[1] = "1" + else: + gt[0] = "0" + gt[1] = "1" else: gt[0] = "0" gt[1] = "0" @@ -125,8 +144,12 @@ def obtain_sites_and_genotypes(input_fns, heterozygous=True): gt[0] = "0" gt[1] = "0" elif dhfc >= 1.25 and dhfc < 1.75: - gt[0] = "0" - gt[1] = "1" + if not heterozygous: + gt[0] = "1" + gt[1] = "1" + else: + gt[0] = "0" + gt[1] = "1" else: gt[0] = "1" gt[1] = "1" @@ -404,7 +427,7 @@ def write_to_output_vcf(output_vcf_fn, chrom_length_dict, sites_interval_trees, output_vcf.write("\n") return 0 -def merge_vcfs(input_fn, output_vcf_fn, chrom_length_dict, reciprocal=0.5, heterozygous=True): +def merge_vcfs(input_fn, output_vcf_fn, chrom_length_dict, reciprocal=0.5, heterozygous=True, genotyping=True): """ Merge VCF files of single samples containing CNVs called by Hecaton, producing a single VCF file containing the merged calls @@ -414,6 +437,8 @@ def merge_vcfs(input_fn, output_vcf_fn, chrom_length_dict, reciprocal=0.5, heter length of chromosomes as values :param heterozygous: Boolean indicating whether calls will be classified into heterozygous or homozygous. All calls will be homozygous if false. + :param genotyping: Boolean indicating whether calls will be genotyped + based on read depth computed by duphold. :param reciprocal: Minimum reciprocal overlap required to collapse CNVs :return: 0 (integer) """ @@ -423,7 +448,7 @@ def merge_vcfs(input_fn, output_vcf_fn, chrom_length_dict, reciprocal=0.5, heter for line in input_file: input_fns.append(line.strip()) # obtain all sites and genotype information - sites_set, samples_dict = obtain_sites_and_genotypes(input_fns, heterozygous) + sites_set, samples_dict = obtain_sites_and_genotypes(input_fns, heterozygous, genotyping) # get interval tree in which sites considered to be the same CNV are collapsed sites_interval_tree = get_interval_tree_from_sites_set(sites_set, chrom_length_dict, @@ -447,7 +472,7 @@ def main(): raise ValueError("Reciprocal overlap must be between 0 and 1") # merge vcf files merge_vcfs(args.input_file, args.output_vcf, chrom_length_dict, - args.reciprocal, args.heterozygous) + args.reciprocal, args.heterozygous, args.no_genotyping) if __name__ == "__main__": main()