Statistics (#29)

gather counts for GVF and VCF files into a report

Author: khetherin

convert_gvf_to_vcf/assistingconverter.py

@@ -69,7 +69,7 @@ def convert_gvf_attributes_to_vcf_values(column9_of_gvf,
     gvf_attribute_dictionary = get_gvf_attributes(column9_of_gvf)
     vcf_info_values = {} # key is info field value; value is value
     vcf_format_values = {} # key is format field value; value is value
-    catching_for_review = []
+    dropped_gvf_attributes = []
     for attrib_key, attrib_value in gvf_attribute_dictionary.items():
         if attrib_key in mapping_attribute_dict:
             field_values = mapping_attribute_dict[attrib_key]
@@ -96,5 +96,5 @@ def convert_gvf_attributes_to_vcf_values(column9_of_gvf,
                     logger.warning(f"Unsupported Field: {field}")
         else:
             logger.info(f"catching attribute keys for review at a later date {attrib_key} {attrib_value}")
-            catching_for_review.append(attrib_key)
+            dropped_gvf_attributes.append(attrib_key)
     return gvf_attribute_dictionary, vcf_info_values, vcf_format_values

convert_gvf_to_vcf/conversionstatistics.py

@@ -0,0 +1,137 @@
+"""
+The purpose of this file is to calculate the (conversion) statistics of the GVF and VCF file.
+"""
+import os
+from datetime import datetime
+import hashlib
+from collections import Counter
+
+class FileStatistics:
+    """
+    The responsibility of this class is to accumulate counts and calculate statistics of a GVF to VCF file conversion.
+    """
+    def __init__(self, gvf_file_path, gvf_pragmas, samples):
+        self.gvf_file_path = gvf_file_path
+        # file name
+        self.gvf_file_name = os.path.basename(self.gvf_file_path)
+        self.gvf_extension = os.path.splitext(self.gvf_file_path)[1]
+        # file metadata
+        self.gvf_metadata_stats = os.stat(self.gvf_file_path)
+        self.gvf_file_size = self.gvf_metadata_stats.st_size
+        gvf_last_modified_timestamp_raw = self.gvf_metadata_stats.st_mtime
+        self.gvf_last_modified_timestamp = datetime.fromtimestamp(gvf_last_modified_timestamp_raw)
+        self.gvf_version = self.find_version(gvf_pragmas)
+        self.gvf_md5 = self.get_file_md5(self.gvf_file_path)
+        # biological samples
+        self.sample_number = len(samples) #
+        self.vcf_sample_number_count = Counter()
+        #TODO: add missing samples in VCF - awaiting bug fix
+        #TODO: number of times a sample is missing in the merged vcf line, count per sample
+        # chromsomes
+        self.gvf_chromosome_count = Counter()
+        self.vcf_chromosome_count = Counter()
+        # # variants
+        self.gvf_feature_line_count = 0
+        self.gvf_sv_count = Counter()
+        self.vcf_data_line_count = 0 # used in property
+        self.vcf_number_of_merges = 0
+        self.vcf_alt_alleles_count = Counter() #used in property
+        # # attribute mapping
+        self.vcf_info_counter = Counter() # used in property
+        self.vcf_format_counter = Counter()
+        self.vcf_variant_region_SOID = Counter()
+        self.vcf_variant_call_SOID = Counter()
+
+    @property
+    def vcf_alt_missing(self):
+        return self.vcf_alt_alleles_count['.']
+
+    @property
+    def vcf_imprecise_variants(self):
+        return self.vcf_info_counter["IMPRECISE"]
+
+    @property
+    def vcf_precise_variants(self):
+        return self.vcf_data_line_count - self.vcf_imprecise_variants
+
+    @staticmethod
+    def get_file_md5(path_to_file):
+        hash_md5 = hashlib.md5()
+        with open(path_to_file, "rb") as f:
+            # Read file in 4KB chunks until you get to an empty byte string
+            for chunk in iter(lambda: f.read(4096), b""):
+                hash_md5.update(chunk)
+        md5 = hash_md5.hexdigest()
+        return md5
+
+    @staticmethod
+    def find_version(header_list):
+        """ Finds the gvf version in the header.
+        :params: header_list
+        :return: gvf version or if not found None
+        """
+        for line in header_list:
+            if "##gvf-version" in line:
+                try:
+                    return line.split('=')[1].strip()
+                except IndexError:
+                    return f"GVF version unknown: {line}"
+        return None
+
+    def __str__(self):
+        summary_string = (f"======GVF======\n"
+                          f"File name = {self.gvf_file_name}\n"
+                          f"File path = {self.gvf_file_path}\n"
+                          f"File extension = {self.gvf_extension}\n"
+                          f"File size = {self.gvf_file_size} bytes\n"
+                          f"Last mod Timestamp = {self.gvf_last_modified_timestamp}\n"
+                          f"Version = {self.gvf_version}\n"
+                          f"md5 string = {self.gvf_md5}\n"
+                          f"gvf chrom = {self.gvf_chromosome_count}\n"
+                          f"vcf chrom = {self.vcf_chromosome_count}\n"
+                          f"gvf_feature_line_count = {self.gvf_feature_line_count}\n"
+                          f"vcf_data_line_count = {self.vcf_data_line_count}\n"
+                          f"vcf_number_of_merges = {self.vcf_number_of_merges}\n"
+                          f"vcf_alt_allele_count = {self.vcf_alt_alleles_count}\n"
+                          f"vcf_info_counter = {self.vcf_info_counter}\n"
+                          f"vcf_format_counter = {self.vcf_format_counter}\n"
+                          f"vcf_sample_count = {self.vcf_sample_number_count}\n"
+                          f"vcf_alt_missing = {self.vcf_alt_missing}\n"
+                          f"vcf_imprecise_variants = {self.vcf_imprecise_variants}\n"
+                          f"vcf_precise_variants = {self.vcf_precise_variants}\n"
+                          f"=====\n"
+                          )
+        return summary_string
+
+    def print_report(self, file_to_write):
+        key_width = 30
+        text_report = f'''
+        GVF
+            {"File name":<{key_width}}\t{self.gvf_file_name:>}
+            {"File path":<{key_width}}\t{self.gvf_file_path:>}
+            {"File extension":<{key_width}}\t{self.gvf_extension:>}
+            {"File size":<{key_width}}\t{self.gvf_file_size:>}
+            {"Last modified":<{key_width}}\t{self.gvf_last_modified_timestamp}
+            {"Version":<{key_width}}\t{self.gvf_version:>}
+            {"md5":<{key_width}}\t{self.gvf_md5:>}
+            
+             {"COUNTS":^{key_width}}
+             {"GVF Chromosome counts:":<}\n\t\t\t\t{self.gvf_chromosome_count}
+             {"Number of GVF feature lines:":<{key_width}}\t{self.gvf_feature_line_count:>}
+             {"Number of samples:":<{key_width}}\t{self.sample_number}
+             {"Number of structural variants seen:":<{key_width}}\n\t\t\t\t{self.gvf_sv_count}
+        VCF
+            {"COUNTS":^{key_width}}
+            {"VCF Chromosome counts:":<{key_width}}\n\t\t\t\t{self.vcf_chromosome_count}
+            {"VCF ALT allele counts:":<{key_width}}\n\t\t\t\t{self.vcf_alt_alleles_count}
+            {"Number of VCF data lines:":<{key_width}}\t{self.vcf_data_line_count:>}
+            {"Number of VCF data line merges:":<{key_width}}\t{self.vcf_number_of_merges:>}
+            {"Number of times each VCF sample has been seen:":<{key_width}}\n\t\t\t\t{self.vcf_sample_number_count}
+            {"Number of times INFO keys seen"}\n\t\t\t\t{self.vcf_info_counter}
+            {"VCF variant region Sequence Ontology ID counts"}\n\t\t\t\t{self.vcf_variant_region_SOID}
+            {"VCF variant call Sequence Ontology ID counts"}\n\t\t\t\t{self.vcf_variant_call_SOID}
+            
+        '''
+        #
+        with open(file_to_write, "w") as stats_file:
+            stats_file.write(text_report)

convert_gvf_to_vcf/convertGVFtoVCF.py

@@ -1,10 +1,14 @@
 import argparse
 import os
+from collections import Counter
 from ebi_eva_common_pyutils.logger import logging_config as log_cfg
+
+from convert_gvf_to_vcf.conversionstatistics import FileStatistics
 from convert_gvf_to_vcf.lookup import Lookup
 from convert_gvf_to_vcf.utils import read_in_gvf_header, read_in_gvf_data
 from convert_gvf_to_vcf.vcfline import VcfLineBuilder
 
+
 logger = log_cfg.get_logger(__name__)
 
 # setting up paths to useful directories
@@ -236,7 +240,7 @@ def get_pragma_tokens(pragma_value, first_delimiter, second_delimiter):
     return pragma_tokens
 
 
-def write_header(vcf_output, pragmas_for_vcf, header_lines_per_type, is_missing_format_value, samples):
+def write_header(vcf_output, pragmas_for_vcf, header_lines_per_type, header_fields, samples):
     logger.info(f"Total number of samples in this VCF: {len(samples)}")
     vcf_header_file = vcf_output + '_header'
     with open(vcf_header_file, "w") as vcf_header_output:
@@ -250,7 +254,6 @@ def write_header(vcf_output, pragmas_for_vcf, header_lines_per_type, is_missing_
 
         # Part 2 of VCF file: Write the VCF header line.
         # Write the header.
-        header_fields = generate_vcf_header_line(is_missing_format=is_missing_format_value, samples=samples)
         vcf_header_output.write(f"{header_fields}\n")
     return vcf_header_file
 
@@ -259,14 +262,13 @@ def convert(gvf_input, vcf_output, assembly):
     logger.info("Running the GVF to VCF converter")
     logger.info(f"The provided input file is: {gvf_input}")
     logger.info(f"The provided output file is: {vcf_output}")
+
     if assembly:
         logger.info(f"The provided assembly file is: {assembly}")
     assembly_file = os.path.abspath(assembly)
     assert os.path.isfile(assembly_file), "Assembly file does not exist"
+    assert os.path.isfile(gvf_input), "GVF file does not exist"
 
-    # Read input file and separate out its components
-    logger.info(f"Reading in the following GVF header from {gvf_input}")
-    gvf_pragmas, gvf_pragma_comments = read_in_gvf_header(gvf_input)
     # Creating lookup object to store important dictionaries and log what has been stored.
     reference_lookup = Lookup(assembly_file)
     logger.info("Creating the reference lookup object.")
@@ -275,51 +277,80 @@ def convert(gvf_input, vcf_output, assembly):
     logger.info(f"Storing the assembly file: {assembly_file}")
     logger.info("Storing the IUPAC ambiguity dictionary.")
 
+    # Read input file and separate out its components
+    logger.info(f"Reading in the following GVF header from {gvf_input}")
+    gvf_pragmas, gvf_pragma_comments = read_in_gvf_header(gvf_input)
+
     # Preparation work:
     # Store the VCF metainformation and ensure preservation of important GVF data.
     # This information will be useful when creating the VCF header.
-    # TODO: refactor function generate_vcf_metainfo
     (
         pragmas_for_vcf,
         samples
     ) = convert_gvf_pragmas_for_vcf_header(gvf_pragmas, gvf_pragma_comments, reference_lookup)
-
-    # TODO: place the all_header_lines_per_type_dict into the reference_lookup.
+    report = FileStatistics(gvf_file_path=gvf_input, gvf_pragmas=gvf_pragmas, samples=samples)
 
     # Create data structure to store all possible outcomes for header lines (for fields ALT, INFO, FILTER, FORMAT)
     all_header_lines_per_type_dict = {
         htype: generate_vcf_header_structured_lines(htype, reference_lookup.mapping_attribute_dict) for htype in
         ["ALT", "INFO", "FILTER", "FORMAT"]
     }
+
     vcf_builder = VcfLineBuilder(all_header_lines_per_type_dict, reference_lookup, samples)
     is_missing_format_value = True
+    # VCF dataline generation
     # Convert each feature line in the GVF file to a VCF object (stores all the data for a line in the VCF file).
     # NOTE: Main Logic lives here.
     vcf_data_file = vcf_output + '_data_lines'
     with open(vcf_data_file, "w") as open_data_lines:
         logger.info("Generating the VCF datalines")
         previous_vcf_line = None
-        for gvf_lines_obj in read_in_gvf_data(gvf_input):
-            current_vcf_line = vcf_builder.build_vcf_line(gvf_lines_obj)
+        for gvf_entry in read_in_gvf_data(gvf_input):
+            # record GVF counts
+            report.gvf_feature_line_count += 1
+            report.gvf_chromosome_count[gvf_entry.seqid] += 1
+            report.gvf_sv_count.update([gvf_entry.feature_type])
+            # create the VCF line object
+            current_vcf_line = vcf_builder.build_vcf_line(gvf_entry)
             # is_missing_format_value will only be true if all the format field are missing.
             is_missing_format_value = is_missing_format_value and current_vcf_line.format_keys == ['.']
             # Each GVF feature has been converted to a VCF object so begin comparing and merging the VCF objects.
             if previous_vcf_line:
                 if current_vcf_line == previous_vcf_line:
                     current_vcf_line.merge(previous_vcf_line, list_of_sample_names=samples)
+                    report.vcf_number_of_merges += 1
                 else:
-                    open_data_lines.write(str(previous_vcf_line) + "\n")
+                    # TODO: address this in next bug fix
+                    assert current_vcf_line > previous_vcf_line, f"File not sorted.\ncurrent_vcf_line.pos {current_vcf_line.pos} is smaller than previous_vcf_line.pos {previous_vcf_line.pos}. See the following line:\n{str(gvf_entry)}"
+                    record_vcf_entry(open_data_lines, previous_vcf_line, report)
             previous_vcf_line = current_vcf_line
+        # Process the final previous_vcf_line
         if previous_vcf_line:
-            open_data_lines.write(str(previous_vcf_line) + "\n")
+            record_vcf_entry(open_data_lines, previous_vcf_line, report)
         else:
             logger.warning("No feature lines were found for this GVF file.")
 
+    # VCF header generation
     header_lines_per_type = vcf_builder.build_vcf_header()
-    vcf_header_file = write_header(vcf_output, pragmas_for_vcf, header_lines_per_type,
-                                   is_missing_format_value, samples)
+    header_fields = generate_vcf_header_line(is_missing_format=is_missing_format_value, samples=samples)
+    vcf_header_file = write_header(vcf_output, pragmas_for_vcf, header_lines_per_type, header_fields, samples)
 
     logger.info(f"Combining the header and data lines to the following VCF output: {vcf_output}")
+    construct_vcf_output(vcf_header_file, vcf_data_file, vcf_output)
+
+    logger.info("Remove the temporary files")
+    cleanup_temp_files([vcf_header_file, vcf_data_file])
+
+    logger.info("Printing the summary of conversion report.")
+    vcf_output_directory = os.path.dirname(vcf_output)
+    stats_summary_file = os.path.join(vcf_output_directory, "summary_stats.txt")
+    report.print_report(stats_summary_file)
+
+    logger.info("GVF to VCF conversion complete")
+
+
+#helper functions for convert
+def construct_vcf_output(vcf_header_file, vcf_data_file, vcf_output):
     with open(vcf_output, "w") as vcf_output:
         with open(vcf_header_file, "r") as vcf_header_fh:
             for line in vcf_header_fh:
@@ -328,12 +359,23 @@ def convert(gvf_input, vcf_output, assembly):
             for line in vcf_data_fh:
                 vcf_output.write(line)
     vcf_output.close()
-    logger.info("Remove the temporary files")
-    if os.path.exists(vcf_header_file):
-        os.remove(vcf_header_file)
-    if os.path.exists(vcf_data_file):
-        os.remove(vcf_data_file)
-    logger.info("GVF to VCF conversion complete")
+
+def record_vcf_entry(open_data_lines, previous_vcf_line, report):
+    # write the VCF line and record counts (use update when iterable)
+    open_data_lines.write(str(previous_vcf_line) + "\n")
+    report.vcf_data_line_count += 1
+    report.vcf_chromosome_count[previous_vcf_line.chrom] += 1
+    report.vcf_alt_alleles_count[previous_vcf_line.alt] += 1
+    report.vcf_info_counter.update(previous_vcf_line.info_dict.keys())
+    report.vcf_variant_region_SOID.update([previous_vcf_line.info_dict.get("VARREGSOID")])
+    report.vcf_variant_call_SOID.update([previous_vcf_line.info_dict.get("VARCALLSOID")])
+    report.vcf_sample_number_count.update(previous_vcf_line.order_sample_names)
+    report.vcf_format_counter.update(previous_vcf_line.vcf_values_for_format)
+
+def cleanup_temp_files(list_of_temp_files):
+    for temp_file in list_of_temp_files:
+        if os.path.exists(temp_file):
+            os.remove(temp_file)
 
 def main():
     # Parse command line arguments

convert_gvf_to_vcf/vcfline.py

@@ -479,6 +479,16 @@ def __eq__(self, other_vcf_line):
             return (self.chrom == other_vcf_line.chrom) and (self.pos == other_vcf_line.pos) and (self.ref == other_vcf_line.ref)
         return False
 
+    def __gt__(self, other_vcf_line):
+        """ This checks ordering by pos within the SAME chromosome for the merging algorithm.
+        Order is not enforced between different chromosomes to avoid cross-chromosome ordering.
+        :param: other_vcf_line: VcfLine object to compare against.
+        :return: boolean - True if position comes after the previous vcf line on the same chromosome.
+        """
+        if isinstance(other_vcf_line, VcfLine):
+            return (self.chrom != other_vcf_line.chrom) or (self.pos > other_vcf_line.pos)
+        return False
+
     @property
     def format_keys(self):
         set_of_format_key = set([format_key