<commit_before>#!/usr/bin/env perl
# Copyright [1999-2015] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute
# 
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# 
#      http://www.apache.org/licenses/LICENSE-2.0
# 
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.



=head1 CONTACT

  Please email comments or questions to the public Ensembl
  developers list at <http://lists.ensembl.org/mailman/listinfo/dev>.

  Questions may also be sent to the Ensembl help desk at
  <helpdesk.org>.

=cut


=head1 NAME
import_vcf.pl - imports variations from a VCF file into an Ensembl variation DB

by <NAME> (<EMAIL>)
=cut

use strict;
use Bio::EnsEMBL::Registry;
use Bio::EnsEMBL::Utils::Sequence qw(reverse_comp);
use Bio::EnsEMBL::Variation::Utils::VEP qw(parse_line get_all_consequences);
use Bio::EnsEMBL::Variation::Utils::Sequence qw(SO_variation_class);

# object types need to imported explicitly to use new_fast
use Bio::EnsEMBL::Variation::Variation;
use Bio::EnsEMBL::Variation::IndividualGenotype;

# use this for remapping
use Bio::EnsEMBL::SimpleFeature;

use Getopt::Long;
use FileHandle;
use Socket;
use IO::Handle;
use Data::Dumper;
use Time::HiRes qw(gettimeofday tv_interval);
use ImportUtils qw(load);
use Digest::MD5 qw(md5_hex);
use Cwd 'abs_path';

use constant DISTANCE => 100_000;
use constant MAX_SHORT => 2**16 -1;

my %Printable = ( "\\"=>'\\', "\r"=>'r', "\n"=>'n', "\t"=>'t', "\""=>'"' );

$| = 1;

my $config = configure();

# open cross-process pipes for fork comms
socketpair(CHILD, PARENT, AF_UNIX, SOCK_STREAM, PF_UNSPEC) or die "ERROR: Failed to open socketpair: $!";
CHILD->autoflush(1);
PARENT->autoflush(1);

$config->{start_time} = [gettimeofday];

# check if we are forking
if(defined($config->{fork})) {
	die "ERROR: Could not find input file\n" unless defined($config->{input_file}) && -e $config->{input_file};
	die "ERROR: Input file is not bgzipped, cannot fork\n" unless $config->{input_file} =~ /\.gz$/;
	die "ERROR: Tabix index file ".$config->{input_file}.".tbi not found, cannot fork\n" unless $config->{input_file}.'.tbi';
	
	run_forks($config);
}
else {
	main($config);
}

my $elapsed = tv_interval($config->{start_time}, [gettimeofday]);

debug($config, "Took $elapsed s");

sub configure {
	
	# COMMAND LINE OPTIONS
	######################
	
	# get command-line options
	my $config = {};
	my $args = scalar @ARGV;
	
	GetOptions(
		$config,
		
		'help|h',
		'input_file|i=s',
		'tmpdir=s',
		'tmpfile=s',
		'config=s',
		
		'progress_update=i',
		'no_progress',
		
		'species=s',
		'registry|r=s',
		'host=s',
		'database|db=s',
		'user=s',
		'password=s',
		'port=i',
		
		'sql=s',
		'coord_system=s',
		
		'source=s',
		'source_desc=s',
		'population|pop=s',
		'pedigree=s',
		'panel=s',
		'gmaf=s',
		'somatic',
		
		'flank=s',
		'gp',
    'remap=s',
		'ind_prefix=s',
		'pop_prefix=s',
		'var_prefix=s',
		
		'disable_keys',
		'tables=s',
		'skip_tables=s',
		'add_tables=s',
		
		'only_existing',
    'no_merge',
		'skip_n',
		'mart_genotypes',
		
		'create_name',
		'chrom_regexp=s',
		'force_no_var',
    'ss_ids',
		
		'fork=i',
		'test=i',
		'backup',
		'move',
		
		'recover',
		'recover_point=s',
		'recover_pos=s',
		'no_recover',
		
		'cache=s',
		'fasta=s',
		
	# die if we can't parse arguments - better to get user to sort out their command line
	# than potentially do the wrong thing
	) or die "ERROR: Failed to parse command line arguments - check the documentation!\n";
	
	
	# print usage message if requested or no args supplied
	if(defined($config->{help}) || !$args) {
		&usage;
		exit(0);
	}
	
	# read config from file?
  read_config_from_file($config, $config->{config}) if defined $config->{config};
	
	# remove recover from config otherwise changes session id
	my $recover;
	if(defined($config->{recover})) {
		$recover = 1;
		delete $config->{recover};
	}
	
	# create session ID from config
	$config->{session_id} = md5_hex(join '|', (map {$_." ".$config->{$_}} sort keys %$config));
	
	debug($config, "Session ID is ", $config->{session_id});
	
	$config->{recover} = 1 if $recover;
	
	# sanity checks
	die("ERROR: Cannot run in test mode using forks\n") if defined($config->{fork}) and defined($config->{test});
	die("ERROR: Cannot manually recover using forks\n") if defined($config->{fork}) and (defined($config->{recover_point}) || defined($config->{recover_pos}));
	
	# recover?
	if(defined $config->{recover}) {
		open IN, (join '/', ($ENV{HOME}, '.import_vcf', $config->{session_id})) or die "ERROR: Could not recover session with ID ".$config->{session_id}."\n";
		my $a = <IN>;
		($config->{recover_point}, $config->{pid}) = split ' ', $a;
		close IN;
		
		die "ERROR: Cannot recover - session with ID ".$config->{session_id}." has been flagged as finished" if $config->{recover_point} eq 'FINISHED';
		
		debug($config, "Found recover point ", $config->{recover_point}, " for session ID ", $config->{session_id}, " PID ", $config->{pid});
	}
	
	# set defaults
	$config->{species}         ||= "human";
	$config->{flank}           ||= 200;
	$config->{port}            ||= 3306;
	$config->{format}            = 'vcf';
	$config->{ind_prefix}      ||= '';
	$config->{pop_prefix}      ||= '';
	$config->{coord_system}    ||= 'chromosome';
	$config->{progress_update} ||= 100;
	$config->{pid}             ||= $$;
	$config->{somatic}         ||= 0;
  
  # check remap arg
  if(defined($config->{remap})) {
    die "ERROR: remap argument incorrect - must be of format --remap [from_assembly],[to_assembly]" unless $config->{remap} =~ /^\S+?\,\S+?$/;
    
    $config->{remap} = [split(',', $config->{remap})];
  }
	
	# recovery not possible if forking
	#$config->{no_recover} = 1 if defined($config->{fork});
	
	# set default list of tables to write to
	my $tables = {
		'variation'                       => 1,
		'variation_feature'               => 1,
		'variation_synonym'               => 1,
		'flanking_sequence'               => 1,
		'allele'                          => 1,
		'population_genotype'             => 1,
		'compressed_genotype_var'         => 1,
		'individual_genotype_multiple_bp' => 0,
		'compressed_genotype_region'      => 0,
		'transcript_variation'            => 0,
		'sample'                          => 1,
		'population'                      => 1,
		'individual'                      => 1,
		'individual_population'           => 1,
		'allele_code'                     => 1,
		'genotype_code'                   => 1,
		'meta_coord'                      => 1,
		'source'                          => 1,
	};
	
	# override this with options if provided
	if(defined($config->{tables})) {
		
		# reset
		$tables->{$_} = 0 foreach keys %$tables;
		
		# set include tables
		foreach my $table(split /\,/, $config->{tables}) {
			$tables->{$table} = 1 if defined($tables->{$table});
		}
	}
	
	if(defined($config->{add_tables})) {
		
		# add tables
		foreach my $table(split /\,/, $config->{add_tables}) {
			$tables->{$table} = 1 if defined($tables->{$table});
		}
	}
	
	if(defined($config->{skip_tables})) {
		
		# set skip tables
		foreach my $table(split /\,/, $config->{skip_tables}) {
			$tables->{$table} = 0 if defined($tables->{$table});
		}
	}
	
	# force some back in
	$tables->{$_} = 1 for qw/source meta_coord/;
	
	$tables->{individual_genotype_multiple_bp} = 1 if defined($config->{mart_genotypes});
	
	# special case for sample, we also want to set the other sample tables to on
	if($tables->{sample}) {
		$tables->{$_} = 1 for qw/population individual individual_population/;
	}
	
	# force population if user wants allele or population_genotype
	if($tables->{allele} || $tables->{population_genotype}) {
		$tables->{$_} = 1 for qw/population sample/;
	}
	
	# force individual tables for individual level data
	if($tables->{compressed_genotype_region} || $tables->{individual_genotype_multiple_bp} || $tables->{compressed_genotype_var}) {
		$tables->{$_} = 1 for qw/population individual individual_population sample variation variation_synonym variation_feature flanking_sequence/;
	}
	
	if($tables->{population_genotype} || $tables->{compressed_genotype_region} || $tables->{compressed_genotype_var}) {
		$tables->{$_} = 1 for qw/genotype_code variation variation_synonym variation_feature flanking_sequence/;
	}
	
	if($tables->{allele} || $tables->{genotype_code}) {
		$tables->{$_} = 1 for qw/allele_code variation variation_synonym variation_feature flanking_sequence/;
	}
	
	# won't be writing to these tables if only_existing mode
	if(defined $config->{only_existing}) {
		$tables->{$_} = 0 for qw/source variation variation_synonym variation_feature flanking_sequence/;
	}
	
	# check that at least one has been set
	die "ERROR: no tables left included\n" unless grep {$tables->{$_}} keys %$tables;
	
	$config->{tables} = $tables;
	
	die "ERROR: tmpdir not specified\n" if !defined $config->{tmpdir} && $tables->{compressed_genotype_region};
	$config->{tmpfile} ||= 'compress.txt';
	$ImportUtils::TMP_DIR  = $config->{tmpdir};
	$ImportUtils::TMP_FILE = $config->{tmpfile};
	
	
	$config->{reg} = 'Bio::EnsEMBL::Registry';
	
	# VEP stuff
	$config->{vep}->{chunk_size}        ||= 50000;
	$config->{vep}->{cache_region_size} ||= 1000000;
	$config->{vep}->{compress}          ||= 'zcat';
	$config->{vep}->{terms}               = 'SO';
	$config->{vep}->{tr_cache}            = {};
	$config->{vep}->{rf_cache}            = {};
	$config->{vep}->{quiet}               = 1;
	$config->{vep}->{original}            = 1;
	$config->{vep}->{no_progress}         = 1;
	$config->{buffer_size}              ||= 100;
	$config->{vep}->{species}             = $config->{species};
	
	if(defined($config->{cache})) {
		die("ERROR: Could not find cache directory ".$config->{cache}) unless -e $config->{cache};		
		$config->{vep}->{dir}             = $config->{cache};
		$config->{vep}->{cache}           = 1;
		$config->{vep}->{offline}         = 1;
	}
	
	
    if(defined($config->{fasta})) {
        die "ERROR: Specified FASTA file/directory not found" unless -e $config->{fasta};
        
        eval q{ use Bio::DB::Fasta; };
        
        if($@) {
            die("ERROR: Could not load required BioPerl module\n");
        }
        
        # try to overwrite sequence method in Slice
        eval q{
            package Bio::EnsEMBL::Slice;
            
            # define a global variable so that we can pull in config hash
            our $config;
            
            {
                # don't want a redefine warning spat out, thanks
                no warnings 'redefine';
                
                # overwrite seq method to read from FASTA DB
                sub seq {
                    my $self = shift;
                    
                    # special case for in-between (insert) coordinates
                    return '' if($self->start() == $self->end() + 1);
                    
                    my $seq;
                    
                    if(defined($config->{fasta_db})) {
                        $seq = $config->{fasta_db}->seq($self->seq_region_name, $self->start => $self->end);
                        reverse_comp(\$seq) if $self->strand < 0;
                    }
                    
                    else {
                        return $self->{'seq'} if($self->{'seq'});
                      
                        if($self->adaptor()) {
                          my $seqAdaptor = $self->adaptor()->db()->get_SequenceAdaptor();
                          return ${$seqAdaptor->fetch_by_Slice_start_end_strand($self,1,undef,1)};
                        }
                    }
                    
                    # default to a string of Ns if we couldn't get sequence
                    $seq ||= 'N' x $self->length();
                    
                    return $seq;
                }
            }
            
            1;
        };
        
        if($@) {
            die("ERROR: Could not redefine sequence method\n");
        }
        
        # copy to Slice for offline sequence fetching
        $Bio::EnsEMBL::Slice::config = $config->{vep};
        
        # spoof a coordinate system
        $config->{vep}->{coord_system} = Bio::EnsEMBL::CoordSystem->new(
            -NAME => 'chromosome',
            -RANK => 1,
        );
        
        debug($config, "Checking/creating FASTA index");
        $config->{vep}->{fasta_db} = Bio::DB::Fasta->new($config->{fasta});
    }
	
	# get terminal width for progress bars
	my $width;
	
	# module may not be installed
	eval q{
		use Term::ReadKey;
	};
	
	if(!$@) {
		my ($w, $h);
		
		# module may be installed, but e.g.
		eval {
			($w, $h) = GetTerminalSize();
		};
		
		$width = $w if defined $w;
	}
	
	$width ||= 60;
	$width -= 12;
	$config->{terminal_width} = $width;
	
	return $config;
}

# reads config from a file
sub read_config_from_file {
    my $config = shift;
    my $file = shift;
    
    open CONFIG, $file or die "ERROR: Could not open config file \"$file\"\n";
    
    while(<CONFIG>) {
        next if /^\#/;
        my @split = split /\s+|\=/;
        my $key = shift @split;
        $key =~ s/^\-//g;
        
        if(defined($config->{$key}) && ref($config->{$key}) eq 'ARRAY') {
            push @{$config->{$key}}, @split;
        }
        else {
            $config->{$key} ||= $split[0];
        }
    }
    
    close CONFIG;
    
    debug($config, "Read configuration from $file") unless defined($config->{quiet});
}


# main sub-routine does most of the code execution
sub main {
	
	my $config = shift;
	
	# log start time
	my $start_time = time();
	
	
	# SET UP VARIABLES
	##################
	
	my (
		%headers,
		$prev_seq_region,
		$genotypes,
		$var_counter,
		@vf_buffer,
	);
	
	
	# CONNECT TO DBS
	################
	
	connect_to_dbs($config);
	
	# populate from SQL?
	if(defined($config->{sql}) && !defined($config->{recover})) {
		sql_populate($config, $config->{sql});
		
		if(defined($config->{test})) {
			debug($config, "(TEST) Adding schema version to meta table");
		}
		else {
			# add schema version to meta
			my $sth = $config->{dbVar}->prepare(qq{
				INSERT INTO meta (species_id, meta_key, meta_value)
				VALUES (NULL, ?, ?)
			});
			$sth->execute('schema_version', $config->{reg}->software_version);
			$sth->finish();
		}
		
		# try and do attrib
		attrib($config);
	}
	
	# get adaptors
	get_adaptors($config);
	
	# insert from files if recovering
	if(defined($config->{recover})) {
		debug($config, "Importing data from recovered session's temporary files");
		
		foreach my $table(qw(allele population_genotype compressed_genotype_region)) {
			debug($config, "Importing data into $table");
			import_tmp_file($config, $table);
		}
	}
	
	# backup
	backup($config) if defined($config->{backup}) || defined($config->{move});
	
	# TEST MODE?
	############
	
	debug($config, "Running in test mode - will read first ", $config->{test}, " lines of input") if defined($config->{test});
	
	# DB PREP
	#########
	
	# get seq_region_id hash
	$config->{seq_region_ids} = get_seq_region_ids($config);
	
	# if failed, try and copy from core DB
	if(!scalar keys %{$config->{seq_region_ids}}) {
		
		copy_seq_region_from_core($config);
		
		# now reload
		$config->{seq_region_ids} = get_seq_region_ids($config);
	}
	
	die("ERROR: seq_region not populated\n") unless scalar keys %{$config->{seq_region_ids}};
	
	# get/set source_id
	die("ERROR: no source specified\n") if !(defined $config->{source}) && !defined($config->{only_existing});
	$config->{source_id} = get_source_id($config) unless defined($config->{only_existing});
	
	# get population object
	if($config->{tables}->{population}) {
		die("ERROR: no population specified\n") unless defined $config->{population} || defined $config->{panel};
		$config->{populations} = population($config);
	}
	
	# disable keys if requested
	if(defined $config->{disable_keys}) {
		debug($config, "Disabling keys");
		
		foreach my $table(qw(allele population_genotype)) {#grep {$config->{tables}->{$_}} keys %$config->{tables}) {
			$config->{dbVar}->do(qq{ALTER TABLE $table DISABLE KEYS;});
		}
	}
	
	# GET INPUT FILE HANDLE
	#######################
	
	# if we forked, already have file handle but need to reopen original file
	# to read data from the column definition headers
	if(defined($config->{forked})) {
		my $tmp_file_handle = get_input_file_handle($config);
		
		while(<$tmp_file_handle>) {
			chomp;
			next if /^##/;
			
			my @split = split /\t/;
			
			# column definition line
			if(/^#/) {
				%headers = %{parse_header($config, \@split)};
				last;
			}
		}
	}
	
	else {
		$config->{in_file_handle} = get_input_file_handle($config);
	}
	
	
	# PEDIGREE FILE
	###############
	
	$config->{pedigree} = pedigree($config) if defined($config->{pedigree});
	
	# MAIN FILE LOOP
	################
	
	my $in_file_handle = $config->{in_file_handle};
	
	my $last_skipped = 0;
	
	# read the file
	while(<$in_file_handle>) {
		chomp;
		
		# header lines
		next if /^##/;
		
		my @split = split /\s+/;
		my $data = {};
		$data->{line} = $_;
		
		# column definition line
		if(/^#/) {
			%headers = %{parse_header($config, \@split)};
			
			# leave a file telling the master process to fork the others
			if(defined($config->{forked})) {
				open TMP, '> '.$ENV{HOME}.'/.import_vcf/'.$config->{pid};
				print TMP '1';
				close TMP;
			}
		}
		
		# data
		else {
			
			$config->{prev_time} ||= [gettimeofday];
			
			# recover?
			if(defined($config->{recover_point}) || defined($config->{recover_pos})) {
				$config->{skipped}->{already_processed}++;
				
				if(defined($config->{recover_point}) && md5_hex($data->{line}) eq $config->{recover_point}) {
					delete $config->{recover_point};
					$config->{recover_check} = 1;
					debug(
						$config,
						"Found recovery point, skipped ",
						$config->{skipped}->{already_processed},
						" already processed variants"
					);
				}
				
				if(defined($config->{recover_pos})) {
					my $rpos = $config->{recover_pos};
					
					if($rpos =~ /\:/) {
						my ($chr, $pos) = split /\:/, $rpos;
						if($data->{line} =~ /^$chr\s+$pos\s+/) {
							delete $config->{recover_pos};
							$config->{recover_check} = 1;
							debug(
								$config,
								"Found recovery point, skipped ",
								$config->{skipped}->{already_processed},
								" already processed variants"
							);
						}
					}
					elsif($data->{line} =~ /^\w+?\s+$rpos\s+/) {
						delete $config->{recover_pos};
						$config->{recover_check} = 1;
						debug(
							$config,
							"Found recovery point, skipped ",
							$config->{skipped}->{already_processed},
							" already processed variants"
						);
					}
				}
				
				next;
			}
			
			# check we're not skipping loads in a row
			if($last_skipped > 100 && $last_skipped =~ /(5|0)00$/) {
				debug($config, "WARNING: Skipped last $last_skipped variants, are you sure this is running OK? Maybe --gp is enabled when it shouldn't be, or vice versa?");
			}
			
			# parse into a hash
			$data->{$_} = $split[$headers{$_}] for keys %headers;
			
			# skip non-variant lines
			if($data->{ALT} eq '.') {
				$config->{skipped}->{non_variant}++;
				next;
			}
			
			# parse info column
			my %info;	
			foreach my $chunk(split /\;/, $data->{INFO}) {
				my ($key, $val) = split /\=/, $chunk;
				$info{$key} = $val;
			}
			
			$data->{info} = \%info;
			
			# skip unwanted chromosomes
			#next if defined($config->{chrom_regexp}) && $data->{'#CHROM'} !~ m/$chrom_regexp/;
			
			# use VEP's parse_line to get a skeleton VF
			($data->{tmp_vf}) = @{parse_line($config, $data->{line})};
			
			if(!defined($data->{tmp_vf})) {
				$config->{skipped}->{could_not_parse}++;
				$last_skipped++;
				next;
			}
      
      # remap?
      if(defined($config->{remap})) {
        my $success = remap($config, $data->{tmp_vf});
        if(!$success) {
          $config->{skipped}->{unable_to_remap}++;
          $last_skipped++;
          next;
        }
      }
			
			if(!defined($config->{seq_region_ids}->{$data->{tmp_vf}->{chr}})) {
				$config->{skipped}->{missing_seq_region}++;
				$last_skipped++;
				next;
			}
			
			# copy seq region ID
			$data->{tmp_vf}->{seq_region_id} = $config->{seq_region_ids}->{$data->{tmp_vf}->{chr}};
			
			# could be a structural variation feature
			next unless $data->{tmp_vf}->isa('Bio::EnsEMBL::Variation::VariationFeature');
      
      # ssIDs as IDs?
      if(defined($config->{ss_ids})) {
        my ($ss_id) = grep {$_ =~ /^\d+$/} split(/\;/, $data->{ID});
        
        if(defined($ss_id)) {
          $data->{SS_ID} = $ss_id;
        }
      }
			
			# sometimes ID has many IDs separated by ";", take the lowest rs number, otherwise the first
			if($data->{ID} =~ /\;/) {
				if($data->{ID} =~ /rs/) {
					$data->{ID} = (sort {(split("rs", $a))[-1] <=> (split("rs", $b))[-1]} split /\;/, $data->{ID})[0];
				}
        elsif(defined($config->{ss_ids}) && defined($data->{SS_ID})) {
          $data->{ID} = 'ss'.$data->{SS_ID};
        }
				else {
					$data->{ID} = (split /\;/, $data->{ID})[0];
				}
			}
      elsif(defined($config->{ss_ids}) && defined($data->{SS_ID})) {
        $data->{ID} = 'ss'.$data->{SS_ID};
      }
			
			# make a var name if none exists
			if(!defined($data->{ID}) || $data->{ID} eq '.' || defined($config->{create_name})) {
				$data->{ID} =
					($config->{var_prefix} ? $config->{var_prefix} : 'tmp').
					'_'.$data->{'#CHROM'}.'_'.$data->{POS}.'_'.$data->{REF}.'_'.$data->{ALT};
				$data->{made_up_name} = 1;
			}
			
			$data->{tmp_vf}->{variation_name} = $data->{ID};
			
			# parse genotypes
			$data->{genotypes} = get_genotypes($config, $data, \@split);
			
			# get variation object
			$data->{variation} = variation($config, $data);
			
			# transcript variation (get cons)
			get_all_consequences($config->{vep}, [$data->{tmp_vf}]) if $config->{tables}->{transcript_variation};
			
			# get variation_feature object
			$data->{vf} = variation_feature($config, $data);
			
			# attach variation to genotypes
			$_->{variation} = $data->{variation} for @{$data->{genotypes}};
			
			# skip variation if no dbID
			if(!defined($data->{variation}->{dbID}) && !defined($config->{test})) {
				$config->{skipped}->{var_not_present}++;
				$last_skipped++;
				next;
			}
			
			# transcript variation (write to DB)
			transcript_variation($config, [$data->{tmp_vf}]) if $config->{tables}->{transcript_variation} && defined($data->{tmp_vf}->dbID);
			
			#if($config->{tables}->{transcript_variation}) {
			#	push @vf_buffer, $data->{vf};
			#	if(scalar @vf_buffer == $config->{buffer_size}) {
			#		transcript_variation($config, \@vf_buffer);
			#		@vf_buffer = ();
			#	}
			#}
			
			# alleles
			allele($config, $data) if $config->{tables}->{allele};
			
			# population genotypes
			population_genotype($config, $data) if $config->{tables}->{population_genotype};
			
			# individual genotypes
			individual_genotype($config, $data) if $config->{tables}->{compressed_genotype_var} || defined($config->{mart_genotypes});
			
			# GENOTYPES BY REGION
			#####################
			
			# multi bp
			#if($config->{tables}->{individual_genotype_multiple_bp} && $force_multi && @{$data->{genotypes}}) {
			#	&multi_bp_genotype($dbVar, $data);
			#}
			
			# compressed by region
			if($config->{tables}->{compressed_genotype_region} && @{$data->{genotypes}} && !defined($config->{test})) {
				my $vf = $data->{vf};
				
				$vf->{seq_region_id} = $vf->slice->get_seq_region_id if !defined($vf->{seq_region_id});
				
				if(defined($vf->{seq_region_id}) && defined($vf->{start})) {
				
					foreach my $gt(@{$data->{genotypes}}) {
						my $sample_id = $gt->individual->dbID;
						
						next if $gt->genotype_string =~ /\./;
						
						# add to compress hash for writing later
						if (!defined $genotypes->{$sample_id}->{region_start}){
							$genotypes->{$sample_id}->{region_start} = $vf->{start};
							$genotypes->{$sample_id}->{region_end} = $vf->{end};
						}
						
						# write previous data?
						#compare with the beginning of the region if it is within the DISTANCE of compression
						if (
							defined($genotypes->{$sample_id}->{genotypes}) &&
							(
								(abs($genotypes->{$sample_id}->{region_start} - $vf->{start}) > DISTANCE()) ||
								(abs($vf->{start} - $genotypes->{$sample_id}->{region_end}) > MAX_SHORT) ||
								(defined($prev_seq_region) && $vf->{seq_region_id} != $prev_seq_region) ||
								($vf->{start} - $genotypes->{$sample_id}->{region_end} - 1 < 0)
							)
						) {
							#snp outside the region, print the region for the sample we have already visited and start a new one
							print_file($config,$genotypes, $prev_seq_region, $sample_id);
							delete $genotypes->{$sample_id}; #and remove the printed entry
							$genotypes->{$sample_id}->{region_start} = $vf->{start};
						}
						
						if ($vf->{start} != $genotypes->{$sample_id}->{region_start}){
							#compress information
							my $blob = pack ("w",$vf->{start} - $genotypes->{$sample_id}->{region_end} - 1);
							$genotypes->{$sample_id}->{genotypes} .=
								escape($blob).
								escape(pack("w", $data->{variation}->dbID || 0)).
								escape(pack("w", $config->{individualgenotype_adaptor}->_genotype_code($gt->genotype)));
						}
						else{
							#first genotype starts in the region_start, not necessary the number
							$genotypes->{$sample_id}->{genotypes} =
								escape(pack("w", $data->{variation}->dbID || 0)).
								escape(pack("w", $config->{individualgenotype_adaptor}->_genotype_code($gt->genotype)));
						}
						
						$genotypes->{$sample_id}->{region_end} = $vf->{start};
					}
				}
			}
			
			$prev_seq_region = $data->{vf}->{seq_region_id};
			$last_skipped = 0;
			
			$var_counter++;<commit_msg>Considers genotype phasing when creating compressed_genotype_region entries<commit_after>#!/usr/bin/env perl
# Copyright [1999-2015] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute
# 
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# 
#      http://www.apache.org/licenses/LICENSE-2.0
# 
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.



=head1 CONTACT

  Please email comments or questions to the public Ensembl
  developers list at <http://lists.ensembl.org/mailman/listinfo/dev>.

  Questions may also be sent to the Ensembl help desk at
  <helpdesk.org>.

=cut


=head1 NAME
import_vcf.pl - imports variations from a VCF file into an Ensembl variation DB

by <NAME> (<EMAIL>)
=cut

use strict;
use Bio::EnsEMBL::Registry;
use Bio::EnsEMBL::Utils::Sequence qw(reverse_comp);
use Bio::EnsEMBL::Variation::Utils::VEP qw(parse_line get_all_consequences);
use Bio::EnsEMBL::Variation::Utils::Sequence qw(SO_variation_class);

# object types need to imported explicitly to use new_fast
use Bio::EnsEMBL::Variation::Variation;
use Bio::EnsEMBL::Variation::IndividualGenotype;

# use this for remapping
use Bio::EnsEMBL::SimpleFeature;

use Getopt::Long;
use FileHandle;
use Socket;
use IO::Handle;
use Data::Dumper;
use Time::HiRes qw(gettimeofday tv_interval);
use ImportUtils qw(load);
use Digest::MD5 qw(md5_hex);
use Cwd 'abs_path';

use constant DISTANCE => 100_000;
use constant MAX_SHORT => 2**16 -1;

my %Printable = ( "\\"=>'\\', "\r"=>'r', "\n"=>'n', "\t"=>'t', "\""=>'"' );

$| = 1;

my $config = configure();

# open cross-process pipes for fork comms
socketpair(CHILD, PARENT, AF_UNIX, SOCK_STREAM, PF_UNSPEC) or die "ERROR: Failed to open socketpair: $!";
CHILD->autoflush(1);
PARENT->autoflush(1);

$config->{start_time} = [gettimeofday];

# check if we are forking
if(defined($config->{fork})) {
	die "ERROR: Could not find input file\n" unless defined($config->{input_file}) && -e $config->{input_file};
	die "ERROR: Input file is not bgzipped, cannot fork\n" unless $config->{input_file} =~ /\.gz$/;
	die "ERROR: Tabix index file ".$config->{input_file}.".tbi not found, cannot fork\n" unless $config->{input_file}.'.tbi';
	
	run_forks($config);
}
else {
	main($config);
}

my $elapsed = tv_interval($config->{start_time}, [gettimeofday]);

debug($config, "Took $elapsed s");

sub configure {
	
	# COMMAND LINE OPTIONS
	######################
	
	# get command-line options
	my $config = {};
	my $args = scalar @ARGV;
	
	GetOptions(
		$config,
		
		'help|h',
		'input_file|i=s',
		'tmpdir=s',
		'tmpfile=s',
		'config=s',
		
		'progress_update=i',
		'no_progress',
		
		'species=s',
		'registry|r=s',
		'host=s',
		'database|db=s',
		'user=s',
		'password=s',
		'port=i',
		
		'sql=s',
		'coord_system=s',
		
		'source=s',
		'source_desc=s',
		'population|pop=s',
		'pedigree=s',
		'panel=s',
		'gmaf=s',
		'somatic',
		
		'flank=s',
		'gp',
    'remap=s',
		'ind_prefix=s',
		'pop_prefix=s',
		'var_prefix=s',
		
		'disable_keys',
		'tables=s',
		'skip_tables=s',
		'add_tables=s',
		
		'only_existing',
    'no_merge',
		'skip_n',
		'mart_genotypes',
		
		'create_name',
		'chrom_regexp=s',
		'force_no_var',
    'ss_ids',
		
		'fork=i',
		'test=i',
		'backup',
		'move',
		
		'recover',
		'recover_point=s',
		'recover_pos=s',
		'no_recover',
		
		'cache=s',
		'fasta=s',
		
	# die if we can't parse arguments - better to get user to sort out their command line
	# than potentially do the wrong thing
	) or die "ERROR: Failed to parse command line arguments - check the documentation!\n";
	
	
	# print usage message if requested or no args supplied
	if(defined($config->{help}) || !$args) {
		&usage;
		exit(0);
	}
	
	# read config from file?
  read_config_from_file($config, $config->{config}) if defined $config->{config};
	
	# remove recover from config otherwise changes session id
	my $recover;
	if(defined($config->{recover})) {
		$recover = 1;
		delete $config->{recover};
	}
	
	# create session ID from config
	$config->{session_id} = md5_hex(join '|', (map {$_." ".$config->{$_}} sort keys %$config));
	
	debug($config, "Session ID is ", $config->{session_id});
	
	$config->{recover} = 1 if $recover;
	
	# sanity checks
	die("ERROR: Cannot run in test mode using forks\n") if defined($config->{fork}) and defined($config->{test});
	die("ERROR: Cannot manually recover using forks\n") if defined($config->{fork}) and (defined($config->{recover_point}) || defined($config->{recover_pos}));
	
	# recover?
	if(defined $config->{recover}) {
		open IN, (join '/', ($ENV{HOME}, '.import_vcf', $config->{session_id})) or die "ERROR: Could not recover session with ID ".$config->{session_id}."\n";
		my $a = <IN>;
		($config->{recover_point}, $config->{pid}) = split ' ', $a;
		close IN;
		
		die "ERROR: Cannot recover - session with ID ".$config->{session_id}." has been flagged as finished" if $config->{recover_point} eq 'FINISHED';
		
		debug($config, "Found recover point ", $config->{recover_point}, " for session ID ", $config->{session_id}, " PID ", $config->{pid});
	}
	
	# set defaults
	$config->{species}         ||= "human";
	$config->{flank}           ||= 200;
	$config->{port}            ||= 3306;
	$config->{format}            = 'vcf';
	$config->{ind_prefix}      ||= '';
	$config->{pop_prefix}      ||= '';
	$config->{coord_system}    ||= 'chromosome';
	$config->{progress_update} ||= 100;
	$config->{pid}             ||= $$;
	$config->{somatic}         ||= 0;
  
  # check remap arg
  if(defined($config->{remap})) {
    die "ERROR: remap argument incorrect - must be of format --remap [from_assembly],[to_assembly]" unless $config->{remap} =~ /^\S+?\,\S+?$/;
    
    $config->{remap} = [split(',', $config->{remap})];
  }
	
	# recovery not possible if forking
	#$config->{no_recover} = 1 if defined($config->{fork});
	
	# set default list of tables to write to
	my $tables = {
		'variation'                       => 1,
		'variation_feature'               => 1,
		'variation_synonym'               => 1,
		'flanking_sequence'               => 1,
		'allele'                          => 1,
		'population_genotype'             => 1,
		'compressed_genotype_var'         => 1,
		'individual_genotype_multiple_bp' => 0,
		'compressed_genotype_region'      => 0,
		'transcript_variation'            => 0,
		'sample'                          => 1,
		'population'                      => 1,
		'individual'                      => 1,
		'individual_population'           => 1,
		'allele_code'                     => 1,
		'genotype_code'                   => 1,
		'meta_coord'                      => 1,
		'source'                          => 1,
	};
	
	# override this with options if provided
	if(defined($config->{tables})) {
		
		# reset
		$tables->{$_} = 0 foreach keys %$tables;
		
		# set include tables
		foreach my $table(split /\,/, $config->{tables}) {
			$tables->{$table} = 1 if defined($tables->{$table});
		}
	}
	
	if(defined($config->{add_tables})) {
		
		# add tables
		foreach my $table(split /\,/, $config->{add_tables}) {
			$tables->{$table} = 1 if defined($tables->{$table});
		}
	}
	
	if(defined($config->{skip_tables})) {
		
		# set skip tables
		foreach my $table(split /\,/, $config->{skip_tables}) {
			$tables->{$table} = 0 if defined($tables->{$table});
		}
	}
	
	# force some back in
	$tables->{$_} = 1 for qw/source meta_coord/;
	
	$tables->{individual_genotype_multiple_bp} = 1 if defined($config->{mart_genotypes});
	
	# special case for sample, we also want to set the other sample tables to on
	if($tables->{sample}) {
		$tables->{$_} = 1 for qw/population individual individual_population/;
	}
	
	# force population if user wants allele or population_genotype
	if($tables->{allele} || $tables->{population_genotype}) {
		$tables->{$_} = 1 for qw/population sample/;
	}
	
	# force individual tables for individual level data
	if($tables->{compressed_genotype_region} || $tables->{individual_genotype_multiple_bp} || $tables->{compressed_genotype_var}) {
		$tables->{$_} = 1 for qw/population individual individual_population sample variation variation_synonym variation_feature flanking_sequence/;
	}
	
	if($tables->{population_genotype} || $tables->{compressed_genotype_region} || $tables->{compressed_genotype_var}) {
		$tables->{$_} = 1 for qw/genotype_code variation variation_synonym variation_feature flanking_sequence/;
	}
	
	if($tables->{allele} || $tables->{genotype_code}) {
		$tables->{$_} = 1 for qw/allele_code variation variation_synonym variation_feature flanking_sequence/;
	}
	
	# won't be writing to these tables if only_existing mode
	if(defined $config->{only_existing}) {
		$tables->{$_} = 0 for qw/source variation variation_synonym variation_feature flanking_sequence/;
	}
	
	# check that at least one has been set
	die "ERROR: no tables left included\n" unless grep {$tables->{$_}} keys %$tables;
	
	$config->{tables} = $tables;
	
	die "ERROR: tmpdir not specified\n" if !defined $config->{tmpdir} && $tables->{compressed_genotype_region};
	$config->{tmpfile} ||= 'compress.txt';
	$ImportUtils::TMP_DIR  = $config->{tmpdir};
	$ImportUtils::TMP_FILE = $config->{tmpfile};
	
	
	$config->{reg} = 'Bio::EnsEMBL::Registry';
	
	# VEP stuff
	$config->{vep}->{chunk_size}        ||= 50000;
	$config->{vep}->{cache_region_size} ||= 1000000;
	$config->{vep}->{compress}          ||= 'zcat';
	$config->{vep}->{terms}               = 'SO';
	$config->{vep}->{tr_cache}            = {};
	$config->{vep}->{rf_cache}            = {};
	$config->{vep}->{quiet}               = 1;
	$config->{vep}->{original}            = 1;
	$config->{vep}->{no_progress}         = 1;
	$config->{buffer_size}              ||= 100;
	$config->{vep}->{species}             = $config->{species};
	
	if(defined($config->{cache})) {
		die("ERROR: Could not find cache directory ".$config->{cache}) unless -e $config->{cache};		
		$config->{vep}->{dir}             = $config->{cache};
		$config->{vep}->{cache}           = 1;
		$config->{vep}->{offline}         = 1;
	}
	
	
    if(defined($config->{fasta})) {
        die "ERROR: Specified FASTA file/directory not found" unless -e $config->{fasta};
        
        eval q{ use Bio::DB::Fasta; };
        
        if($@) {
            die("ERROR: Could not load required BioPerl module\n");
        }
        
        # try to overwrite sequence method in Slice
        eval q{
            package Bio::EnsEMBL::Slice;
            
            # define a global variable so that we can pull in config hash
            our $config;
            
            {
                # don't want a redefine warning spat out, thanks
                no warnings 'redefine';
                
                # overwrite seq method to read from FASTA DB
                sub seq {
                    my $self = shift;
                    
                    # special case for in-between (insert) coordinates
                    return '' if($self->start() == $self->end() + 1);
                    
                    my $seq;
                    
                    if(defined($config->{fasta_db})) {
                        $seq = $config->{fasta_db}->seq($self->seq_region_name, $self->start => $self->end);
                        reverse_comp(\$seq) if $self->strand < 0;
                    }
                    
                    else {
                        return $self->{'seq'} if($self->{'seq'});
                      
                        if($self->adaptor()) {
                          my $seqAdaptor = $self->adaptor()->db()->get_SequenceAdaptor();
                          return ${$seqAdaptor->fetch_by_Slice_start_end_strand($self,1,undef,1)};
                        }
                    }
                    
                    # default to a string of Ns if we couldn't get sequence
                    $seq ||= 'N' x $self->length();
                    
                    return $seq;
                }
            }
            
            1;
        };
        
        if($@) {
            die("ERROR: Could not redefine sequence method\n");
        }
        
        # copy to Slice for offline sequence fetching
        $Bio::EnsEMBL::Slice::config = $config->{vep};
        
        # spoof a coordinate system
        $config->{vep}->{coord_system} = Bio::EnsEMBL::CoordSystem->new(
            -NAME => 'chromosome',
            -RANK => 1,
        );
        
        debug($config, "Checking/creating FASTA index");
        $config->{vep}->{fasta_db} = Bio::DB::Fasta->new($config->{fasta});
    }
	
	# get terminal width for progress bars
	my $width;
	
	# module may not be installed
	eval q{
		use Term::ReadKey;
	};
	
	if(!$@) {
		my ($w, $h);
		
		# module may be installed, but e.g.
		eval {
			($w, $h) = GetTerminalSize();
		};
		
		$width = $w if defined $w;
	}
	
	$width ||= 60;
	$width -= 12;
	$config->{terminal_width} = $width;
	
	return $config;
}

# reads config from a file
sub read_config_from_file {
    my $config = shift;
    my $file = shift;
    
    open CONFIG, $file or die "ERROR: Could not open config file \"$file\"\n";
    
    while(<CONFIG>) {
        next if /^\#/;
        my @split = split /\s+|\=/;
        my $key = shift @split;
        $key =~ s/^\-//g;
        
        if(defined($config->{$key}) && ref($config->{$key}) eq 'ARRAY') {
            push @{$config->{$key}}, @split;
        }
        else {
            $config->{$key} ||= $split[0];
        }
    }
    
    close CONFIG;
    
    debug($config, "Read configuration from $file") unless defined($config->{quiet});
}


# main sub-routine does most of the code execution
sub main {
	
	my $config = shift;
	
	# log start time
	my $start_time = time();
	
	
	# SET UP VARIABLES
	##################
	
	my (
		%headers,
		$prev_seq_region,
		$genotypes,
		$var_counter,
		@vf_buffer,
	);
	
	
	# CONNECT TO DBS
	################
	
	connect_to_dbs($config);
	
	# populate from SQL?
	if(defined($config->{sql}) && !defined($config->{recover})) {
		sql_populate($config, $config->{sql});
		
		if(defined($config->{test})) {
			debug($config, "(TEST) Adding schema version to meta table");
		}
		else {
			# add schema version to meta
			my $sth = $config->{dbVar}->prepare(qq{
				INSERT INTO meta (species_id, meta_key, meta_value)
				VALUES (NULL, ?, ?)
			});
			$sth->execute('schema_version', $config->{reg}->software_version);
			$sth->finish();
		}
		
		# try and do attrib
		attrib($config);
	}
	
	# get adaptors
	get_adaptors($config);
	
	# insert from files if recovering
	if(defined($config->{recover})) {
		debug($config, "Importing data from recovered session's temporary files");
		
		foreach my $table(qw(allele population_genotype compressed_genotype_region)) {
			debug($config, "Importing data into $table");
			import_tmp_file($config, $table);
		}
	}
	
	# backup
	backup($config) if defined($config->{backup}) || defined($config->{move});
	
	# TEST MODE?
	############
	
	debug($config, "Running in test mode - will read first ", $config->{test}, " lines of input") if defined($config->{test});
	
	# DB PREP
	#########
	
	# get seq_region_id hash
	$config->{seq_region_ids} = get_seq_region_ids($config);
	
	# if failed, try and copy from core DB
	if(!scalar keys %{$config->{seq_region_ids}}) {
		
		copy_seq_region_from_core($config);
		
		# now reload
		$config->{seq_region_ids} = get_seq_region_ids($config);
	}
	
	die("ERROR: seq_region not populated\n") unless scalar keys %{$config->{seq_region_ids}};
	
	# get/set source_id
	die("ERROR: no source specified\n") if !(defined $config->{source}) && !defined($config->{only_existing});
	$config->{source_id} = get_source_id($config) unless defined($config->{only_existing});
	
	# get population object
	if($config->{tables}->{population}) {
		die("ERROR: no population specified\n") unless defined $config->{population} || defined $config->{panel};
		$config->{populations} = population($config);
	}
	
	# disable keys if requested
	if(defined $config->{disable_keys}) {
		debug($config, "Disabling keys");
		
		foreach my $table(qw(allele population_genotype)) {#grep {$config->{tables}->{$_}} keys %$config->{tables}) {
			$config->{dbVar}->do(qq{ALTER TABLE $table DISABLE KEYS;});
		}
	}
	
	# GET INPUT FILE HANDLE
	#######################
	
	# if we forked, already have file handle but need to reopen original file
	# to read data from the column definition headers
	if(defined($config->{forked})) {
		my $tmp_file_handle = get_input_file_handle($config);
		
		while(<$tmp_file_handle>) {
			chomp;
			next if /^##/;
			
			my @split = split /\t/;
			
			# column definition line
			if(/^#/) {
				%headers = %{parse_header($config, \@split)};
				last;
			}
		}
	}
	
	else {
		$config->{in_file_handle} = get_input_file_handle($config);
	}
	
	
	# PEDIGREE FILE
	###############
	
	$config->{pedigree} = pedigree($config) if defined($config->{pedigree});
	
	# MAIN FILE LOOP
	################
	
	my $in_file_handle = $config->{in_file_handle};
	
	my $last_skipped = 0;
	
	# read the file
	while(<$in_file_handle>) {
		chomp;
		
		# header lines
		next if /^##/;
		
		my @split = split /\s+/;
		my $data = {};
		$data->{line} = $_;
		
		# column definition line
		if(/^#/) {
			%headers = %{parse_header($config, \@split)};
			
			# leave a file telling the master process to fork the others
			if(defined($config->{forked})) {
				open TMP, '> '.$ENV{HOME}.'/.import_vcf/'.$config->{pid};
				print TMP '1';
				close TMP;
			}
		}
		
		# data
		else {
			
			$config->{prev_time} ||= [gettimeofday];
			
			# recover?
			if(defined($config->{recover_point}) || defined($config->{recover_pos})) {
				$config->{skipped}->{already_processed}++;
				
				if(defined($config->{recover_point}) && md5_hex($data->{line}) eq $config->{recover_point}) {
					delete $config->{recover_point};
					$config->{recover_check} = 1;
					debug(
						$config,
						"Found recovery point, skipped ",
						$config->{skipped}->{already_processed},
						" already processed variants"
					);
				}
				
				if(defined($config->{recover_pos})) {
					my $rpos = $config->{recover_pos};
					
					if($rpos =~ /\:/) {
						my ($chr, $pos) = split /\:/, $rpos;
						if($data->{line} =~ /^$chr\s+$pos\s+/) {
							delete $config->{recover_pos};
							$config->{recover_check} = 1;
							debug(
								$config,
								"Found recovery point, skipped ",
								$config->{skipped}->{already_processed},
								" already processed variants"
							);
						}
					}
					elsif($data->{line} =~ /^\w+?\s+$rpos\s+/) {
						delete $config->{recover_pos};
						$config->{recover_check} = 1;
						debug(
							$config,
							"Found recovery point, skipped ",
							$config->{skipped}->{already_processed},
							" already processed variants"
						);
					}
				}
				
				next;
			}
			
			# check we're not skipping loads in a row
			if($last_skipped > 100 && $last_skipped =~ /(5|0)00$/) {
				debug($config, "WARNING: Skipped last $last_skipped variants, are you sure this is running OK? Maybe --gp is enabled when it shouldn't be, or vice versa?");
			}
			
			# parse into a hash
			$data->{$_} = $split[$headers{$_}] for keys %headers;
			
			# skip non-variant lines
			if($data->{ALT} eq '.') {
				$config->{skipped}->{non_variant}++;
				next;
			}
			
			# parse info column
			my %info;	
			foreach my $chunk(split /\;/, $data->{INFO}) {
				my ($key, $val) = split /\=/, $chunk;
				$info{$key} = $val;
			}
			
			$data->{info} = \%info;
			
			# skip unwanted chromosomes
			#next if defined($config->{chrom_regexp}) && $data->{'#CHROM'} !~ m/$chrom_regexp/;
			
			# use VEP's parse_line to get a skeleton VF
			($data->{tmp_vf}) = @{parse_line($config, $data->{line})};
			
			if(!defined($data->{tmp_vf})) {
				$config->{skipped}->{could_not_parse}++;
				$last_skipped++;
				next;
			}
      
      # remap?
      if(defined($config->{remap})) {
        my $success = remap($config, $data->{tmp_vf});
        if(!$success) {
          $config->{skipped}->{unable_to_remap}++;
          $last_skipped++;
          next;
        }
      }
			
			if(!defined($config->{seq_region_ids}->{$data->{tmp_vf}->{chr}})) {
				$config->{skipped}->{missing_seq_region}++;
				$last_skipped++;
				next;
			}
			
			# copy seq region ID
			$data->{tmp_vf}->{seq_region_id} = $config->{seq_region_ids}->{$data->{tmp_vf}->{chr}};
			
			# could be a structural variation feature
			next unless $data->{tmp_vf}->isa('Bio::EnsEMBL::Variation::VariationFeature');
      
      # ssIDs as IDs?
      if(defined($config->{ss_ids})) {
        my ($ss_id) = grep {$_ =~ /^\d+$/} split(/\;/, $data->{ID});
        
        if(defined($ss_id)) {
          $data->{SS_ID} = $ss_id;
        }
      }
			
			# sometimes ID has many IDs separated by ";", take the lowest rs number, otherwise the first
			if($data->{ID} =~ /\;/) {
				if($data->{ID} =~ /rs/) {
					$data->{ID} = (sort {(split("rs", $a))[-1] <=> (split("rs", $b))[-1]} split /\;/, $data->{ID})[0];
				}
        elsif(defined($config->{ss_ids}) && defined($data->{SS_ID})) {
          $data->{ID} = 'ss'.$data->{SS_ID};
        }
				else {
					$data->{ID} = (split /\;/, $data->{ID})[0];
				}
			}
      elsif(defined($config->{ss_ids}) && defined($data->{SS_ID})) {
        $data->{ID} = 'ss'.$data->{SS_ID};
      }
			
			# make a var name if none exists
			if(!defined($data->{ID}) || $data->{ID} eq '.' || defined($config->{create_name})) {
				$data->{ID} =
					($config->{var_prefix} ? $config->{var_prefix} : 'tmp').
					'_'.$data->{'#CHROM'}.'_'.$data->{POS}.'_'.$data->{REF}.'_'.$data->{ALT};
				$data->{made_up_name} = 1;
			}
			
			$data->{tmp_vf}->{variation_name} = $data->{ID};
			
			# parse genotypes
			$data->{genotypes} = get_genotypes($config, $data, \@split);
			
			# get variation object
			$data->{variation} = variation($config, $data);
			
			# transcript variation (get cons)
			get_all_consequences($config->{vep}, [$data->{tmp_vf}]) if $config->{tables}->{transcript_variation};
			
			# get variation_feature object
			$data->{vf} = variation_feature($config, $data);
			
			# attach variation to genotypes
			$_->{variation} = $data->{variation} for @{$data->{genotypes}};
			
			# skip variation if no dbID
			if(!defined($data->{variation}->{dbID}) && !defined($config->{test})) {
				$config->{skipped}->{var_not_present}++;
				$last_skipped++;
				next;
			}
			
			# transcript variation (write to DB)
			transcript_variation($config, [$data->{tmp_vf}]) if $config->{tables}->{transcript_variation} && defined($data->{tmp_vf}->dbID);
			
			#if($config->{tables}->{transcript_variation}) {
			#	push @vf_buffer, $data->{vf};
			#	if(scalar @vf_buffer == $config->{buffer_size}) {
			#		transcript_variation($config, \@vf_buffer);
			#		@vf_buffer = ();
			#	}
			#}
			
			# alleles
			allele($config, $data) if $config->{tables}->{allele};
			
			# population genotypes
			population_genotype($config, $data) if $config->{tables}->{population_genotype};
			
			# individual genotypes
			individual_genotype($config, $data) if $config->{tables}->{compressed_genotype_var} || defined($config->{mart_genotypes});
			
			# GENOTYPES BY REGION
			#####################
			
			# multi bp
			#if($config->{tables}->{individual_genotype_multiple_bp} && $force_multi && @{$data->{genotypes}}) {
			#	&multi_bp_genotype($dbVar, $data);
			#}
			
			# compressed by region
			if($config->{tables}->{compressed_genotype_region} && @{$data->{genotypes}} && !defined($config->{test})) {
				my $vf = $data->{vf};
				
				$vf->{seq_region_id} = $vf->slice->get_seq_region_id if !defined($vf->{seq_region_id});
				
				if(defined($vf->{seq_region_id}) && defined($vf->{start})) {
				
					foreach my $gt(@{$data->{genotypes}}) {
						my $sample_id = $gt->individual->dbID;
						
						next if $gt->genotype_string =~ /\./;
						
						# add to compress hash for writing later
						if (!defined $genotypes->{$sample_id}->{region_start}){
							$genotypes->{$sample_id}->{region_start} = $vf->{start};
							$genotypes->{$sample_id}->{region_end} = $vf->{end};
						}
						
						# write previous data?
						#compare with the beginning of the region if it is within the DISTANCE of compression
						if (
							defined($genotypes->{$sample_id}->{genotypes}) &&
							(
								(abs($genotypes->{$sample_id}->{region_start} - $vf->{start}) > DISTANCE()) ||
								(abs($vf->{start} - $genotypes->{$sample_id}->{region_end}) > MAX_SHORT) ||
								(defined($prev_seq_region) && $vf->{seq_region_id} != $prev_seq_region) ||
								($vf->{start} - $genotypes->{$sample_id}->{region_end} - 1 < 0)
							)
						) {
							#snp outside the region, print the region for the sample we have already visited and start a new one
							print_file($config,$genotypes, $prev_seq_region, $sample_id);
							delete $genotypes->{$sample_id}; #and remove the printed entry
							$genotypes->{$sample_id}->{region_start} = $vf->{start};
						}
						
						if ($vf->{start} != $genotypes->{$sample_id}->{region_start}){
							#compress information
							my $blob = pack ("w",$vf->{start} - $genotypes->{$sample_id}->{region_end} - 1);
							$genotypes->{$sample_id}->{genotypes} .=
								escape($blob).
								escape(pack("w", $data->{variation}->dbID || 0)).
								escape(pack("w", $config->{individualgenotype_adaptor}->_genotype_code($gt->genotype, $gt->phased)));
						}
						else{
							#first genotype starts in the region_start, not necessary the number
							$genotypes->{$sample_id}->{genotypes} =
								escape(pack("w", $data->{variation}->dbID || 0)).
								escape(pack("w", $config->{individualgenotype_adaptor}->_genotype_code($gt->genotype, $gt->phased)));
						}
						
						$genotypes->{$sample_id}->{region_end} = $vf->{start};
					}
				}
			}
			
			$prev_seq_region = $data->{vf}->{seq_region_id};
			$last_skipped = 0;
			
			$var_counter++;
{'fraction_non_alphanumeric': 0.18420650095602295, 'fraction_numerical': 0.00497131931166348, 'mean_word_length': 2.8006685560642395, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 853, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4711772', 'n_tokens_mistral': 20872, 'n_tokens_neox': 17571, 'n_words': 3853}