<commit_before>#!/usr/bin/env perl

use FindBin;
use lib "$FindBin::RealBin/../lib";

use strict;
use warnings;
use Bio::Perl;
use Data::Dumper;
use Getopt::Long;
use File::Basename;
use File::Spec;
use threads;
use Thread::Queue;
use Schedule::SGELK;
#use ExtUtils::Command;


my ($name,$scriptsdir,$suffix)=fileparse($0);
$scriptsdir=File::Spec->rel2abs($scriptsdir);

# Logging
my $logmsgFh;
sub logmsg {
  local $0=basename $0;
  my $FH = *STDOUT; 
  my $msg="$0: ".(caller(1))[3].": @_\n";

  # print the message to the logfile if it's not the same as stdout
  #print join("\t","===",fileno($logmsgFh),fileno($FH))."\n";
  if(defined($logmsgFh) && fileno($logmsgFh) > 0){
    print $logmsgFh $msg;
  }
  print $FH $msg;
}
local $SIG{'__DIE__'} = sub { local $0=basename $0; my $e = $_[0]; $e =~ s/(at [^\s]+? line \d+\.$)/\nStopped $1/; my $msg="$0: ".(caller(1))[3].": ".$e; logmsg($msg);  die(); };
END{
  close $logmsgFh if(defined($logmsgFh) && fileno($logmsgFh) > 0);
}

my $sge=Schedule::SGELK->new(-verbose=>1,-numnodes=>20,-numcpus=>8);
exit(main());

sub main{
  # start with the settings that are on by default, and which can be turned off by, e.g., --noclean
  my $settings={trees=>1,msa=>1, matrix=>1, clean=>1};
  GetOptions($settings,qw(ref=s bamdir=s logdir=s vcfdir=s tmpdir=s readsdir=s asmdir=s msadir=s help numcpus=s numnodes=i allowedFlanking=s keep min_alt_frac=s min_coverage=i trees! queue=s qsubxopts=s msa! matrix! mapper=s snpcaller=s msa-creation=s clean)) or die $!;
  # Lyve-SET
  $$settings{allowedFlanking}||=0;
  $$settings{keep}||=0;
  $$settings{min_alt_frac}||=0.75;
  $$settings{min_coverage}||=10;
  # modules defaults
  $$settings{mapper}||="smalt";
  $$settings{snpcaller}||="varscan";
  $$settings{'msa-creation'}||="lyve-set-lowmem";
  # queue stuff
  $$settings{vcfToAlignment_xopts}||="-l mem_free=100G -q highmem.q";
  $$settings{numcpus}||=1;
  $$settings{numnodes}||=20;
  $$settings{qsubxopts}||="";
  $$settings{queue}||="";
  # Some things need to just be lowercase to make things easier downstream
  $$settings{$_}=lc($$settings{$_}) for(qw(msa-creation snpcaller mapper));

  # A warning about --noclean
  if(!$$settings{clean}){
    die "Warning: the --noclean option has been removed in Lyve-SET v0.9.2.  You should clean your reads before using Lyve-SET.";
  }

  ##########################################################
  ### Other defaults: reference genome; default directories#
  ##########################################################
  my $project=shift(@ARGV) || '.'; # by default the user is in the project directory
  die usage($settings) if($$settings{help});
  logmsg "Project was defined as $project"; # No period at the end of sentence to avoid ambiguous '..' in the logmsg
  # Checks to make sure that this is a project
  #system("set_manage.pl '$project'"); die if $?;
  # Set the defaults if they aren't set already
  $$settings{ref}||="$project/reference/reference.fasta";
  # Check SET directories' existence and set their defaults
  for my $param (qw(vcfdir bamdir msadir readsdir tmpdir asmdir logdir)){
    my $b=$param;
    $b=~s/dir$//;  # e.g., vcfdir => vcf
    $$settings{$param}||="$project/$b";
    die "ERROR: Could not find $param under $$settings{$param}/\n  mkdir $$settings{$param} to resolve this problem.\n".usage($settings) if(!-d $$settings{$param});
    $$settings{$param}=File::Spec->rel2abs($$settings{$param});
  }

  # SGE params
  $$settings{workingdir}=$$settings{logdir};
  for (qw(workingdir numnodes numcpus keep qsubxopts queue)){
    $sge->set($_,$$settings{$_}) if($$settings{$_});
  }

  # open the log file now that we know the logdir is there
  open($logmsgFh,">$$settings{logdir}/launch_set.log") or die "ERROR: could not open log file $$settings{logdir}/launch_set.log";
  logmsg "======\nOpened logfile at $$settings{logdir}/launch_set.log\n=====";

  # Check the reference parameter
  die "ERROR: reference file was not given\n".usage($settings) if(!defined($$settings{ref}));
  die "ERROR: Could not find the reference file at $$settings{ref}\n".usage($settings) if(!-f $$settings{ref});
  my $ref=$$settings{ref};

  simulateReads($settings);
  indexReference($ref,$settings);
  mapReads($ref,$$settings{readsdir},$$settings{bamdir},$settings);
  variantCalls($ref,$$settings{bamdir},$$settings{vcfdir},$settings);

  my $pooled;
  if($$settings{matrix}){
    $pooled=variantsToMatrix($ref,$$settings{bamdir},$$settings{vcfdir},$$settings{msadir},$settings);
  } else {
    logmsg "The matrix was not requested; wrapping up";
    return 0;
  }

  if($$settings{msa}){
    pooledToAlignment($pooled,$settings);
  } else {
    logmsg "The alignment was not requested; wrapping up";
    return 0;
  }

  if($$settings{trees}){
    logmsg "Launching set_processMsa.pl";
    $sge->pleaseExecute("set_processMsa.pl --auto --msaDir '$$settings{msadir}' --numcpus $$settings{numcpus} 2>&1 | tee $$settings{logdir}/set_processMsa.log ",{numcpus=>$$settings{numcpus},jobname=>"set_processMsa.pl"});
    $sge->wrapItUp();
  } else {
    logmsg "The phylogeny was not requested; wrapping up";
    return 0;
  }

  return 0;
}

sub simulateReads{
  my($settings)=@_;
  logmsg "Simulating reads from any assemblies";
  my @asm=glob("$$settings{asmdir}/*.{fasta,fna,ffn,fa,mfa,mfasta}");
  return if(!@asm);

  my $exec=`which wgsim`;
  die "ERROR: Could not find wgsim which is part of samtools" if $?;
  chomp($exec);

  for my $asm (@asm){
    my $b=fileparse $asm;
    my $outFastq="$$settings{readsdir}/$b.wgsim.fastq";
    my $outGz="$outFastq.gz";
    my $out1="$$settings{tmpdir}/$b.1.fastq";
    my $out2="$$settings{tmpdir}/$b.2.fastq";
    next if(-e $outGz && !-e $outFastq); 
    logmsg "I did not find $outFastq simulated reads. Simulating now with $exec.";
    
    # Four successive jobs that depend on each other:
    #   1. Simulate into paired ends
    #   2. Shuffle
    #   3. Gzip
    #   4. Remove simulated unshuffled reads (depends on the shuffle step and not gzip)
    $sge->pleaseExecute("$exec -d 400 -N 100000 -1 250 -2 250 -e 0.0 -r 0.0 -R 0.0000 -h '$asm' $out1 $out2",{jobname=>"wgsim$b",numcpus=>1});
    $sge->pleaseExecute("run_assembly_shuffleReads.pl $out1 $out2 > $outFastq",{jobname=>"shuffle$b",qsubxopts=>"-hold_jid wgsim$b",numcpus=>1});
    $sge->pleaseExecute("gzip -v '$outFastq' 2>&1",{jobname=>"gzip$b",qsubxopts=>"-hold_jid shuffle$b",numcpus=>1});
    $sge->pleaseExecute("rm -v $out1 $out2",{jobname=>"rmTmp$b",qsubxopts=>"-hold_jid shuffle$b",numcpus=>1});
  }
  $sge->wrapItUp();
}

sub indexReference{
  my($ref,$settings)=@_;
  logmsg "Indexing the reference for read mapping";

  # sanity check: see if the reference has reserved characters
  my $in=Bio::SeqIO->new(-file=>$ref);
  while(my $seq=$in->next_seq){
    my $defline=$seq->id." ".$seq->desc;
    die "Offending character found in the defline\n $defline\n $1" if($defline=~/([\:])/);
  }

  logmsg "Indexing with $$settings{mapper}";
  if($$settings{mapper} eq 'smalt'){
    return $ref if(-e "$ref.sma" && -e "$ref.smi");
    system("smalt index -k 13 -s 2 $ref $ref 2>&1");
    die if $?;
  } elsif($$settings{mapper} eq 'snap'){
    return $ref if(-d "$ref.snap" && -e "$ref.snap/GenomeIndex");
    system("snap index $ref $ref.snap -s 16");
    die if $?;
  }
  return $ref;
}

sub mapReads{
  my($ref,$readsdir,$bamdir,$settings)=@_;
  logmsg "Mapping reads with $$settings{mapper}";<commit_msg>integrated masking for phages<commit_after>#!/usr/bin/env perl

use FindBin;
use lib "$FindBin::RealBin/../lib";

use strict;
use warnings;
use Bio::Perl;
use Data::Dumper;
use Getopt::Long;
use File::Basename;
use File::Spec;
use threads;
use Thread::Queue;
use Schedule::SGELK;
#use ExtUtils::Command;


my ($name,$scriptsdir,$suffix)=fileparse($0);
$scriptsdir=File::Spec->rel2abs($scriptsdir);

# Logging
my $logmsgFh;
sub logmsg {
  local $0=basename $0;
  my $FH = *STDOUT; 
  my $msg="$0: ".(caller(1))[3].": @_\n";

  # print the message to the logfile if it's not the same as stdout
  #print join("\t","===",fileno($logmsgFh),fileno($FH))."\n";
  if(defined($logmsgFh) && fileno($logmsgFh) > 0){
    print $logmsgFh $msg;
  }
  print $FH $msg;
}
local $SIG{'__DIE__'} = sub { local $0=basename $0; my $e = $_[0]; $e =~ s/(at [^\s]+? line \d+\.$)/\nStopped $1/; my $msg="$0: ".(caller(1))[3].": ".$e; logmsg($msg);  die(); };
END{
  close $logmsgFh if(defined($logmsgFh) && fileno($logmsgFh) > 0);
}

my $sge=Schedule::SGELK->new(-verbose=>1,-numnodes=>20,-numcpus=>8);
exit(main());

sub main{
  # start with the settings that are on by default, and which can be turned off by, e.g., --noclean
  my $settings={trees=>1,msa=>1, matrix=>1, clean=>1};
  GetOptions($settings,qw(ref=s bamdir=s logdir=s vcfdir=s tmpdir=s readsdir=s asmdir=s msadir=s help numcpus=s numnodes=i allowedFlanking=s keep min_alt_frac=s min_coverage=i trees! queue=s qsubxopts=s msa! matrix! mapper=s snpcaller=s msa-creation=s clean)) or die $!;
  # Lyve-SET
  $$settings{allowedFlanking}||=0;
  $$settings{keep}||=0;
  $$settings{min_alt_frac}||=0.75;
  $$settings{min_coverage}||=10;
  # modules defaults
  $$settings{mapper}||="smalt";
  $$settings{snpcaller}||="varscan";
  $$settings{'msa-creation'}||="lyve-set-lowmem";
  # queue stuff
  $$settings{vcfToAlignment_xopts}||="-l mem_free=100G -q highmem.q";
  $$settings{numcpus}||=1;
  $$settings{numnodes}||=20;
  $$settings{qsubxopts}||="";
  $$settings{queue}||="";
  # Some things need to just be lowercase to make things easier downstream
  $$settings{$_}=lc($$settings{$_}) for(qw(msa-creation snpcaller mapper));

  # A warning about --noclean
  if(!$$settings{clean}){
    die "Warning: the --noclean option has been removed in Lyve-SET v0.9.2.  You should clean your reads before using Lyve-SET.";
  }

  ##########################################################
  ### Other defaults: reference genome; default directories#
  ##########################################################
  my $project=shift(@ARGV) || '.'; # by default the user is in the project directory
  die usage($settings) if($$settings{help});
  logmsg "Project was defined as $project"; # No period at the end of sentence to avoid ambiguous '..' in the logmsg
  # Checks to make sure that this is a project
  #system("set_manage.pl '$project'"); die if $?;
  # Set the defaults if they aren't set already
  $$settings{ref}||="$project/reference/reference.fasta";
  # Check SET directories' existence and set their defaults
  for my $param (qw(vcfdir bamdir msadir readsdir tmpdir asmdir logdir)){
    my $b=$param;
    $b=~s/dir$//;  # e.g., vcfdir => vcf
    $$settings{$param}||="$project/$b";
    die "ERROR: Could not find $param under $$settings{$param}/\n  mkdir $$settings{$param} to resolve this problem.\n".usage($settings) if(!-d $$settings{$param});
    $$settings{$param}=File::Spec->rel2abs($$settings{$param});
  }

  # SGE params
  $$settings{workingdir}=$$settings{logdir};
  for (qw(workingdir numnodes numcpus keep qsubxopts queue)){
    $sge->set($_,$$settings{$_}) if($$settings{$_});
  }

  # open the log file now that we know the logdir is there
  open($logmsgFh,">$$settings{logdir}/launch_set.log") or die "ERROR: could not open log file $$settings{logdir}/launch_set.log";
  logmsg "======\nOpened logfile at $$settings{logdir}/launch_set.log\n=====";

  # Check the reference parameter
  die "ERROR: reference file was not given\n".usage($settings) if(!defined($$settings{ref}));
  die "ERROR: Could not find the reference file at $$settings{ref}\n".usage($settings) if(!-f $$settings{ref});
  my $ref=$$settings{ref};

  simulateReads($settings);
  maskReference($ref,$settings);
  indexReference($ref,$settings);
  mapReads($ref,$$settings{readsdir},$$settings{bamdir},$settings);
  variantCalls($ref,$$settings{bamdir},$$settings{vcfdir},$settings);

  my $pooled;
  if($$settings{matrix}){
    $pooled=variantsToMatrix($ref,$$settings{bamdir},$$settings{vcfdir},$$settings{msadir},$settings);
  } else {
    logmsg "The matrix was not requested; wrapping up";
    return 0;
  }

  if($$settings{msa}){
    pooledToAlignment($pooled,$settings);
  } else {
    logmsg "The alignment was not requested; wrapping up";
    return 0;
  }

  if($$settings{trees}){
    logmsg "Launching set_processMsa.pl";
    $sge->pleaseExecute("set_processMsa.pl --auto --msaDir '$$settings{msadir}' --numcpus $$settings{numcpus} 2>&1 | tee $$settings{logdir}/set_processMsa.log ",{numcpus=>$$settings{numcpus},jobname=>"set_processMsa.pl"});
    $sge->wrapItUp();
  } else {
    logmsg "The phylogeny was not requested; wrapping up";
    return 0;
  }

  return 0;
}

sub simulateReads{
  my($settings)=@_;
  logmsg "Simulating reads from any assemblies";
  my @asm=glob("$$settings{asmdir}/*.{fasta,fna,ffn,fa,mfa,mfasta}");
  return if(!@asm);

  my $exec=`which wgsim`;
  die "ERROR: Could not find wgsim which is part of samtools" if $?;
  chomp($exec);

  for my $asm (@asm){
    my $b=fileparse $asm;
    my $outFastq="$$settings{readsdir}/$b.wgsim.fastq";
    my $outGz="$outFastq.gz";
    my $out1="$$settings{tmpdir}/$b.1.fastq";
    my $out2="$$settings{tmpdir}/$b.2.fastq";
    next if(-e $outGz && !-e $outFastq); 
    logmsg "I did not find $outFastq simulated reads. Simulating now with $exec.";
    
    # Four successive jobs that depend on each other:
    #   1. Simulate into paired ends
    #   2. Shuffle
    #   3. Gzip
    #   4. Remove simulated unshuffled reads (depends on the shuffle step and not gzip)
    $sge->pleaseExecute("$exec -d 400 -N 100000 -1 250 -2 250 -e 0.0 -r 0.0 -R 0.0000 -h '$asm' $out1 $out2",{jobname=>"wgsim$b",numcpus=>1});
    $sge->pleaseExecute("run_assembly_shuffleReads.pl $out1 $out2 > $outFastq",{jobname=>"shuffle$b",qsubxopts=>"-hold_jid wgsim$b",numcpus=>1});
    $sge->pleaseExecute("gzip -v '$outFastq' 2>&1",{jobname=>"gzip$b",qsubxopts=>"-hold_jid shuffle$b",numcpus=>1});
    $sge->pleaseExecute("rm -v $out1 $out2",{jobname=>"rmTmp$b",qsubxopts=>"-hold_jid shuffle$b",numcpus=>1});
  }
  $sge->wrapItUp();
}

sub maskReference{
  my($ref,$settings)=@_;
  
  # Don't do this if there is already an unmasked file, aka if it's already done
  return if(-e "$ref.unmasked"); 

  logmsg "Masking the reference genome for phages";
  system("set_findPhages.pl --numcpus $$settings{numcpus} $ref > $$settings{tempdir}/reference.masked.fasta");
  die if $?;

  system("mv -v $ref $ref.unmasked"); die if $?;
  system("mv -v $$settings{tempdir}/reference.masked.fasta $ref");
  die if $?;
}

sub indexReference{
  my($ref,$settings)=@_;
  logmsg "Indexing the reference for read mapping";

  # sanity check: see if the reference has reserved characters
  my $in=Bio::SeqIO->new(-file=>$ref);
  while(my $seq=$in->next_seq){
    my $defline=$seq->id." ".$seq->desc;
    die "Offending character found in the defline\n $defline\n $1" if($defline=~/([\:])/);
  }

  logmsg "Indexing with $$settings{mapper}";
  if($$settings{mapper} eq 'smalt'){
    return $ref if(-e "$ref.sma" && -e "$ref.smi");
    system("smalt index -k 13 -s 2 $ref $ref 2>&1");
    die if $?;
  } elsif($$settings{mapper} eq 'snap'){
    return $ref if(-d "$ref.snap" && -e "$ref.snap/GenomeIndex");
    system("snap index $ref $ref.snap -s 16");
    die if $?;
  }
  return $ref;
}

sub mapReads{
  my($ref,$readsdir,$bamdir,$settings)=@_;
  logmsg "Mapping reads with $$settings{mapper}";
{'fraction_non_alphanumeric': 0.20196499446938643, 'fraction_numerical': 0.013143340490597956, 'mean_word_length': 4.715879509111193, 'pattern_counts': {'":': 4, '<': 3, '<?xml version=': 0, '>': 92, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '197468', 'n_tokens_mistral': 5845, 'n_tokens_neox': 5580, 'n_words': 1354}