reStraining

#!/usr/bin/env perl
use warnings;
use strict;
use Getopt::Long;
use FindBin qw($Bin);
use lib "$Bin/../lib";
use Cwd;

## This program is Copyright (C) 2018-23, Felix Krueger (fkrueger@altoslabs.com)

## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.

## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.

## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.


### This script filters the latest VCF file for various SNPs versus the GRCm39
### mouse genome build and writes high confidence SNPs into a folder called 'SNPs_Sanger';
### it then creates a new genome called MGP_strains_N-masked.

### It has been tested with the latest version of the SNP file 'mgp_REL2021_snps.vcf.gz'
### at the time of writing (02 April 2023). I can't guarantee that it will work with any other VCF file

# Updating 02 April 2023 to change over to GRCm39 and v8 annotations


## Reading in a BAM or SAM file
my $pipeline_version = '0.4.0';
my $parent_dir = getcwd();
my %strains;
my ($vcf_file,$genome_folder,$nmasking,$genome_build,$v7) = process_commandline ();

my %snps; # storing all filtered SNPs

warn "\nSummarising reStraining parameters - Genome Preparation\n";

warn "="x50,"\n";

warn "Processing SNPs from VCF file:\t\t$vcf_file\n";
warn "Reference genome:\t\t\t$genome_folder\n";
warn "\n";

### Strain information
print "Generating N-masking output for all of the following strains:\n";
print "="x37,"\n";  
foreach my $ind(sort {$a<=>$b} keys %strains){
      print "$strains{$ind}->{strain}  //  ";
}
print "\n";
print "="x37,"\n\n";

if ($v7){
    warn "Using v7 file $vcf_file\n";
}

### Dealing with chromosomes - At the moment we are trying to limit the analysis on chromosome 1 [mainly for time and memory reasons]
#my @chroms = (1); # this is to just use chromosome 1 
my @chroms = detect_chroms(); # this alternate method will detect all chromosomes from the VCF header and use those instead.

print "Using the following chromosome(s):\t";
# print "Using the following chromosomes (detected from VCF file >>$vcf_file<<):\n";
print join ("\t",@chroms),"\n\n";


### Storing the entire genome sequence
my %chromosomes; # genomic sequence
read_genome_into_memory($parent_dir);

### Determining and Filtering homozygous high-confidence SNPs for the strain in question
filter_relevant_SNP_calls_from_VCF(); # the last number is the strain identity, here the first strain
warn "Finished filtering and writing out SNPs\n\n";

### Create modified genome
my $new_n_total = 0;
my $new_snp_total = 0;
my $already_total = 0;
my $low_confidence = 0;


# Writing a genome generation report file
my $report = "reStrainingOrder_genome_preparation_report.txt";
open (REPORT,'>',$report) or die "Failed to write to file $report: $!\n";

for my $chr (@chroms) {
    create_modified_chromosome($chr);
}

if ($nmasking){
    warn "\n\nSummary\n$new_n_total Ns were newly introduced into the N-masked genome for all current MGP strains in total\n";
    print REPORT "\nSummary\n$new_n_total Ns were newly introduced into the N-masked genome for all current MGP strains in total\n";
}
warn "\n";
close REPORT;

warn "All done. Genome(s) are now ready to be indexed with your favourite aligner!\nFYI, aligners known to work with reStrainingOrder are Bowtie2, Tophat, STAR, HISAT2, HiCUP and Bismark (STAR and HISAT2 require disabling soft-clipping, please check the manual for details)\n\n";      



#############################################################
### SUBROUTINES
#############################################################


sub create_modified_chromosome {
    
    my ($chr) = @_; 
    warn "Processing chromosome $chr\n";
    
    unless ($chromosomes{$chr}){
	    warn "\nThe chromosome name given in the VCF file was '$chr' and was not found in the reference genome.\nA rather common mistake might be that the VCF file was downloaded from Ensembl (who use chromosome names such as 1, 2, X, MT)\nbut the genome from UCSC (who use chromosome names such as chr1, chr2, chrX, chrM)\n";
        warn "The chromosome names in the reference genome folder were:\n";
    	foreach my $c (sort keys %chromosomes){
    	    warn "$c\n";
    	}
	   die "[FATAL ERROR] Please ensure that the same version of the genome is used for both VCF annotations and reference genome (FastA files). Exiting...\n\n";
    }

    my $n_sequence = $chromosomes{$chr};

    warn "Getting SNPs for chromosome $chr...\n";
    my @snps = @{read_snps($chr)};
    warn "done.\n";
    
    unless (@snps){
	   @snps = ();
	   warn "Clearing SNP array...\n"
    }
    
    my $count = 0;

    my $lastPos = 0;

    my $already = 0;
    my $warn = 0;
    my $new_n = 0;
    my $new_snp = 0;

    foreach my $snp (@snps) {
    	# Apply the SNP
    	++$count;
    	# warn "$snp->[0]\t$snp->[1]/$snp->[2]\n";
    	if ($snp->[0] == $lastPos) {
    	    # Duplicate SNP
    	    next;
    	}
          
    	$lastPos = $snp->[0];
    	
    	# Check if the reference base is the same as the SNP base
    	#if (substr ($sequence,$snp->[0]-1,1) eq $snp->[2]) {
    	    # warn "Skipping $snp->[0] $snp->[1]/$snp->[2] since the ref and SNP base are the same\n";
    	 #   ++$already;
    	 #   next;
    	#}
          
    	# Check the reference base is correct
    	if (substr ($n_sequence,$snp->[0]-1,1) ne $snp->[1]) {
    	    warn "Skipping $snp->[0] $snp->[1]/$snp->[2] since the reference base didn't match\n";
    	    $warn++;
    	    next;
    	}
          
    	### Ref/Alt bases are matching, so we can proceed to changing the ref base for the SNP base or Ns (N-masking)
    	
    	### N-masking 
    	if ($nmasking){ # default
    	    my $return =  substr($n_sequence,($snp->[0])-1,1,'N');  # Replacing the base with 'N'
    	    unless ($return){
    		  warn "Replacing failed...\n";
    	    }
    	    ++$new_n;
    	}
	
    }
  
    $new_n_total += $new_n;
    $new_snp_total += $new_snp;
     
    $already_total += $already;
    if ($nmasking){ 
	   write_SNP_chromosome($chr,$n_sequence,1);
    }
    
    warn "$count SNPs total for chromosome $chr\n";
    if ($nmasking){ # default
	   warn "$new_n positions on chromosome $chr were changed to 'N'\n";
	   print REPORT "$new_n positions on chromosome $chr were changed to 'N'\n";
    }
    warn "\n";

}


sub write_SNP_chromosome {

    my ($chr,$sequence,$nm) = @_; # $nm will discriminate between N-masking and full sequence output
    if ($nm){
	   warn "Writing modified chromosome (N-masking)\n";
    }
    else{
	   die "Not N-masking is not an option!\n";
    }
    
    my $type;
    my $outfile;
    
    if ($nm){
	   $type = 'N-masked';
	   $outfile = "chr${chr}.N-masked.fa";
    }
    
    # warn "Starting sequence is ".length($sequence)." bp\n";
    if ($nm){
	    warn "Writing N-masked output to: ${parent_dir}/MGP_strains_${type}/$outfile\n";
	    unless (-d "${parent_dir}/MGP_strains_${type}/"){ # creating the output directory if required
	       mkdir "${parent_dir}/MGP_strains_${type}/";
	    }      
    	open (OUT,'>',"${parent_dir}/MGP_strains_${type}/${outfile}") or die "Failed to write to file ${parent_dir}/MGP_strains_${type}/${outfile}: $!\n\n";
    	print OUT ">$chr\n";
    }
    
    my $pos = 0;

    # Writing out chromosome files with 100 characters per line
    while ($pos < length($sequence)-100) {
	   print OUT substr($sequence,$pos,100),"\n";
	   $pos += 100;
    }
    print OUT substr($sequence,$pos),"\n"; # rest
    close OUT or die $!;
   
}


sub read_snps {
    
    my ($chr) = @_;
    my @snps = ();
    my $file = "${parent_dir}/SNPs_directory/chr$chr.txt";

    ### If the SNP folder doesn't exist we can be certain that something is going wrong
    unless (-d "${parent_dir}/SNPs_directory"){
	   die "Folder >>${parent_dir}/SNPs_directory<< did not exist, I'm afraid something seems to have gone wrong...\n\n";
    }

    ### not sure but I think for some chromosomes there might not be any SNP files, e.g. chr MT or chrY. In this case the sequence is written out again unmodified
    unless (-e $file) {
	   warn "Couldn't find SNP file for chromosome '$chr' '$file' didn't exist. Skipping...\n";
	   return \@snps;
    }
    warn "Reading SNPs from file $file\n";
    
    open (IN,$file) or die $!;
    $_ = <IN>; # chromosome name header line
    while (<IN>) {
	    $_ =~ s/\r//; # Windows line endings...
	    chomp;
	    # warn "$_\n"; sleep(1);
	    next unless ($_);
	    my (undef,$pos,$ref_allele,$snp_allele) = split(/\t/);
	    # warn  "$pos , $ref_allele , $snp_allele\n"; sleep(1);
    	push @snps,[$pos,$ref_allele,$snp_allele];
    }

    # sorting snps
    @snps = sort {$a->[0] <=> $b->[0]} @snps;

    return \@snps;
    
    close IN or warn "Failed to close filehandle IN for file $file: $!\n\n";

}

###

sub filter_relevant_SNP_calls_from_VCF{
    
    warn "Now reading from VCF file $vcf_file\n";
	# warn "STRAIN is: $strain\n\n"; sleep(10);
    if ($vcf_file =~ /gz$/){
	   open (IN,"gunzip -c $vcf_file |") or die "Failed to open file '$vcf_file': $!\n";
    } 
    else{
	   open (IN, $vcf_file) or die "Failed to read Input VCF file '$vcf_file': $!\n";
    }
    
    my $matrix_file;
    
    if ($v7){
        $matrix_file = "MGPv7_SNP_matrix_chr1.txt.gz";
    } 
    else{
        $matrix_file = "MGPv8_SNP_matrix_chr1.txt.gz";
    }

    ### Printing the SNP matrix to this file, which will be important for the reStrainingOrder process
    warn "Writing a SNP matrix for chromosome 1 for all strains to $matrix_file\n\n";
    open (THEMATRIX,"| gzip -c - > $matrix_file") or die "Failed to write to file '$matrix_file': $!\n";

    print THEMATRIX join ("\t","Chromosome","Position","REF","ALT");
     foreach my $index (sort {$a <=> $b} keys %strains){
        # warn "$index\t$strains{$index}->{strain}\n";
        print THEMATRIX "\t$strains{$index}->{strain}";
    }
    print THEMATRIX "\n";
    

    my %all_SNPs; # storing filtered SNPs
    my $count = 0;
    my $matrix_count = 0;
    my $positions_discarded = 0;
    
    my $indel_pos = 0;
    my $positions_skipped = 0;
    my %fhs;
    
    my $format_index;    # required to get extract entries from FORMAT field
	my $info_index;      # required to look at INFO field, e.g. for INDELs
	my $gt_index;        # required to get GENOTYPE
	my $fi_index;        # required to get FILTER value
	
    my $dir = "SNPs_directory";
    unless (-d $dir){
		warn "Folder '$dir' doesn't exist. Creating it for you...\n\n";
		mkdir $dir or die "Failed to created directory $dir\n: $!\n\n";
	}
    
    # Opening filehandles for the SNP files
    for my $chr (@chroms) {
    	my $filename = "SNPs_directory/chr".$chr.'.txt';
    	open (my $fh,'>',$filename) or die "Couldn't open filehandle $!\n";
    	$fhs{$chr} = $fh;
        print {$fhs{$chr}} ">$chr\n";
    }  

    my $last_chr;

    while (<IN>){
    	$_ =~ s/(\r|\n)//g; # removing end of line characters
        
        next if ($_ =~ /^\#\#/); # filters out header information lines
    	# warn "$_\n"; sleep(1);
        if ($_ =~ /^\#CHROM/){ # Table Header
    	
        	# my ($name) = (split /\t/)[$strain_index];
			# warn "Analysing SNP fields for name >$name<\n";
            # sleep(5);
			my @format_fields = split /\t/;
			my $field_index = 0;
			foreach my $field (@format_fields){
				# warn "$field_index\t$field\n";#
				if ($field eq "FORMAT"){#
					$format_index = $field_index;
				}
				if ($field eq "INFO"){#
					$info_index = $field_index;
				}
				$field_index++;
			}
			
			if (defined $format_index){
				# warn "Using FORMAT field index: $format_index\n";
                # sleep(1);
			}
			else{
				die "Failed to extract index of field 'FORMAT'. Hmmm...";
			}

			if (defined $info_index){
				# warn "Using INFO field index: $info_index\n";
                # sleep(1);
			}
			else{
				die "Failed to extract index of field 'INFO'. Hmmm...";
			}

            #my ($name) = (split /\t/)[$strain_index];
    	    #warn "Analysing SNP fields for name >$name<\n";
    	    next;	
    	}
    	$count++;
    	if ($count%1000000 ==0){
    	    warn "processed $count lines\n";
    	}
	    # last if ($count == 1000);
		# warn "$_\n"; sleep(1);
        # foreach my $index(0..$#strains){
        # next if ($index <= 8); 
        # $strains{$strains[$index]} = $index; 
        # warn "$index\t$strains[$index]\n";
        #   }
        # }
    
	    # This is the format of the VCF (v5) file
        #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  129P2_OlaHsd    129S1_SvImJ     129S5SvEvBrd    AKR_J   A_J     BALB_cJ BTBR_T+_Itpr3tf_J       BUB_BnJ C3H_HeH C3H_HeJ C57BL_10J       C57BL_6NJ       C57BR_cdJ      C57L_J  C58_J   CAST_EiJ        CBA_J   DBA_1J  DBA_2J  FVB_NJ  I_LnJ   KK_HiJ  LEWES_EiJ       LP_J    MOLF_EiJ        NOD_ShiLtJ      NZB_B1NJ        NZO_HlLtJ       NZW_LacJ        PWK_PhJ RF_J    SEA_GnJ SPRET_EiJ     ST_bJ    WSB_EiJ ZALENDE_EiJ
        
        # This is the format of the v7 VCF file:
        #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	129P2_OlaHsd	129S1_SvImJ	129S5SvEvBrd	A_J	AKR_J	B10.RIII	BALB_cByJ	BALB_cJ	BTBR_T+_Itpr3tf_J	BUB_BnJ	C3H_HeH	C3H_HeJ	C57BL_10J	C57BL_10SnJ	C57BL_6NJ	C57BR_cdJ	C57L_J	C58_J	CAST_EiJ	CBA_J	CE_J	CZECHII_EiJ	DBA_1J	DBA_2J	FVB_NJ	I_LnJ	JF1_MsJ	KK_HiJ	LEWES_EiJ	LG_J	LP_J	MA_MyJ	MOLF_EiJ	NOD_ShiLtJ	NON_LtJ	NZB_B1NJ	NZO_HlLtJ	NZW_LacJ	PL_J	PWK_PhJ	QSi3	QSi5	RF_J	RIIIS_J	SEA_GnJ	SJL_J	SM_J	SPRET_EiJ-ST_bJ	SWR_J	WSB_EiJ	ZALENDE_EiJ
        
        # 10 04 2023: This is the format of the v8 VCF file ("mgp_REL2021_snps.vcf.gz"):
        ### #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	129P2_OlaHsd	129S1_SvImJ	129S5SvEvBrd	A_J	AKR_J	B10.RIII	BALB_cByJ	BALB_cJ	BTBR_T+_Itpr3tf_J	BUB_BnJ	C3H_HeH
        ###	C3H_HeJ	C57BL_10J	C57BL_10SnJ	C57BL_6NJ	C57BR_cdJ	C57L_J	C58_J	CAST_EiJ	CBA_J	CE_J	CZECHII_EiJ	DBA_1J	DBA_2J	FVB_NJ	I_LnJ	JF1_MsJ	KK_HiJ	LEWES_EiJ	LG_J
        ### LP_J	MAMy_J	MOLF_EiJ	NOD_ShiLtJ	NON_LtJ	NZB_B1NJ	NZO_HlLtJ	NZW_LacJ	PL_J	PWK_PhJ	QSi3	QSi5	RF_J	RIIIS_J	SEA_GnJ	SJL_J	SM_J	SPRET_EiJ	ST_bJ	SWR_J	WSB_EiJ	ZALENDE_EiJ

        my ($chr,$pos,undef,$ref,$alt,undef,undef,$info,$format,@strains) = (split /\t/);

		#  warn "genotype: $gt\nfilter:   $fi\n";
		#  warn "$fi\n"; sleep(1);
        unless (defined $last_chr){
            $last_chr = $chr; # only once
            warn "Processing first chromosome (chr $chr)\n";
        }
        
        unless ($chr eq $last_chr){
            warn "Reached a new chromosome (chr $chr)\n";
            $last_chr = $chr; 
        }
        # warn "$chr , $pos , $ref , $alt\n"; sleep(1);
        unless ($chr eq '1'){
            warn "Reached a new chromosome ($chr), but we are currently limiting filtering to chromosome 1\n";
            last;
        }

		# skipping if the reference base is not well defined in Black6
		if ($ref =~ /[^ATCG]/){ # reference base contained any non A, C, T, G characters or more than one base
			# warn "ref was: $ref; skipping\n";
            ++$positions_skipped;
			next;
		}
		
        # The v7 file contains both SNPs as wella s INDELS. We are only interested in SNPs here, so removing all INDELs
        if ($v7){
            if ($info =~ /INDEL/){
                # warn "This variant is an INDEL:\n$info\nRemoving...\n"; sleep(1);
                ++$indel_pos;
                next;
            }
        }

		# skipping if the SNP is not well defined in the strain of interest
		# if ($alt =~ /[^ATCG,]/){ # Alt base contained any non A, C, T, G characters or commas which separate several different variants
        # SNPs may also have more than 1 separate ALT bases. To get the maximum out of the SNP analysis, one 
        # should probably include these positions. For reStrainingOrder it should be enough to only look at 
        # positions that have one REF and one ALT base.        
    
		if ($alt =~ /[^ATCG]/){ # Alt base contained any non A, C, T, G characters or commas which separate several different variants
            # warn "SNP was: $alt; skipping\n";
			++$positions_skipped;
			# next;
            # warn "ALT base not clearly defined:\t$_\n"; sleep(1);
            next;
		}
       
        # Rather than looping through all of the strains to find positions that are unique to one strain or another
        # (the first method), we are generating a matrix of all possible SNPs for currently Chr 1 for all strains
        my @straincalls;
        
        # The first 8 fields are irrelevant: #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  
		foreach my $index (0..$#strains){
        	# print "INDEX: $index\tSTRAIN: $strains{$index}->{strain}\t\t$strains[$index]\n"; sleep(1);
            my ($gt,$gq,$dp,$mq0f,$gp,$pl,$an,$mq,$dv,$dp4,$sp,$sgb,$pv4,$fi,$adf,$adr,$ad);
        
            # for v5: $strain is in the form: GT:GQ:DP:MQ0F:GP:PL:AN:MQ:DV:DP4:SP:SGB:PV4:FI
	        # for v7: $strain is in the form: GT:PL:DP:ADF:ADR:AD:GQ:FI
            # for v8: $strain is in the form: GT:PL:DP:AD:GQ:FI   # Now te default (as of 10 April 2023)

            if ($v7){ # the v7 version, currently experimental
 		        ($gt,$pl,$dp,$adf,$adr,$ad,$gq,$fi) = split/:/,$strains[$index];
    	    }
            else{ # mgp v5, default
                ($gt,$pl,$dp,$ad,$gq,$fi) = split/:/,$strains[$index]; 
            }

            #warn "genotype: $gt\nfilter:   $fi\n";
        	#warn "$fi\n"; sleep(1);
        	# $gt is the Genotype:
        	# print "CHR: $chr\tPOS: $pos\tREF: $ref\tALT: $alt\tGT: $gt\tFILTER: $fi\t$strains{$index}\t$strains[$index]\n"; # sleep(1);
        	# '.'   = no genotype call was made
        	# '0/0' = genotype is the same as the reference genome
        	# '1/1' = homozygous alternative allele; can also be '2/2', '3/3', etc. if more than one alternative allele is present
        	# '0/1' = heterozygous genotype; can also be '1/2', '0/2', etc.
        	
        	# $fi is FILTER, 1 for high confidence SNP, or 0 for low confidence
        	
        	### We are only looking for 1/1, (2/2 or 3/3 calls not in the first instance) and filter for high confidence as well
        	
            my $pos_fails = 0;
            my $strain_call = 0; # this is a boolean flag for SNP present yes/no
            
            # if (!defined $fi){
            #     warn "$1 was not extracted properly, please fix!\n";
            # }
	
            # Looking at the Filter tag first
            if ($fi eq 1){
            	# warn "FI was high confidence ($fi)\n";
                # fine
            }
            else{
                # warn "FI was low confidence ($fi)\n";
                $pos_fails = 1;
            }
                
        	# Filtering for genotype. If the genotype is the same as the reference we will move
            # on to the next position and/or strain
        	if ($gt eq '0/0'){
        	    # warn "same as reference, next...\n";
        	}
			
			elsif ($gt eq '1/1'){
        	    unless($pos_fails){
                    $strain_call = 1;  # This is the only time at which we are setting the $strain_call to TRUE
                }
        	    #if ($alt =~ /[^ATCG]/){ 	 
            	#	warn "homozygous alternative allele: >$alt<\n";
            	#	my ($new_alt) = (split (/,/,$alt))[0];
            	#	 warn "New ALT is $new_alt\n";
            	#	 warn "genotype: $gt\nfilter:   $fi\n\n"; sleep(1);
            	#	$alt = $new_alt;  
        	    #}
        	}
        	elsif ($gt eq '2/2'){
                warn "$gt\n";
        	    die "Should never happen\n";
                
                # ++$homozygous;
        	    if ($alt =~ /[^ATCG]/){ 	 
            		# warn "homozygous alternative allele: >$alt<\n";
            		my ($new_alt) = (split (/,/,$alt))[1];
            		# warn "New ALT is $new_alt\n\n"; sleep(1);
            		$alt = $new_alt;    
        	    }
        	}
        	elsif ($gt eq '3/3'){
                warn "$gt\n";
                die "Should never happen\n";
        	    # ++$homozygous;
        	    if ($alt =~ /[^ATCG]/){
        		    # warn "homozygous alternative allele: >$alt<\n";
        		    my ($new_alt) = (split (/,/,$alt))[2];                                                                                                        
         		    # warn "New ALT is $new_alt\n\n"; sleep(1);                                                                             
        		    $alt = $new_alt;   
        	    }
        	}
        	else{
        	    # this could be positions without genotype call ./. or any form of heterozygous call. 
			}
    
            # Adding the boolean strain call to a temporary array
            push @straincalls, $strain_call;
            # warn "$SNP\n";sleep(1);
            # if (defined $unique_snp){
                #warn "SNP was already defined by a previous strain. Failing this position out...\n";
            #    $pos_fails++;
		    #		last;
            #   }
            # else{
          
				# Output example
				# Variation ID    Chromosome name Position on Chromosome (bp)     Strand  Allele
				# rs2020560       10      98212004        1       A/T
				# my $SNP = join ("\t",$strains{$index}->{strain},$chr,$pos,$ref,$alt,$strain_call);
                # warn "CALL: $SNP\n~~~~~~~~~~~~~~~~~~\n"; sleep(1);
            #    ++$unique_single_strain;       
				#$unique_strain = $strains{$index}->{strain};
               # $unique_snp = $SNP;
              #  $unique_index = $index;
                #warn  "SNP was not present yet, storing it (for strain $strains{$index}->{strain})\n";
            #}

            # if (exists $snps{$location} ){
            #   warn "SNP $snps{$location} position was present already\n";
            # }
            # else{
            #   $snps{$location} = $SNP;
            #   warn  "High confidence SNP for strain $strains{$index}\n";
            # }
        	# print {$fhs{$chr}} join ("\t",$count,$chr,$pos,'1',join ("\/",$ref,$alt),$strain),"\n";
    	}

        # Ensure that at least one strain had a high confidence SNP call
        my $should_stay = 0;
        foreach my $call(@straincalls){
            # warn "$call"; 
            if ($call eq 1){
                ++$should_stay;
                last;
            }
        }
        # sleep(1);
        # Now that we have processed all strains we can print this position to THE MATRIX
        if ($should_stay){  
            my $snp = join ("\t",$chr,$pos,$ref,$alt,@straincalls);

            # Only printing out positions for chromosome 1 to THE MATRIX
            # Will print for Ensembl and UCSC version (changed 22 01 2021)
	        if ($chr eq 1 or $chr eq "chr1"){
                print THEMATRIX "$snp\n";	
                ++$matrix_count;
            }

            my $location = join (":",$chr,$pos);      
            # $snps{$location}->{snp} = $snp;
            #print {$fhs{$chr}} "$snps{$location}->{snp}\n"; # also writing into the SNP folder
       
            print {$fhs{$chr}} "$snp\n"; # also writing into the SNP folder
           
        }
        else{
            ++$positions_discarded;
        }
        # $snps{$location}->{strain} = $unique_strain;
        # $strains{$unique_index}->{count}++;
        # warn "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n";
    }

    # Writing a report file
    my $report = "reStrainingOrder_SNP_filtering_report.txt";
    open (REPORT,'>',$report) or die "Failed to write to file $report: $!\n";
    
    warn "\nSNP position summary for all MGP strains (based on mouse genome build $genome_build)\n";
    warn "="x77,"\n\n";
    warn "Positions read in total:\t$count\n";
    if ($v7){
		warn "$indel_pos\tPositions were INDELs (and hence skipped)\n";
	}
    warn "Positions skipped because the REF/ALT bases were not well defined:\t$positions_skipped\n";
    warn "Positions discarded as no strain had a high confidence call:\t$positions_discarded\n\n";
    warn "Positions printed to THE MATRIX in total:\t$matrix_count\n";
    # warn "$homozygous\tSNP were homozygous. Of these:\n";
    	
    print REPORT "SNP position summary for all MGP strains (based on mouse genome build $genome_build)\n";
    print REPORT "="x75,"\n\n";
    print REPORT "Positions read in total:\t$count\n";
    if ($v7){
		print REPORT "$indel_pos\tPositions were INDELs (and hence skipped)\n";
	} 
    print REPORT "Positions skipped because the REF/ALT bases were not well defined:\t$positions_skipped\n";
    print REPORT "Positions discarded as no strain had a high confidence call:\t$positions_discarded\n\n";
    print REPORT "Positions printed to THE MATRIX in total:\t$matrix_count\n";
    
    close REPORT or warn "Failed to close filehandle REPORT\n";

    #foreach my $location (keys %snps){
    #    my ($chr) = (split /\t/,$snps{$location}->{snp})[1];
        # wparn "CHR: $chr\n"; 
	 #   print ALLSNP "$snps{$location}->{snp}\n";
     #   print {$fhs{$chr}} "$snps{$location}->{snp}\n"; # also writing into the SNP folder
    #}
   
}


sub detect_chroms{
    
    my %chrom; # detecting the chromosomes from the VCF file
    my @chrom;
    
    if ($vcf_file =~ /gz$/){
	open (DETECT,"gunzip -c $vcf_file |") or die "Failed to open file '$vcf_file': $!\n";
    }
    else{
	open (DETECT, $vcf_file) or die "Failed to read Input VCF file '$vcf_file': $!\n";
    }
    
    # warn "Detecting chromosomes from file '$vcf_file'\n\n";
    while (<DETECT>){
	$_ =~ s/(\r|\n)//g; # removing end of line characters
	last unless ($_ =~ /^\#/);
	
	if ($_ =~ /^\#\#contig/){ # filters header lines
	    # warn "$_\n"; # sleep(1);
	    $_ =~ /ID=(.+),/;
	    my $chr = $1;
	    # warn "Identified chromosome $chr\n";
	    unless (exists $chrom{$chr}){
		$chrom{$chr}++;
	    }
	}
    }
    
    foreach my $chr(keys %chrom){
	# warn "$chr\n";	
	push @chrom, $chr;
    }
    # close DETECT or warn "Failed to close filehandle DETECT: $!\n";
    return @chrom;
    
}

sub detect_strains{
    
    my $vcf_file = shift;
    my %strains; # detecting the available strains from the VCF file
    
    if ($vcf_file =~ /gz$/){
	open (STRAIN,"gunzip -c $vcf_file |") or die "Failed to open file '$vcf_file': $!\n";
    }
    else{
	open (STRAIN, $vcf_file) or die "Failed to read Input VCF file '$vcf_file': $!\n";
    }
    
    # warn "Detecting strains from file '$vcf_file'\n\n";
    while (<STRAIN>){
	$_ =~ s/(\r|\n)//g; # removing end of line characters
	next if ($_ =~ /^\#\#/);	
	if ($_ =~ /^\#CHROM/){ # header line listing all different strains
	    # warn "$_\n"; sleep(1);
	}
	last unless ($_ =~ /^\#/); # everything from now on are the SNPs themselves
	
	my @strains = split (/\t/);

    # warn "Here ist the \@strains ARRRAY:\n", join ("\n",@strains),"\n"; 
    # The first 8 fields are irrelevant: #CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  
    for(0..8){
        my $removed = shift @strains;
    }

        foreach my $index(0..$#strains){
            $strains{$index}->{strain} = $strains[$index];
            $strains{$index}->{count} = 0;
            # warn "$index\t$strains[$index]\n";
        }
    }
    
    # close STRAIN or warn "Failed to close filehandle STRAIN: $!\n";
    return %strains;

}

sub read_genome_into_memory{
    
    ## working directoy
    my $cwd = shift;
    
    ## reading in and storing the specified genome in the %chromosomes hash
    chdir ($genome_folder) or die "Can't move to $genome_folder: $!";
    warn "Now reading in and storing sequence information of the genome specified in: $genome_folder\n\n";

    my @chromosome_filenames =  <*.fa>;

    ### if there aren't any genomic files with the extension .fa we will look for files with the extension .fasta
    unless (@chromosome_filenames){
      @chromosome_filenames =  <*.fasta>;
    }

    unless (@chromosome_filenames){
      die "The specified reference genome folder $genome_folder does not contain any sequence files in FastA format (with .fa or .fasta file extensions)\n";
    }

    my $SQ_count = 0;

    foreach my $chromosome_filename (@chromosome_filenames){

	open (CHR_IN,$chromosome_filename) or die "Failed to read from sequence file $chromosome_filename $!\n";
	### first line needs to be a fastA header
	my $first_line = <CHR_IN>;
	chomp $first_line;
	$first_line =~ s/\r//;
	### Extracting chromosome name from the FastA header
	my $chromosome_name = extract_chromosome_name($first_line);
	my $sequence;

	while (<CHR_IN>){
	  chomp;
	  $_ =~ s/\r//; # removing carriage returns if present
	  if ($_ =~ /^>/){
	
	    ### storing the previous chromosome in the %chromosomes hash, only relevant for Multi-Fasta-Files (MFA)
	    if (exists $chromosomes{$chromosome_name}){
	      print "chr $chromosome_name (",length $sequence ," bp)\n";
	      die "Exiting because chromosome name already exists. Please make sure all chromosomes have a unique name!\n";
	    }
	    else {
	      if (length($sequence) == 0){
		warn "Chromosome $chromosome_name in the multi-fasta file $chromosome_filename did not contain any sequence information!\n";
	      }
	      print "chr $chromosome_name (",length $sequence ," bp)\n";
	      $chromosomes{$chromosome_name} = $sequence;
	    }
	    ### resetting the sequence variable
	    $sequence = '';
	    ### setting new chromosome name
	    $chromosome_name = extract_chromosome_name($_);
	  }
	  else{
	    $sequence .= uc$_;
	  }
	}
	
 	### Processing last chromosome of a multi Fasta File or the only entry in case of single entry FastA files

	if (exists $chromosomes{$chromosome_name}){
	    print "chr $chromosome_name (",length $sequence ," bp)\t";
	    die "Exiting because chromosome name already exists. Please make sure all chromosomes have a unique name.\n";
	}
	else{
	    if (length($sequence) == 0){
		warn "Chromosome $chromosome_name in the file $chromosome_filename did not contain any sequence information!\n";
	    }

	    print "chr $chromosome_name (",length $sequence ," bp)\n";
	    $chromosomes{$chromosome_name} = $sequence;
	}
    }
    print "\n";
    chdir $cwd or die "Failed to move to directory $cwd\n";

}

sub extract_chromosome_name {
    ## Bowtie seems to extract the first string after the inition > in the FASTA file, so we are doing this as well
    my $fasta_header = shift;
    if ($fasta_header =~ s/^>//){
	my ($chromosome_name) = split (/\s+/,$fasta_header);
	return $chromosome_name;
    }
    else{
	die "The specified chromosome ($fasta_header) file doesn't seem to be in FASTA format as required!\n";
    }
}


###################################################################


sub process_commandline{
    my $help;
    my $version;
    my $vcf_file;
    my $genome_folder;
    my $v7_MGP_file;
    my $nmasking = 1; # This is the default
    
    my $command_line = GetOptions ('help'        => \$help,
				   'versions'                    => \$version,
				   'vcf_file=s'                  => \$vcf_file,
				   'nmasking'                    => \$nmasking,
				   'genome|reference_genome=s'   => \$genome_folder,
                   'v7_VCF'                      => \$v7_MGP_file,
	);
  
    ### EXIT ON ERROR if there were errors with any of the supplied options
    unless ($command_line){
	    die "Please respecify command line options\n";
    }
    
    ### HELPFILE
    if ($help){
	    print_helpfile();
	    exit;
    }
    
    if ($version){
	print << "VERSION";
	
                          reStrainingOrder - Genome Preparation
			           version: $pipeline_version
                             Copyright 2018-23, Felix Krueger
                                 Altos Bioinformatics
                   https://github.com/FelixKrueger/reStrainingOrder

VERSION
		 ;
   	exit 1;
    }

    if (defined $vcf_file){
    	unless(-e $vcf_file){
    	    die "Input VCF file '$vcf_file' doesn't exist in the folder. Please check filenames and try again!\n\n";
    	}
        
        ### March 2022: Now also accepting the v7 combined SNP and INDEL file
		if ($v7_MGP_file){
			my $tmp_vcf_name = $vcf_file;
			$tmp_vcf_name =~ s/.*\///; # removing path information

			if($tmp_vcf_name eq 'mgp_REL2005_snps_indels.vcf.gz') {
				warn "Using v7 MGP file 'mgp_REL2005_snps_indels.vcf.gz'. This file used to be available from:\n";
				warn "ftp://ftp-mouse.sanger.ac.uk/REL-2004-v7-SNPs_Indels/mgp_REL2005_snps_indels.vcf.gz\n\n";
				warn "PLEASE NOTE:\n============\nThis file is currently not marked as Current_SNPs, so consider this approach experimental for the time being (as of 21 03 2022)\n\n";	
				sleep(1);
			}
			else{
				warn "Version 7 input file selected. Input VCF file '$vcf_file' doesn't appear to be file 'mgp_REL2005_snps_indels.vcf.gz'. If something goes wrong, I won't accept responsibility....!\n\n";
				sleep(1);
			}
		}
    }
    else{
       die "\nYou need to provide a VCF file detailing SNPs positions with '--vcf_file your.file' (e.g.: '--vcf mgp_REL2021_snps.vcf.gz'). Please respecify!\n\n";
    }

    my $strain_index;
	%strains = detect_strains($vcf_file);

	### GENOME FOLDER
    
    unless ($genome_folder){
	    warn "Reference genome folder was not specified! Please use --reference_genome </genome/folder/>\n";
    	# print_helpfile();
    	exit;
    }
    my $parent_dir = getcwd();    

    ### checking that the genome folder exists  
    unless ($genome_folder =~/\/$/){
	   $genome_folder =~ s/$/\//;
    }
    
    if (chdir $genome_folder){
    	my $absolute_genome_folder = getcwd; ## making the genome folder path absolute
    	unless ($absolute_genome_folder =~/\/$/){
    	    $absolute_genome_folder =~ s/$/\//;
    	}
    	warn "Reference genome folder provided is $genome_folder\t(absolute path is '$absolute_genome_folder)'\n\n";
    	$genome_folder = $absolute_genome_folder;
    }
    else{
	   die "Failed to move to genome folder > $genome_folder <: $!\n\n'reStraining --help' for more details\n\n";
    }
    chdir $parent_dir or die "Failed to move back to parent directory $parent_dir\n\n";

    unless (defined $genome_build){
    	$genome_build = 'GRCm39';
    }

    return ($vcf_file,$genome_folder,$nmasking,$genome_build,$v7_MGP_file);

}


sub print_helpfile{
  print <<EOF

  SYNOPSIS:

reStraining is designed to read in a variant call file from the Mouse Genomes Project (e.g. this latest
file: 'mgp_REL2021_snps.vcf.gz') and generate a new multi-strain genome version where all high-confidence
SNPs in any of the deeply sequenced straind (currently ~51) are masked by the ambiguity nucleobase 'N' (N-masking).

The resulting .fa files are ready to be indexed with your favourite aligner. Proved and tested aligners include Bowtie2,
STAR, HISAT2 and Bismark. Please note that STAR and HISAT2 may require you to disable soft-clipping.

For more information on the reStraining genome preparation process, please refer to the User Guide at:
https://github.com/FelixKrueger/reStrainingOrder/blob/master/Docs/README.md#running-restraining



  USAGE:    reStraining --vcf_file <file> --reference_genome /path/to/genome/


--vcf_file <file>             Mandatory file specifying SNP information for mouse strains from the Mouse Genomes Project
                              (https://www.mousegenomes.org/). The file used and approved for genome build GRCm39 is 
                              called 'mgp_REL2021_snps.vcf.gz'. Please note that future versions of this file or entirely
                              different VCF files might not work out-of-the-box but may require some tweaking.
                              SNP calls are read from the VCF files, and high confidence SNPs are written into
                              a folder in the current working directory called SNPs_directory/chr<chromosome>.txt,
                              in the following format:

                              Chromosome  Position   REF   ALT   [0 or 1 for each strain]

                              Where 0 means REF base, and 1 means ALT base.

--v7_VCF                      This will use the file 'mgp_REL2005_snps_indels.vcf.gz' for backward compativility.
                              This file contains both SNP and INDEL information, but INDELs are skipped. The v7
                              file contains a number of additional strains, and was released in May2020.

--reference_genome <PATH>     The path to the reference genome, typically the strain 'Black6' (C57BL/6J). The mouse genomes
                              project uses the Ensembl genome build GRCm39, and will likely fail in conjunction with genomes
                              from NCBI or UCSF. reStraining expects one or more FastA files in this folder
                              (file extension: .fa or .fasta). 
                              
                              The mouse genome files can be downloaded from Ensembl via this command:
                              wget ftp://ftp.ensembl.org/pub/release-109/fasta/mus_musculus/dna/*dna.chromosome.*


--help                        Displays this help information and exits.

--version                     Displays version information and exits.


                                                             Last modified: 11 04 2023

EOF
    ;
  exit 1;
}