#!/usr/bin/perl
# Author: Martin Dahlo
#
# Usage: perl polyAlign.pl -a <aligner> -i <index/reference genome> -o <output directory> -u <user file> -r <read file(s for bioscope mate pair, comma sep.)> -q <qual file(s for bioscope mate pair, comma sep.)> (-j <jump database file for mosaik>)
# ex.
# perl 

use warnings;
use strict;
use Cwd;
use File::Path;
use Getopt::Std;
use Time::Local;

=pod
This script will take the arguments that most align programs want, and format
them in the correct way for each align program. No more hassle with reading
the manual for each program and trying to figure out how to make the program
run. Only to be used with default parameters. If you are skilled enough to
know how to tweak aligner parameters, you won't need this script in the first
place :)

Pseudo code:
* Get all the arguments
* Write a submission file for SLURM using the given arguments
(* Submit the file to the queue system)
=cut


### global settings ###
my $scripts = "/bubo/home/h18/dahlo/glob/work/uppmaxScripts"; # why is this not in the path in perl, when it is in my path outside perl?

# define messages
my $usage = "Usage: perl polyAlign.pl -a <aligner> -i <index/reference genome> -o <output directory> -u <user file> -r <read file(s for bioscope mate pair, comma sep.)> -q <qual file(s for bioscope mate pair, comma sep.)> (-j <jump database file for mosaik>)\n";
my $userFileExp = "ERROR: User file (-u) not specified.\n\nPlease create your own user file with information that is needed to write a submission file to the SLURM system on the cluster.\n\nExample of a user file:\n\n# Uppmax project id\nproj=b2010033\n\n# Type of job (node/core?)\ntype=node\n\n# Output file (WARNING, relative to the output directory!)\nout=run.out\n\n# Error file (WARNING, relative to the output directory!)\nerr=run.err\n\n# Email-type\nemailType=ALL\n\n# Email address\nemailAdd=martin.dahlo\@scilifelab.uu.se\n\n# Requested time\nreqTime=24:00:00\n\n# Job name\njobName=polyAlign\n";
my $validAligners = "ERROR: Aligner (-a) not specified!\nValid aligners are:\nbowtie\t\tBowtie aligner\nbwa\t\tBwa aligner\nbioscopeSE\tBioscope single end (just ordinary) run\nbioscopeMP\tBioscope mate pair run\nbioscopeWTSE\tBioscope whole transcriptome run, single end\nbioscopePE\tBioscope paired end run\nmaq\t\tMaq aligner\nmosaik\t\tMosaik aligner\n\n";


# create a timestamp
my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst)=localtime(time);
$year += 1900; # to make it y2k compatible (seriously, how lame is that? :)  )
# add zero to all one-digit numbers  (9 -> 09, 10 -> 10)
$year = sprintf("%.2d", $year);
$mon = sprintf("%.2d", $mon+1); # adjust to real month, not array index
$mday = sprintf("%.2d", $mday);
$hour = sprintf("%.2d", $hour);
$min = sprintf("%.2d", $min);
$sec = sprintf("%.2d", $sec);
my $stamp = "$year$mon$mday$hour$min$sec";
#print "$stamp\n$year\t$mon\t$mday\t$hour\t$min\t$sec\n";

# get number of cores on the node
open(NR,"cat /proc/cpuinfo | grep processor | wc -l |") || die "Failed: $!\n";
my $cores = <NR>;
chomp($cores);


# get common options
our ($opt_a,$opt_u,$opt_o,$opt_i,$opt_r,$opt_q,$opt_j,$opt_r1,$opt_r2,$opt_q1,$opt_q2);
getopt('auoirqj');
my $aligner = $opt_a or die $validAligners;
my $userFile = $opt_u or die $userFileExp;
$userFile = absPath($userFile);
my $out = $opt_o or die "ERROR: Output directiory (-o) not specified!\n";
$out = absPath($out);
$out = remSlash($out);
my $reads = $opt_r or die "ERROR: Read file (-r) not specified\n";
$reads = absPath($reads);

# extract file name and path
my @filename = split(/\//,$reads);
my $fileName = pop(@filename);
my $readPath = join("/",@filename);

my $qual = $opt_q or die "ERROR: Quality file (-q) not specified\n";
$qual = absPath($qual);


# open file handles
mkpath($out);
open SBATCH, ">", "$out/$stamp.sbatch" or die $!;
open USER, "<", $userFile or die $!;

# read the user file
my %user;
while (<USER>) {
    chomp;                  # no newline
    s/#.*//;                # no comments
    s/^\s+//;               # no leading white
    s/\s+$//;               # no trailing white
    next unless length;     # anything left?
    my ($var, $value) = split(/\s*=\s*/, $_, 2); # separate with a equal sign
    $user{$var} = $value;
} 












if($aligner eq 'bowtie'){
  # print a warning about the sam files
  print "There seems to be something wrong with the SAM files produced in this program.\nSamtools is used to convert the bowtie alignment, but there is something wrong with the result.\nThe bowtie alignment itself (.map) seems to be fine though.\n";
  
  # get option
  my $index = $opt_i or die "ERROR: Full path and prefix of Bowtie index (-i) not specified\n";
  $index =~ /(.+)\/(\S+)$/;  # extract the last non-whitespace after the last /  ie. the bare index prefix
  my $indexPath = $1;
  $index = $2;
  
  # write submission file
  print SBATCH <<EOF;
#!/bin/bash -l
#SBATCH -A $user{'proj'}
#SBATCH -p $user{'type'}
#SBATCH -o $out/$user{'out'}
#SBATCH -e $out/$user{'err'}
#SBATCH --mail-type=$user{'emailType'}
#SBATCH --mail-user=$user{'emailAdd'}
#SBATCH -t $user{'reqTime'}
#SBATCH -J $aligner.$stamp.$user{'jobName'}

module load bioinfo-tools
module load bowtie

export BOWTIE_INDEXES=$indexPath

# align the reads to the index using $cores threads
bowtie -C --sam -p $cores $index -f $reads -Q $qual $out/$fileName.aligned.sam

# convert the sam file to bam
samtools view -bS -o $out/$fileName.aligned.bam $out/$fileName.aligned.sam
EOF

  # submit and remove the file
  #system("sbatch $out/$stamp.sbatch");
  #system("rm $out/$stamp.sbatch");
  print "To start the alignment, please submit your file to the queue system by typing:\nsbatch $out/$stamp.sbatch\n";



















}elsif($aligner eq 'bwa'){
  
  # get the index
  my $index = $opt_i or die "ERROR: Reference genome fasta file (-i) not specified\n";
  $index = absPath($index);
  $index =~ /.+\/(\S+)$/;  # extract the last non-whitespace after the last /  ie. the bare index prefix and the path
  my $indexName = $1;
  
  # pick out the read and qual prefix (retarded program..)
  my $prefix;
  my $readName;
  if($fileName =~ m/(.+_)([FR]3).csfasta/){
    $prefix = $1;
    $readName = "$prefix$2";
  }else{
    die "The authors of maq thought it would be neat if the reads and quality file had to be named whatevernameyouwant_F3.csfasta and whatevernameyouwant_F3_QV.qual\nSo if you get this message, it means your files where not named it this manner. So rename them and try again.\n";
  }
  
  # write submission file
  print SBATCH <<EOF;
#!/bin/bash -l
#SBATCH -A $user{'proj'}
#SBATCH -p $user{'type'}
#SBATCH -o $out/$user{'out'}
#SBATCH -e $out/$user{'err'}
#SBATCH --mail-type=$user{'emailType'}
#SBATCH --mail-user=$user{'emailAdd'}
#SBATCH -t $user{'reqTime'}
#SBATCH -J $aligner.$stamp.$user{'jobName'}

module load bioinfo-tools
module load bwa
module load samtools


# convert the solid reads to fastq format, using a short prefix (reads) to save disc space
solid2fastq.pl $readPath/$prefix $out/reads

# unzip the weird zip file (huge!)  According to Jonatan this is not needed. No time to test it, so i will leave it in..
gunzip $out/reads.single.fastq.gz

# align the reads to the index file, using $cores threads
bwa aln -t $cores -c $index $out/reads.single.fastq > $out/$fileName.sai

# convert the result to sam format
bwa samse $index $out/$fileName.sai $out/reads.single.fastq > $out/$fileName.aligned.sam

# convert the sam file to bam
samtools view -bS -o $out/$fileName.aligned.bam $out/$fileName.aligned.sam

EOF
  
  # submit and remove the file
  #system("sbatch $out/$stamp.sbatch");
  #system("rm $out/$stamp.sbatch");
  print "To start the alignment, please submit your file to the queue system by typing:\nsbatch $out/$stamp.sbatch\n";
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
# Bioscope single end
}elsif($aligner eq 'bioscopeSE'){
  
  #### SETTINGS TO BE USED (more detailed settings can be edited further down in this file)
  my $runName = "Run1";
  my $sampleName = "Sample1";
  my $type = "F3";
  my $index = $opt_i or die "ERROR: Path to reference genome fasta file (-i) not specified\n";
  $index = absPath($index);
  my $mappingQualFilterCutoff = 0;
  
  ### Generate the folder tree, if it does not exist
  mkpath("$out/config");
  mkpath("$out/intermediate");
  mkpath("$out/log");
  mkpath("$out/output");
  mkpath("$out/temp");
  mkpath("$out/tmp");
  
  # take the first 30 chars of the file name and use as a shorter namn (is there no easier way?)
  my $shortName;
  if(length($fileName) > 30){
    $shortName = substr $fileName,0,30;
  }else{
    $shortName = substr $fileName,0,length($fileName);
  }
  
  
  ### Create config files
  
  # open file handles
  open PLAN, ">", "$out/config/analysis.plan" or die $!;
  open MAP, ">", "$out/config/Mapping.ini" or die $!;
  
  
  # write analysis.plan
  print PLAN <<EOF;

$out/config/Mapping.ini

EOF


  # write Mapping.ini
  print MAP <<EOF;
pipeline.cleanup.middle.files=1
pipeline.cleanup.temp.files=1


#    Global Parameters 

run.name=$runName
sample.name=$sampleName
reference=$index
output.dir=$out/output
tmp.dir=$out/tmp
intermediate.dir=$out/intermediate
log.dir=$out/log
scratch.dir=/scratch/solid

######################################
##  mapping.run
######################################
mapping.run=1
mapping.output.dir=$out/output/$type/s_mapping
mapping.stats.output.file=$out/output/$type/s_mapping/mapping-stats.txt
read.length=50
mapping.tagfiles=$reads
primer.set=$type
mapping.run.classic=false
matching.max.hits=10
mapping.mismatch.penalty=-2.0
mapping.qual.filter.cutoff=$mappingQualFilterCutoff
clear.zone=5
mismatch.level=4

######################################
##  ma.to.bam.run
######################################
ma.to.bam.run=1
mapping.output.dir=$out/output/$type/s_mapping
ma.to.bam.qual.file=$qual
ma.to.bam.reference=$index
ma.to.bam.output.dir=$out/output/maToBam
ma.to.bam.match.file=$out/output/$type/s_mapping/$fileName.ma
ma.to.bam.output.file.name=$out/output/maToBam/$fileName.ma.bam
ma.to.bam.intermediate.dir=$out/intermediate/maToBam
ma.to.bam.temp.dir=$out/output/temp/maToBam
ma.to.bam.pas.file=
ma.to.bam.output.filter=primary
ma.to.bam.correct.to=reference
ma.to.bam.clear.zone=5
ma.to.bam.mismatch.penalty=-2.0
ma.to.bam.library.type=fragment
ma.to.bam.library.name=lib1
ma.to.bam.slide.name=
ma.to.bam.description=
ma.to.bam.sequencing.center=freetext
ma.to.bam.tints=agy
ma.to.bam.base.qv.max=40
ma.to.bam.temporary.dir=

######################################
##  bioscope-reports.run
######################################
bioscope-reports.run=1
ma.to.bam.output.dir=$out/output/maToBam
reports.output.dir=$out/output/reports
coverage.output.file.1=$out/output/reports/coverage
histogram.output.file.1=$out/output/reports/histogram
annotation.output.file.1=$out/output/reports/annotation
report.input.bam.file=$out/output/maToBam/$fileName.ma.bam
histogram.file.name=coverage-histogram.txt

######################################
##  position.errors.run
######################################
position.errors.run=1
ma.to.bam.output.dir=$out/output/maToBam
position.errors.output.dir=$out/output/position-errors
position.errors.output.file=positionErrors.txt
show.first.position.as.base=1

######################################
##  dibayes.run
######################################
dibayes.run=0
dibayes.output.prefix=diBayes
reference=$index
dibayes.output.dir=$out/output
dibayes.working.dir=$out/tmp/dibayes
dibayes.log.dir=$out/log/dibayes
maximal.read.length=50
call.stringency=high
het.skip.high.coverage=0
input.file.info=$out/output/maToBam/$fileName.ma.bam:0:$out/output/position-errors/$fileName.ma_F3_positionErrors.txt:
detect.2.adjacent.snps=0
reads.no.indel=1
reads.only.unique=0
reads.include.no.mate=0
write.fasta=1
write.consensus=1
compress.consensus=0
input.file.type=BAM
cleanup.tmp.file=1


EOF

  # execute the run
    print SBATCH <<EOF;
#!/bin/bash -l
source /etc/profile
module load bioinfo-tools
module load bioscope

# run bioscope in the background, and let it work its magic
nohup sh bioscope.sh -A $user{'proj'} -l $out/log/log.log $out/config/analysis.plan
EOF

  #system("sh $out/$stamp.sbatch &");
  #system("rm $out/$stamp.sbatch");
  print "To start the alignment, please submit your file to the queue system by typing:\nsh $out/$stamp.sbatch &\n";



















# Bioscope mate pair run
}elsif($aligner eq 'bioscopeMP'){
  
  #### SETTINGS ####  some are set here
  my $mappingQualFilterCutoff =25;
  

  # option stuff
  my @mpReads = split(/,/,$reads);
  my @mpQuals = split(/,/,$qual);
  
  # make sure there is 2 of each
  $#mpReads == 1 ? my $doNothing1=1 : die "2 and only 2 sets of reads/quals should be used with bioscopeMP. Comma separated (-r path/read1,/path/read2 -q qual1,qual2)\n";
  $#mpQuals == 1 ? my $doNothing2=1 : die "2 and only 2 sets of reads/quals should be used with bioscopeMP. Comma separated (-r path/read1,/path/read2 -q qual1,qual2)\n";
  
  # put them in a hash to define which is F3 and R3
  my %reads;
  my %quals;
  my %fileNames;
  my %readPaths;
  for(my $i = 0; $i <= $#mpReads; $i++){
    
    if($mpReads[$i] =~ m/^.*F3.*$/){ # F3
      $reads{'F3'} = absPath($mpReads[$i]);
      $quals{'F3'} = absPath($mpQuals[$i]);
      
      # extract the last non-whitespace after the last /  ie. the bare index prefix and the path
      $reads{'F3'} =~ /(.+)\/(\S+)$/;
      $readPaths{'F3'} = $1;
      $fileNames{'F3'} = $2;
      
    }elsif($mpReads[$i] =~ m/^.*R3.*$/){ # R3
      $reads{'R3'} = absPath($mpReads[$i]);
      $quals{'R3'} = absPath($mpQuals[$i]);
      
      # extract the last non-whitespace after the last /  ie. the bare index prefix and the path
      $reads{'R3'} =~ /(.+)\/(\S+)$/;
      $readPaths{'R3'} = $1;
      $fileNames{'R3'} = $2;
      
    }else{ # neither, something is wrong!!!11
      die "ERROR: at least one of the read files name does NOT include 'F3' or 'R3'. Unable to determin which is which!\n";
    }
  }

  
  # get index
  my $index = $opt_i or die "ERROR: Path to reference genome fasta file (-i) not specified\n";
  $index = absPath($index);
  
  ######################################################################
  #Should check for file permissions to write reference properties file#
  #                        todo maybe? :)                              #
  ######################################################################
  
  
  ### Generate the folder tree, if it does not exist
  mkpath("$out/config");
  
  # take the first 30 chars of the file name and use as a shorter namn (is there no easier way?)
  my $shortName;
  if(length($fileName) > 30){
    $shortName = substr $fileName,0,30;
  }else{
    $shortName = substr $fileName,0,length($fileName);
  }
  
  
  ### Create config files
  
  # open file handles
  open PLAN, ">", "$out/config/analysis.plan" or die $!;
  open F3, ">", "$out/config/F3.ini" or die $!;
  open R3, ">", "$out/config/R3.ini" or die $!;
  open MP, ">", "$out/config/matePair.ini" or die $!;
  
  
  # write analysis.plan
  print PLAN <<EOF;
$out/config/F3.ini
$out/config/R3.ini

$out/config/matePair.ini

EOF


  # write F3.ini
  print F3 <<EOF;

pipeline.cleanup.middle.files=1
pipeline.cleanup.temp.files=1


#    Global Parameters 

run.name=F3run
sample.name=$fileNames{'F3'}
reference=$index
output.dir=$out/output
tmp.dir=$out/tmp
intermediate.dir=$out/intermediate
log.dir=$out/log
scratch.dir=/scratch/solid
primer.set=F3

######################################
##  mapping.run
######################################

mapping.run=1
mapping.output.dir=$out/output/F3/s_mapping
mapping.stats.output.file=$out/output/F3/s_mapping/mapping-stats.txt
mapping.tagfiles=$reads{'F3'}
read.length=50
mapping.run.classic=false
matching.max.hits=10
mapping.mismatch.penalty=-2.0
mapping.qual.filter.cutoff=$mappingQualFilterCutoff
clear.zone=5
mismatch.level=4
mapping.tagfiles.f3=$reads{'F3'}


######################################
##  ma.to.bam.run
######################################
ma.to.bam.run=1
mapping.output.dir=$out/output/F3/s_mapping
ma.to.bam.qual.file=$quals{'F3'}
ma.to.bam.reference=$index
ma.to.bam.output.dir=$out/output/maToBam
ma.to.bam.match.file=$out/output/F3/s_mapping/$fileNames{'F3'}.ma
ma.to.bam.output.file.name=$out/output/maToBam/$fileNames{'F3'}.ma.bam
ma.to.bam.intermediate.dir=$out/intermediate/maToBam
ma.to.bam.temp.dir=$out/output/temp/maToBam
ma.to.bam.pas.file=
ma.to.bam.output.filter=primary
ma.to.bam.correct.to=reference
ma.to.bam.clear.zone=5
ma.to.bam.mismatch.penalty=-2.0
ma.to.bam.library.type=fragment
ma.to.bam.library.name=lib1
ma.to.bam.slide.name=
ma.to.bam.description=
ma.to.bam.sequencing.center=freetext
ma.to.bam.tints=agy
ma.to.bam.base.qv.max=40
ma.to.bam.temporary.dir=

######################################
##  bioscope-reports.run
######################################
bioscope-reports.run=1
ma.to.bam.output.dir=$out/output/maToBam
reports.output.dir=$out/output/reports
coverage.output.file.1=$out/output/reports/coverage
histogram.output.file.1=$out/output/reports/histogram
annotation.output.file.1=$out/output/reports/annotation
report.input.bam.file=$out/output/maToBam/$fileNames{'F3'}.ma.bam
histogram.file.name=coverage-histogram.txt

######################################
##  position.errors.run
######################################
position.errors.run=1
ma.to.bam.output.dir=$out/output/maToBam
position.errors.output.dir=$out/output/position-errors
position.errors.output.file=positionErrors.txt
show.first.position.as.base=1

######################################
##  dibayes.run
######################################
dibayes.run=0
dibayes.output.prefix=diBayes
reference=$index
dibayes.output.dir=$out/output/dibayes
dibayes.working.dir=$out/tmp/dibayes
dibayes.log.dir=$out/log/dibayes
maximal.read.length=50
call.stringency=high
het.skip.high.coverage=0
input.file.info=$out/output/maToBam/$fileNames{'F3'}.ma.bam:0:$out/output/position-errors/$fileNames{'F3'}.ma_F3_positionErrors.txt:
detect.2.adjacent.snps=0
reads.no.indel=1
reads.only.unique=0
reads.include.no.mate=0
write.fasta=1
write.consensus=1
compress.consensus=0
input.file.type=BAM
cleanup.tmp.file=1


EOF



  # write R3.ini
  print R3 <<EOF;

pipeline.cleanup.middle.files=1
pipeline.cleanup.temp.files=1


#    Global Parameters 

run.name=R3run
sample.name=$fileNames{'R3'}
reference=$index
output.dir=$out/output
tmp.dir=$out/tmp
intermediate.dir=$out/intermediate
log.dir=$out/log
scratch.dir=/scratch/solid
primer.set=R3

######################################
##  mapping.run
######################################

mapping.run=1
mapping.output.dir=$out/output/R3/s_mapping
mapping.stats.output.file=$out/output/R3/s_mapping/mapping-stats.txt
mapping.tagfiles=$reads{'R3'}
read.length=50
mapping.run.classic=false
matching.max.hits=10
mapping.mismatch.penalty=-2.0
mapping.qual.filter.cutoff=$mappingQualFilterCutoff
clear.zone=5
mismatch.level=4
mapping.tagfiles.r3=$reads{'R3'}


######################################
##  ma.to.bam.run
######################################
ma.to.bam.run=1
mapping.output.dir=$out/output/R3/s_mapping
ma.to.bam.qual.file=$quals{'R3'}
ma.to.bam.reference=$index
ma.to.bam.output.dir=$out/output/maToBam
ma.to.bam.match.file=$out/output/R3/s_mapping/$fileNames{'R3'}.ma
ma.to.bam.output.file.name=$out/output/maToBam/$fileNames{'R3'}.ma.bam
ma.to.bam.intermediate.dir=$out/intermediate/maToBam
ma.to.bam.temp.dir=$out/output/temp/maToBam
ma.to.bam.pas.file=
ma.to.bam.output.filter=primary
ma.to.bam.correct.to=reference
ma.to.bam.clear.zone=5
ma.to.bam.mismatch.penalty=-2.0
ma.to.bam.library.type=fragment
ma.to.bam.library.name=lib1
ma.to.bam.slide.name=
ma.to.bam.description=
ma.to.bam.sequencing.center=freetext
ma.to.bam.tints=agy
ma.to.bam.base.qv.max=40
ma.to.bam.temporary.dir=

######################################
##  bioscope-reports.run
######################################
bioscope-reports.run=1
ma.to.bam.output.dir=$out/output/maToBam
reports.output.dir=$out/output/reports
coverage.output.file.1=$out/output/reports/coverage
histogram.output.file.1=$out/output/reports/histogram
annotation.output.file.1=$out/output/reports/annotation
report.input.bam.file=$out/output/maToBam/$fileNames{'R3'}.ma.bam
histogram.file.name=coverage-histogram.txt

######################################
##  position.errors.run
######################################
position.errors.run=1
ma.to.bam.output.dir=$out/output/maToBam
position.errors.output.dir=$out/output/position-errors
position.errors.output.file=positionErrors.txt
show.first.position.as.base=1

######################################
##  dibayes.run
######################################
dibayes.run=0
dibayes.output.prefix=diBayes
reference=$index
dibayes.output.dir=$out/output/dibayes
dibayes.working.dir=$out/tmp/dibayes
dibayes.log.dir=$out/log/dibayes
maximal.read.length=50
call.stringency=high
het.skip.high.coverage=0
input.file.info=$out/output/maToBam/$fileNames{'R3'}.ma.bam:0:$out/output/position-errors/$fileNames{'R3'}.ma_R3_positionErrors.txt:
detect.2.adjacent.snps=0
reads.no.indel=1
reads.only.unique=0
reads.include.no.mate=0
write.fasta=1
write.consensus=1
compress.consensus=0
input.file.type=BAM
cleanup.tmp.file=1

EOF



  # write MP.ini
  print MP <<EOF;

pipeline.cleanup.middle.files=1
pipeline.cleanup.temp.files=1


#    Global Parameters 

run.name=MPrun
sample.name=$fileNames{'F3'}
reference=$index
output.dir=$out/output
tmp.dir=$out/tmp
intermediate.dir=$out/intermediate
log.dir=$out/log
scratch.dir=/scratch/solid

######################################
##  pairing.run
######################################
pairing.run=1
mate.pairs.tagfile.dirs=$out/output/F3/s_mapping,$out/output/R3/s_mapping
pairing.output.dir=$out/output/pairing
pairing.color.qual.file.path.1=$quals{'F3'}
pairing.color.qual.file.path.2=$quals{'R3'}
mapping.output.dir=
bam.file.name=
unmapped.bam.file.name=
indel.evidence.list=indel-evidence-list.pas
indel.preset.parameters=1,3,4,5
max.base.qv=40
mate.pairs.rescue.level=0
mates.stats.report.name=pairingStats.stats
indel.max.hits=10
matching.max.hits=100
mapping.mismatch.penalty=-2.0
pairing.anchor.length=25
indel.min.non-matched.length=10
indel.max.mismatches=5
use.template.rescue.file=1
pairing.indel.max.mismatch.tag1=5
pairing.indel.max.mismatch.tag2=5
pair.uniqueness.threshold=10.0
max.insert.estimate=20000
min.insert.estimate=0
primer.set=F3,R3
pairing.mark.duplicates=1
pairing.correct.to=reference
pairing.tints=agy
pairing.library.name=
pairing.output.filter=primary

######################################
##  bioscope-reports.run
######################################
bioscope-reports.run=0
pairing.output.dir=$out/output/pairing
reports.output.dir=$out/output/reports
coverage.output.file.1=
coverage.output.file.2=
histogram.output.file.1=
histogram.output.file.2=
annotation.output.file.1=
annotation.output.file.2=
report.input.bam.file=
histogram.file.name=coverage-histogram.txt

######################################
##  position.errors.run
######################################
position.errors.run=0
pairing.output.dir=$out/output/pairing
position.errors.output.dir=$out/output/position-errors
position.errors.output.file=positionErrors.txt
show.first.position.as.base=1

EOF




  # execute the run
    print SBATCH <<EOF;
#!/bin/bash -l
source /etc/profile
module load bioinfo-tools
module load bioscope

# run bioscope and let it work its magic
nohup sh bioscope.sh -A $user{'proj'} -l $out/log/log.log $out/config/analysis.plan &
EOF

  #system("sh $out/$stamp.sbatch");
  #system("rm $out/$stamp.sbatch");
  print "To start the alignment, please submit your file to the queue system by typing:\nsh $out/$stamp.sbatch\n";















# Bioscope whole transcriptome
}elsif($aligner eq 'bioscopeWTSE'){
  die "NOT FUNCTIONAL YET!\n";
  #### SETTINGS TO BE USED (more detailed settings can be edited further down in this file)
  my $runName = "Run1";
  my $sampleName = "Sample1";
  my $type = "F3";
  my $index = $opt_i or die "ERROR: Path to reference genome fasta file (-i) not specified\n";
  $index = absPath($index);
  my $mappingQualFilterCutoff = 0;
  my $extraPath = "/bubo/home/h18/dahlo/glob/work/projects/10ad_wholeTranscriptiome/files/cll_refs"; # path to the gtf and filter files
  
  ### Generate the folder tree, if it does not exist
  mkpath("$out/config");
  mkpath("$out/intermediate");
  mkpath("$out/log");
  mkpath("$out/output");
  mkpath("$out/temp");
  mkpath("$out/tmp");
  
  # take the first 30 chars of the file name and use as a shorter namn (is there no easier way?)
  my $shortName;
  if(length($fileName) > 30){
    $shortName = substr $fileName,0,30;
  }else{
    $shortName = substr $fileName,0,length($fileName);
  }
  
  
  ### Create config files
  
  # open file handles
  open PLAN, ">", "$out/config/analysis.plan" or die $!;
  open MAP, ">", "$out/config/WT_single_read.ini" or die $!;
  
  
  # write analysis.plan
  print PLAN <<EOF;

$out/config/WT_single_read.ini

EOF


  # write WT_single_read.ini
  print MAP <<EOF;
pipeline.cleanup.middle.files=1
pipeline.cleanup.temp.files=1

#    Global Parameters 

run.name=$runName
sample.name=$sampleName
reference=$index
output.dir=$out/output
tmp.dir=$out/tmp
intermediate.dir=$out/intermediate
log.dir=$out/log
scratch.dir=/scratch/solid

######################################
##  wt.genomic.mapping.run
######################################
wt.genomic.mapping.run=1
wt.merge.genomic.map.result.dir=$out/intermediate/s_mapping/genomic_map
wt.genomic.mapping.reference=$index
matching.max.hits=10
read.length=50
mapping.tagfiles=$reads
wt.tmp.dir=$out/tmp

######################################
##  wt.filter.mapping.run
######################################
wt.filter.mapping.run=1
wt.merge.filtered.map.result.dir=$out/intermediate/s_mapping/filter_map
wt.filter.mapping.reference=$index
matching.max.hits=10
mapping.tagfiles=$reads
wt.tmp.dir=$out/tmp

######################################
##  wt.spljunctionextractor.run
######################################
wt.spljunctionextractor.run=1
wt.splext.output.dir=$out/intermediate/spljunctionextraction
wt.splext.reference.file=$index
wt.splext.genegtf.file=$extraPath/refseqHumanFeb2009.gtf

######################################
##  wt.junction.mapping.run
######################################
wt.junction.mapping.run=1
wt.splext.output.dir=$out/intermediate/spljunctionextraction
wt.merge.junction.map.result.dir=$out/intermediate/s_mapping/junction_map
matching.max.hits=10
mapping.tagfiles=$reads
wt.tmp.dir=$out/tmp

######################################
##  wt.merge.run
######################################
wt.merge.run=1
wt.merge.genomic.map.result.dir=$out/intermediate/s_mapping/genomic_map
wt.merge.output.dir=$out/output/mapping
wt.merge.output.bam.file=wt.sr.bam
wt.merge.reference.file=$index
wt.merge.filter.reference.file=$extraPath/human_filter_reference.fasta
wt.merge.qual.file=$qual
wt.merge.tmpdir=$out/tmp/merge
wt.merge.filtered.map.result.dir=$out/intermediate/s_mapping/filter_map
wt.merge.junction.map.result.dir=$out/intermediate/s_mapping/junction_map
wt.merge.junction.reference.file=$out/intermediate/spljunctionextraction/junction.fasta
wt.merge.score.clear.zone=5
wt.merge.min.junction.overhang=8
wt.merge.num.alignments.to.store=1

######################################
##  wt.counttag.run
######################################
wt.counttag.run=1
wt.counttag.input.bam.file=$out/output/mapping/wt.sr.bam
wt.counttag.output.dir=$out/output/counttag
wt.counttag.output.file.name=countagresult.txt
wt.counttag.exon.reference=$extraPath/refseqHumanFeb2009.gtf
read.length=50
wt.counttag.score.clear.zone=5
wt.counttag.alignment.filter.mode=primary
wt.counttag.min.alignment.score=0
wt.counttag.min.mapq=10

EOF


  # write Mapping.ini
  print MAP <<EOF;

EOF

  # execute the run
    print SBATCH <<EOF;
#!/bin/bash -l
source /etc/profile
module load bioinfo-tools
module load bioscope

# run bioscope in the background, and let it work its magic
nohup sh bioscope.sh -A $user{'proj'} -l $out/log/log.log $out/config/analysis.plan
EOF

  #system("sh $out/$stamp.sbatch &");
  #system("rm $out/$stamp.sbatch");

















# Bioscope paired end run
}elsif($aligner eq 'bioscopePE'){
  
  #### SETTINGS ####  some are set here
  my $mappingQualFilterCutoff =0;
  

  # option stuff
  my @mpReads = split(/,/,$reads);
  my @mpQuals = split(/,/,$qual);
  
  # make sure there is 2 of each
  $#mpReads == 1 ? my $doNothing1=1 : die "2 and only 2 sets of reads/quals should be used with bioscopePE. Comma separated (-r path/read1,/path/read2 -q qual1,qual2)\n";
  $#mpQuals == 1 ? my $doNothing2=1 : die "2 and only 2 sets of reads/quals should be used with bioscopePE. Comma separated (-r path/read1,/path/read2 -q qual1,qual2)\n";
  
  # put them in a hash to define which is F3 and F5
  my %reads;
  my %quals;
  my %fileNames;
  my %readPaths;
  for(my $i = 0; $i <= $#mpReads; $i++){
    
    if($mpReads[$i] =~ m/^.*F3.*$/){ # F3
      $reads{'F3'} = absPath($mpReads[$i]);
      $quals{'F3'} = absPath($mpQuals[$i]);
      
      # extract the last non-whitespace after the last /  ie. the bare index prefix and the path
      $reads{'F3'} =~ /(.+)\/(\S+)$/;
      $readPaths{'F3'} = $1;
      $fileNames{'F3'} = $2;
      
    }elsif($mpReads[$i] =~ m/^.*F5.*$/){ # F5
      $reads{'F5'} = absPath($mpReads[$i]);
      $quals{'F5'} = absPath($mpQuals[$i]);
      
      # extract the last non-whitespace after the last /  ie. the bare index prefix and the path
      $reads{'F5'} =~ /(.+)\/(\S+)$/;
      $readPaths{'F5'} = $1;
      $fileNames{'F5'} = $2;
      
    }else{ # neither, something is wrong!!!11
      die "ERROR: at least one of the read files name does NOT include 'F3' or 'F5'. Unable to determin which is which!\n";
    }
  }

  
  # get index
  my $index = $opt_i or die "ERROR: Path to reference genome fasta file (-i) not specified\n";
  $index = absPath($index);
  
  ######################################################################
  #Should check for file permissions to write reference properties file#
  #                        todo maybe? :)                              #
  ######################################################################
  
  
  ### Generate the folder tree, if it does not exist
  mkpath("$out/config");
  
  # take the first 30 chars of the file name and use as a shorter namn (is there no easier way?)
  my $shortName;
  if(length($fileName) > 30){
    $shortName = substr $fileName,0,30;
  }else{
    $shortName = substr $fileName,0,length($fileName);
  }
  
  
  ### Create config files
  
  # open file handles
  open PLAN, ">", "$out/config/analysis.plan" or die $!;
  open F3, ">", "$out/config/F3.ini" or die $!;
  open F5, ">", "$out/config/F5.ini" or die $!;
  open PE, ">", "$out/config/pairedEnd.ini" or die $!;
  
  
  # write analysis.plan
  print PLAN <<EOF;
=$out/config/F3.ini
=$out/config/F5.ini

$out/config/pairedEnd.ini

EOF


  # write F3.ini
  print F3 <<EOF;

pipeline.cleanup.middle.files=1
pipeline.cleanup.temp.files=1


#    Global Parameters 

run.name=F3run
sample.name=$fileNames{'F3'}
reference=$index
output.dir=$out/output
tmp.dir=$out/tmp
intermediate.dir=$out/intermediate
log.dir=$out/log
scratch.dir=/scratch/solid

######################################
##  mapping.run
######################################
mapping.run=1
mapping.output.dir=$out/output/F3/s_mapping
mapping.stats.output.file=$out/output/F3/s_mapping//mapping-stats.txt
mapping.tagfiles=$reads{'F3'}
read.length=50
mapping.run.classic=false
matching.max.hits=100
mapping.mismatch.penalty=-2.0
mapping.qual.filter.cutoff=$mappingQualFilterCutoff
clear.zone=5
mismatch.level=5
mapping.tagfiles.f3=$reads{'F3'}

EOF



  # write F5.ini
  print F5 <<EOF;

pipeline.cleanup.middle.files=1
pipeline.cleanup.temp.files=1


#    Global Parameters 

run.name=F5run
sample.name=$fileNames{'F5'}
reference=$index
output.dir=/$out/output
tmp.dir=/$out/tmp
intermediate.dir=/$out/intermediate
log.dir=/$out/log
scratch.dir=/scratch/solid

######################################
##  mapping.run
######################################
mapping.run=1
mapping.output.dir=/$out/output/F5/s_mapping
mapping.stats.output.file=/$out/output/F5/s_mapping//mapping-stats.txt
mapping.tagfiles=$reads{'F5'}
read.length=35
mapping.run.classic=false
matching.max.hits=100
mapping.mismatch.penalty=-2.0
mapping.qual.filter.cutoff=$mappingQualFilterCutoff
clear.zone=5
mismatch.level=5
mapping.tagfiles.f5=$reads{'F5'}

EOF



  # write PE.ini
  print PE <<EOF;

pipeline.cleanup.middle.files=1
pipeline.cleanup.temp.files=1


#    Global Parameters 

run.name=PErun
sample.name=$fileNames{'F3'}
reference=$index
output.dir=/$out/output
tmp.dir=/$out/tmp
intermediate.dir=/$out/intermediate
log.dir=/$out/log
scratch.dir=/scratch/solid

######################################
##  paired-end-pairing.run
######################################
paired-end-pairing.run=1
mate.pairs.tagfile.dirs=/$out/output/F3/s_mapping,/$out/output/F5/s_mapping
pairing.output.dir=/$out/output/pairing
primer.set=F3,F5-BC
pairing.color.qual.file.path.1=$quals{'F3'}
pairing.color.qual.file.path.2=$quals{'F5'}
mapping.output.dir=
bam.file.name=
unmapped.bam.file.name=
indel.evidence.list=indel-evidence-list.pas
indel.preset.parameters=0
max.base.qv=40
mate.pairs.rescue.level=4
mates.stats.report.name=pairingStats.stats
indel.max.hits=10
matching.max.hits=100
mapping.mismatch.penalty=-2.0
pairing.anchor.length=25
indel.min.non-matched.length=10
indel.max.mismatches=5
use.template.rescue.file=1
pairing.indel.max.mismatch.tag1=5
pairing.indel.max.mismatch.tag2=2
pair.uniqueness.threshold=10.0
max.insert.estimate=20000
min.insert.estimate=0
pairing.mark.duplicates=0
pairing.correct.to=reference
pairing.tints=agy
pairing.library.name=
pairing.output.filter=primary

EOF




  # execute the run
    print SBATCH <<EOF;
#!/bin/bash -l
source /etc/profile
module load bioinfo-tools
module load bioscope

# run bioscope and let it work its magic
nohup sh bioscope.sh -A $user{'proj'} -l $out/log/log.log $out/config/analysis.plan &
EOF

  #system("sh $out/$stamp.sbatch");
  #system("rm $out/$stamp.sbatch");
  print "To start the alignment, please submit your file to the queue system by typing:\nsh $out/$stamp.sbatch &\n";












# MAQ works like it should, it's just that it really eats memory.
# It will probably work with smaller genomes and data sets,
# but a full run with human genome will eat up 24Gb ram like nothing..
# Solution: split the read file into smaller pieces and use
# maq merge on the resulting files. I will focus on other things instead.
}elsif($aligner eq 'maq'){
  
  # print warning. Feel free to remove this if you have fixed the splitting and merging. Or bought a really cool computer :)
  print "Do not use for larget genomes and/or large data sets! Read warning in\npolyAlign.pl:313 for more info.\n(Should be ok to run smaller things)\n";
  
  my $index = $opt_i or die "ERROR: Reference genome fasta file (-i) not specified\n";
  $index = absPath($index);
  $index =~ /.+\/(\S+)$/;  # extract the last non-whitespace after the last /  ie. the bare index prefix and the path
  my $indexName = $1;
  
  # pick out the read and qual prefix (retarded program..)
  my $prefix;
  my $readName;
  if($fileName =~ m/(.+_)([FR]3).csfasta/){
    $prefix = $1;
    $readName = "$prefix$2";
  }else{
    die "The authors of maq thought it would be neat if the reads and quality file had to be named whatevernameyouwant_F3.csfasta and whatevernameyouwant_F3_QV.qual\nSo if you get this message, it means your files where not named it this manner. So rename them and try again.\n";
  }
  
  
  # write submission file
  print SBATCH <<EOF;
#!/bin/bash -l
#SBATCH -A $user{'proj'}
#SBATCH -p $user{'type'}
#SBATCH -o $out/$user{'out'}
#SBATCH -e $out/$user{'err'}
#SBATCH --mail-type=$user{'emailType'}
#SBATCH --mail-user=$user{'emailAdd'}
#SBATCH -t $user{'reqTime'}
#SBATCH -J $aligner.$stamp.$user{'jobName'}
#SBATCH -C fat

module load bioinfo-tools
module load maq
module load samtools
module load bwa

# convert the solid reads to fastq format
solid2fastq.pl $readPath/$prefix $out/$readName

# unzip the weird zip file (huge!)
gunzip $out/$readName.single.fastq.gz

# NOTE TO SELF (and others): -n $cores really screwed things up.. Avoid!

# convert the now-unzipped reads to binary fastq
maq fastq2bfq $out/$readName.single.fastq $out/$readName.bfq

# convert the basespace reference index to colorspace fasta
maq fasta2csfa $index > $out/$indexName.csfa

# convert the colorspace fasta reference index to colorspace BINARY fasta
maq fasta2bfa $out/$indexName.csfa $out/$indexName.csbfa

# convert the basespace reference index to binary fasta
maq fasta2bfa $index $out/$indexName.bfa

# map the reads using both basespace and colorspace
maq map -c $out/aln.cs.map $out/$indexName.csbfa $out/$readName.bfq 2> $out/aln.log

# convert the result to sam
maq2sam-long $out/aln.cs.map > $out/aln.sam

EOF
  
  # submit and remove the file
  #system("sbatch $out/$stamp.sbatch");
  #system("rm $out/$stamp.sbatch");
  print "To start the alignment, please submit your file to the queue system by typing:\nsbatch $out/$stamp.sbatch\n";
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  

  
}elsif($aligner eq 'mosaik'){

  # prepare variables
  my $jump = $opt_j or die "ERROR: Jump database path and prefix (-j) not specified\n";
  $jump = absPath($jump);
  my $index = $opt_i or die "ERROR: Full path and name of reference fasta file (-i) not specified\n";
  $index =~ /(.+)\/(\S+)$/;  # extract the last non-whitespace after the last /  ie. the bare file name
  my $indexPath = $1;
  $index = $2;
  
  # create the reference directory
  mkpath("$out/ref");
  
  # write submission file
  print SBATCH <<EOF;
#!/bin/bash -l

#SBATCH -A $user{'proj'}
#SBATCH -p $user{'type'}
#SBATCH -o $out/$user{'out'}
#SBATCH -e $out/$user{'err'}
#SBATCH --mail-type=$user{'emailType'}
#SBATCH --mail-user=$user{'emailAdd'}
#SBATCH -t $user{'reqTime'}
#SBATCH -J $aligner.$stamp.$user{'jobName'}
#SBATCH -C fat


module add bioinfo-tools
module add mosaik-aligner


# convert reference genome to mosaik format, basespace
MosaikBuild -fr $indexPath/$index -oa $out/ref/$index.fa.dat

# convert reference genome to mosaik format, colorspace
MosaikBuild -cs -fr $indexPath/$index -oa $out/ref/$index.fa.cs.dat

# convert the reads to a read archive mosaik format
MosaikBuild -fr $reads -fq $qual -out $out/$fileName.dat -st solid

# align the read archive to the reference genome
MosaikAligner -in $out/$fileName.dat -out $out/$fileName.aligned.dat -ia $out/ref/$index.fa.cs.dat -ibs $out/ref/$index.fa.dat -hs 15 -mm 4 -mhp 100 -act 20 -rur $out/$fileName.unmapped.fasta -p $cores -j $jump

# sort the result
MosaikSort -in $out/$fileName.aligned.dat -out $out/$fileName.sorted.dat -afl -consed 

# covert it to a bam file
MosaikText -in $out/$fileName.sorted.dat -bam $out/$fileName.sorted.bam

# do some.. coverage stuff
MosaikCoverage -in $out/$fileName.sorted.dat -ia $out/ref/$index.fa.cs.dat -u -od $out -cg
EOF

  # submit and remove the file
  #system("sbatch $out/$stamp.sbatch");
  #system("rm $out/$stamp.sbatch");
  print "To start the alignment, please submit your file to the queue system by typing:\nsbatch $out/$stamp.sbatch\n";
  
  
  
  
}else{
  die $validAligners;
}














### Methods

# Removes the / sign, if any, at the end of directory names
sub remSlash{
  my $str = shift;
  if($str =~ m/\/$/){ # if the last char is a /
    chop($str);
  }
  return $str;
}

# Adjust relative path to absolute paths
sub absPath{
  my $str = shift;
  if($str !~ m/^\//){ # if the first char is not a / ie. a relative path
    $str = &Cwd::cwd()."/$str"; # attach the current work directory infront of the relative path
  }
  return $str;  
}
