<commit_before>#!/usr/bin/perl

=head1 NAME

 intersect_genome_snps.pl
This script extracts the SNP intersection from 2 subsets of mpileup.intersect files, preferably generated with extract_snp_freq.pl.
A report will be generated with the number of similar SNPs removed from the target mpileup file, and the number of non-matchiing removed SNPs.
The criteria for extracting SNPs from the target file is a position which is unique, and does not have a SNP in the reference file.
This should hopefully remove SNPs that are shared between 2 genomes, and will leave potential SNPs that came from the target intersect file   

=cut

=head1 SYPNOSIS

 intersect_genome_snps.pl [-h] --target <input_mpileup_file> --intersect <input_intersecting_mpileup file(s , comma separated) >  -o <output_file>

=head1 DESCRIPTION

 This script reads 2 files generated by extract_snp_freq.pl: subset of mpileup file from a target genome and a reference genomes with sought SNPs, and anohter subset mpileup file of the same target genome and SNPs that should be filtered (for example a target file of an inbred genome filtered from SNPs from a close relative, and a second SNP file from anohter potential ancestor(s) that should be subtracted from the first set)  
 and a varscan file of a reference genome ( for example a possible donor for the target SNPs) .

Both files should have an extra column (#7) with the allele frequency from the target genome.
This allows to plot the SNPs of the target genome X reference genome minus the SNPs from a second reference genome.
This subtracting tage can be run multiple times, if you wish to look at each genome pair separately, each time removing non-specific SNPs of the reference species compared to one other genome.


=head2 I<Flags:>

=over

=item --target | t

B<input mpileup file>             input mpileup file of the target genome (mandatory)

=item --intersect | i

B<input mpileup file (s)>           input mpileup file(s) of a second genome (or more, comma separated lists are supported), to be removed from file #1 above 

    
=item -o

B<output_file>            output file (mandatory)


=item -h

B<help>                   print this help doc

=back

=head1 AUTHOR

<NAME><<EMAIL>>

=cut



use strict;
use warnings;
use File::Slurp;

use Getopt::Long;
use Pod::Usage;
use Carp qw /croak/;

my ( $target_mpileup, @intersect_mpileup, $out, $help); 

GetOptions (
    "target|t=s"    => \$target_mpileup,    # string
    "intersect|i=s" => \@intersect_mpileup, # string
    "out=s"    => \$out,
    "help"  => \$help)   # flag
    or  pod2usage(-verbose  => 2);

if ($help || !$target_mpileup || !@intersect_mpileup || !$out)  { pod2usage(-verbose  => 2); } 

# list file is varscan with chr. in column 1 and position in clumn 2

open (TARGET, "<$target_mpileup") || die "Cannot open target mpileup file $target_mpileup.\n";

@intersect_mpileup = split(/,/,join(',',@intersect_mpileup));
my $err = $out . ".err" ;
my %intersect;

foreach my $file (@intersect_mpileup) {
    #remove trailing and leading spaces
    $file =~ s/^\s(.?)\s*/$1/;
    open (INTERSECT, "<$file") || croak "Cannot open intersect mpileup file $file.\n";
    print STDERR "Intersect file: $file \n\n";

    # parse the intersecting file 
    while ( <INTERSECT> ) {
	my $var_line = $_;
	chomp($var_line);
	my @var = split ("\t", $var_line);
	my $chr = $var[0];
	my $pos = $var[1];
	my $var_allele = $var[6];
	my $existing_allele = $intersect{$chr . ":" . $pos } ;
	my $allele_list = $var_allele;
	if ($existing_allele) { $allele_list = join("\t" , ( $existing_allele, $var_allele ) ) ; } 
	$intersect{$chr . ":" . $pos } = $allele_list;
	print STDERR " Intersect : $chr \t $pos \t allele = $allele_list \n";
    }
}


#print STDERR " Target : $chr \t $pos \t allele = $target_allele ";                                                                                                 

while ( <TARGET> ) {
    my $line = $_;
    chomp($line);
    my @tar = split ("\t", $line);
    my $chr = $tar[0];
    my $pos = $tar[1];
    my $target_allele = $tar[6];

    my $var_allele = $intersect{$chr . ":" . $pos };
    # if an allele exists in the same position in the intersect file, discard it. 
    # print in the error file both alleles for counting which alleles are identical and which are different 
    if ($var_allele) {
	write_file($err, {append => 1 } , ( ( join ("\t" , ( $chr, $pos, $target_allele, $var_allele) ) ) , "\n" ) ) ;
    } else {
	print STDERR " Target : $chr \t $pos \t allele = $target_allele \n";                                                                                        
	write_file($out, { append => 1} , ( $line,"\n")  );
    }
}
<commit_msg>look at the mpileup of the intersecting genomes instead of varscan files, for refining unique SNPs in teh target genome originating from the reference genome<commit_after>#!/usr/bin/perl

=head1 NAME

 intersect_genome_snps.pl
This script extracts the SNP intersection from 2 subsets of mpileup.intersect files, preferably generated with extract_snp_freq.pl.
A report will be generated with the number of similar SNPs removed from the target mpileup file, and the number of non-matchiing removed SNPs.
The criteria for extracting SNPs from the target file is a position which is unique, and does not have a SNP in the reference file.
This should hopefully remove SNPs that are shared between 2 genomes, and will leave potential SNPs that came from the target intersect file   

=cut

=head1 SYPNOSIS

 intersect_genome_snps.pl [-h] --target <input_mpileup_file> --intersect <input_intersecting_mpileup file(s , comma separated) >  -o <output_file>

=head1 DESCRIPTION

 This script reads 2 files generated by extract_snp_freq.pl: subset of mpileup file from a target genome and a reference genomes with sought SNPs, and anohter subset mpileup file of the same target genome and SNPs that should be filtered (for example a target file of an inbred genome filtered from SNPs from a close relative, and a second SNP file from anohter potential ancestor(s) that should be subtracted from the first set)  
 and a varscan file of a reference genome ( for example a possible donor for the target SNPs) .

Both files should have an extra column (#7) with the allele frequency from the target genome.
This allows to plot the SNPs of the target genome X reference genome minus the SNPs from a second reference genome.
This subtracting tage can be run multiple times, if you wish to look at each genome pair separately, each time removing non-specific SNPs of the reference species compared to one other genome.


=head2 I<Flags:>

=over

=item --target | t

B<input mpileup file>             input mpileup file of the target genome (mandatory)

=item --mpileup | m

B<input mpileup file (s)>           input mpileup file(s) of a second genome (or more, comma separated lists are supported), to be removed from file #1 above 

=item --cov | c

B<coverage cutoff> cutoff for ignoring SNPs in positions that are not covered well in the intersecting genome(s) 
   
=item -o

B<output_file>            output file (mandatory)


=item -h

B<help>                   print this help doc

=back

=head1 AUTHOR

<NAME><<EMAIL>>

=cut



use strict;
use warnings;
use File::Slurp;

use Getopt::Long;
use Pod::Usage;
use Carp qw /croak/;

my ( $target_mpileup, @intersect_freq, @intersect_mpileup, $out, $help); 
my $intersect_coverage = 1;

GetOptions (
    "target|t=s"    => \$target_mpileup,    # string
    "intersect|i=s" => \@intersect_freq, 
    "mpileup|m=s"   => \@intersect_mpileup,
    "cov|c=i"       => \$intersect_coverage,
    "out|o=s"    => \$out,
    "help"  => \$help)   # flag
    or  pod2usage(-verbose  => 2);

if ($help || !$target_mpileup || !@intersect_mpileup  || !$out)  { pod2usage(-verbose  => 2); } 


open (TARGET, "<$target_mpileup") || die "Cannot open target mpileup file $target_mpileup.\n";

@intersect_freq = split(/,/,join(',',@intersect_freq));
@intersect_mpileup = split(/,/,join(',',@intersect_mpileup));
my $err = $out . ".err" ;

## this was when we lookied ony at varscan files. Now looking at mpileup because it does not really matter if the intersecting genomes have shared SNPs
### We want to look at coverage threshold, and if the reference SNP is unique to that genome in each position (i.e. if any oth the intersecting genomes has the same nucleotide in that position, the SNP is no longer considered unique 

#############################
#my %intersect;

#foreach my $file (@intersect_freq) {
    #remove trailing and leading spaces
#    $file =~ s/^\s(.?)\s*/$1/;
#    open (INTERSECT, "<$file") || croak "Cannot open intersect mpileup file $file.\n";
#    print STDERR "Intersect file: $file \n\n";

    # parse the intersecting file 
#    while ( <INTERSECT> ) {#
#	my $var_line = $_;#
#	chomp($var_line);
#	my @var = split ("\t", $var_line);
#	my $chr = $var[0];
#	my $pos = $var[1];
#	my $var_allele = $var[6];
#	my $existing_allele = $intersect{$chr . ":" . $pos } ;
#	my $allele_list = $var_allele;
#	if ($existing_allele) { $allele_list = join("\t" , ( $existing_allele, $var_allele ) ) ; } 
#	$intersect{$chr . ":" . $pos } = $allele_list;
#	print STDERR " Intersect : $chr \t $pos \t allele = $allele_list \n";
 #   }
#}

################################################################################
#print STDERR " Target : $chr \t $pos \t allele = $target_allele ";
my %error;
my %target;
my ($chr , $prev_chr, $pos);

print STDERR " TARGET file $target_mpileup \n\n";

#parse the target mpileup subset file containing the candidate SNPs from the genome of interest
while ( <TARGET> ) {
    my $line = $_;
    chomp($line);
    my @tar = split ("\t", $line);
    $chr = $tar[0];
    $pos = $tar[1];
    my $target_allele = $tar[6];
    no warnings 'uninitialized' ;
    if ( $chr ne $prev_chr ) { print STDERR "$chr\n"; }
    $prev_chr = $chr;

    #############Write to out files only after checking the intersecting genomes 


    ##my $var_allele = $intersect{$chr . ":" . $pos };
    # if an allele exists in the same position in the intersect file, discard it. 
    # print in the error file both alleles for counting which alleles are identical and which are different 
    ##if ($var_allele) {
##	$error{ $chr . ":" . $pos }->{target_allele} = $target_allele;
##	$error{ $chr . ":" . $pos }->{var_allele} = $var_allele;
	
	#write_file($err, {append => 1 } , ( ( join ("\t" , ( $chr, $pos, $target_allele, $var_allele) ) ) , "\n" ) ) ;
  ##  } else {
##    print STDERR " Target : $chr \t $pos \t allele = $target_allele \n";                                                                              
    $target{ $chr . ":" . $pos }->{line}=  $line ;
    $target{ $chr . ":" . $pos }->{target_allele}=  $target_allele ;

	#write_file($out, { append => 1} , ( $line,"\n")  );
    ##}
}

#look into the mpileup files of the intersecting genomes for surther filtering by
#looking at minimum coverage required, and if the allele in the target file is unique (must be unique when compared to all intersecting files)

foreach my $mpileup_file (@intersect_mpileup) {
    open (MPILEUP, "<$mpileup_file") || croak "Cannot open intersect mpileup file $mpileup_file.\n";
    my %mpile;
    print STDERR "Finding target-mpileup file matches in file $mpileup_file \n";
    #############
    while ( <MPILEUP> ) {
	my $line = $_;
	chomp($line);
	my @mp = split ("\t", $line);
	$chr = $mp[0];
	if ( $chr ne $prev_chr ) { print STDERR "$chr\n"; }  ;
	$prev_chr = $chr;

	$pos = $mp[1];
	my $allele = $mp[2];
	my $coverage = $mp[3];
	if (exists $target{ $chr . ":" . $pos } ) {
	    my $target_allele = $target{ $chr . ":" . $pos }->{target_allele} ;
	    if ( $target_allele && ( $coverage < $intersect_coverage || $allele eq  $target_allele ) ) {
		delete $target{ $chr . ":" . $pos } ;
		write_file( $err , { append => 1 }, ( join("\t" , ( $chr, $pos, $mpileup_file , $target_allele , $allele , $coverage ) ) , "\n" ) ) ;
	    }
	}
    }
}
    #############
###print STDERR "Finding target-mpileup file matches in file $mpileup_file \n";
    
#    foreach my $key ( sort keys %target ) { #
#	( $chr , $pos ) = split (":" , $key) ;
#	if ( $chr ne $prev_chr ) { print STDERR "$chr\n"; }
#	$prev_chr = $chr;

	#my $match = $mpile{ $key } ;
	#my $allele = $match->{allele};
	#my $coverage = $match->{coverage};
#	my $target_allele = $target{ $key }->{target_allele};
#	my $match_string = $chr . ":" . $pos . "-" . $pos ;
#	my $fasta = "/export/home/naama/ITAG2.3/ITAG_pre2.3_genomic.fasta";
#	my $mline = `samtools mpileup -r $match_string -C 50  -f $fasta  $mpileup_file ` ;# || die "Cannot execute command samtools mpileup \n\n";
#	if ($mline) {
#	    my @fields = split("\t" , $mline) ;
#	    my $allele = $fields[2];
#	    my $coverage = $fields[3];
	    # if the coverage for this position in the mpileup file is lower than the cutoff coverage ( -c )
	    # and/or the target allele is not unique ( similar to the mpileup allele in the intersect file )
	    # remove this position from the target file 
#	    if ( $coverage < $intersect_coverage || $allele eq  $target_allele ) {
#		print STDERR "DELETING $key : coverage = $coverage , allele = $allele target = $target_allele\n";
#		write_file( $err , { append => 1 }, ( join("\t" , ( $chr, $pos, $mpileup_file , $target_allele , $allele , $coverage ) ) , "\n" ) ) ; 
#		delete  $target{ $key } ;
#	    }
#	}
 #   }
#}

# need to write to output file at the end, after looping through all genomes
foreach my $key (sort  keys  %target ) { 
    write_file($out, { append => 1} , ( $target{ $key }->{line} ,"\n")  );
}


{'fraction_non_alphanumeric': 0.13780684991915332, 'fraction_numerical': 0.0047038071439071, 'mean_word_length': 3.2535167239762424, 'pattern_counts': {'":': 15, '<': 38, '<?xml version=': 0, '>': 61, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3915737', 'n_tokens_mistral': 4710, 'n_tokens_neox': 4439, 'n_words': 1623}