alignReads

#!/usr/bin/env perl
require '/home/apa/apa_routines.pm';
use Data::Dumper;
use Getopt::Long;
use Pod::Usage;
use strict;


## TO DO: a million things.
## INCLUDING: better handling for tabular align_summary.txt filenames, test ability to go without $anno, MAKE $PARAMS LIVE, 
## PRIORITY: ensure strand assignments are correct with HTseq-count, with BAMs coming from various aligners...


## Performs genome alignment; returns mapped/unmapped read sets, mapping, quantitation, qc, and tracks


## Dependencies
my $bin = '/home/apa/local/bin';
my $bam2fq = "$bin/bam2fq";
my $bam2track = "$bin/bam2track";
my $readCount = "$bin/readCount";
my $fastqUnaligned = "$bin/fastqUnaligned";
my $mergeAlignSummary = "$bin/mergeAlignSummary";  
my $collapseFastqPair = "$bin/collapseFastqPair";
my $uncollapseFastqPair = "$bin/uncollapseFastqPair";
my $readGroupByLength = "$bin/readGroupByLength";
my $bamQuantitate = "$bin/bamQuantitate";
my $shortstack = "$bin/scriptutils/ShortStack/ShortStack-smallRNApipeline";
my $tophat = 'tophat';
my $bowtie2 = 'bowtie2';
my $bowtie1 = 'bowtie1';
my $bowtie_build = 'bowtie-build';  # bowtie1 only
my $star = 'STAR';
my $bwa = 'bwa';
my $samtools = 'samtools';  # REQUIRES SAMTOOLS >= 1.3.1
my $java = 'java';
my $picard = '/n/apps/CentOS7/bin/picard.jar';
my $htseq_count = 'htseq-count';


## Inputs
my $runmode;     # build type to align to, must be one of: genome, transcriptome, funcRNA, repeats, or ribosomes (future: add miRNA?)
my $geno;        # genome build label
my $anno;        # annot build label (optional)
my $outdir;      # output path
my $clobber;     # output path, if clobbering
my $fastq1;      # end 1 fastq
my $fastq2;      # end 2 fastq, if any
my $aligner;     # which aligner to use   ## CURRENTLY: tophat, shortstack   ## FUTURE: bowtie1, bowtie2, star, bwa
my $repeatfa;    # If 'repeats', use $geno/repeats build.  Else, use this custom repeats fasta; the bowtie1 index files found by replacing .(fa|fasta) with .*.ebwt SHOULD exist (else will create).
my $local;       # if bowtie2, use --local
my $params;      # quoted string containing non-default args to pass to aligner  *** PARAMS MUST BE COMPLETE -- ONLY THESE PARAMS WILL BE USED ***
my $cuffparams;  # quoted string containing non-default args to pass to cufflinks, if running cufflinks  *** PARAMS MUST BE COMPLETE -- ONLY THESE PARAMS WILL BE USED ***
my $cores = 1;   # N cores
my $mem = '1G';  # memory for 'samtools sort'
my $M;           # alternate M for RPKM; will get N input reads if not specified
my $runfqr;      # 1|0 generate fastqReads file for inputs?
my $quant;       # 1|0 quantitate gene models?
my $txome;       # 1|0 align to transcriptome only?
my $unique;      # 1|0 produce a multiread-free bam version?
my $stranded;    # 1|0 are reads stranded?  If so, will also generate stranded counts and bigWigs
my $switch;      # 1|0 switch strands for quantitation (e.g. Truseq directional protocol)
my $smallRNA;    # 1|0 are reads from a smallRNA prep? (alters default alignment params)
my $runrqc;      # 1|0 run RChipQC.R on final bam file?
my $runrsm;      # 1|0 run Picard's CollectRnaSeqMetrics on bam? (HOW TO INDICATE 'PLUS'?)
my $runcuff;     # 1|0 run Cufflinks? (HOW TO PASS CUSTOM ARGUMENTS?)
my $bw;          # 1|0 generate bigWigs?
my $tdf;         # 1|0 generate TDFs? (bigWig for IGV)
my $rlsplit;     # 1|0 generate a read-length read-grouped bam and split into individual read-length bams? (in separate folder)
my $quant;       # 1|0 quantitate features after alignment?  (if $rlsplit, will quantitate these as well)

############## FIXME: should these be $params?
my $htseq_python = '2.7.5';  # python version where desired htseq version is installed
my $alnmax = 100;     # maximum allowed alignments before a read is ignored; default 100
my $ranmax = int($alnmax/10);  # ShortStack --ranmax argument; default int($alnmax/10)
$alnmax = 100 unless $alnmax;   # default maximum alignments per read (before read is ignored)
$ranmax = int($alnmax/10) unless $ranmax;  # ShortStack --ranmax param


## MUST MAKE PRIMARY BAM CREATION DEPENDENT ON ALIGNER/CUSTOM-ARG-GREP RESULTS


## Get, test options
my $cmdline = "$0 @ARGV";
$runmode = shift @ARGV;   # DO BEFORE GetOptions -- MUST BE FIRST ARGUMENT
GetOptions(
    "g=s"=>\$geno, "a=s"=>\$anno, "o=s"=>\$outdir, "oc=s"=>\$clobber, "c=i"=>\$cores, "m=s"=>\$mem, 
    "fq1=s"=>\$fastq1, "fq2=s"=>\$fastq2, "Mval=i"=>\$M, "aln=s"=>\$aligner, "params=s"=>\$params, 
    "transcriptome"=>\$txome, "quant"=>\$quant, "stranded"=>\$stranded, "switch"=>\$switch, 
    "unique"=>\$unique, "smallRNA"=>\$smallRNA, "bw"=>\$bw, "tdf"=>\$tdf,  "rl-split"=>\$rlsplit, 
    "fqr"=>\$runfqr, "rqc"=>\$runrqc, "rsm"=>\$runrsm, "cufflinks"=>\$runcuff, "cuff-params=s"=>\$cuffparams
);

my @OKmodes = qw/ genome transcriptome funcRNA repeats ribosomes /;
my %OKmodes = map {($_=>1)} @OKmodes;
my $OKmodes = join(', ', map {"'$_'"} @OKmodes);
die "Run mode value '$runmode' not recognized: must be one of $OKmodes!\n" unless $OKmodes{$runmode};

$aligner = lc($aligner);
my $STAR_len;
if ($aligner =~ /^star[-_](\d+)/) {
    $STAR_len = $1;
    $aligner = 'star';
}
my $Gdir = "/n/data1/genomes/indexes/$geno";
my $Gpref = "$Gdir/$geno";
my $Adir = "$Gdir/$anno";
my $Apref = "$Adir/$geno.$anno";
my $gtf = "$Apref.cuff.gtf";
my $gtfi = "$gtf.index/$geno.$anno.cuff";
my $Sdir = $anno ? "$Adir/STAR_${STAR_len}bp" : "$Gdir/STAR_${STAR_len}bp";
die "$0: Expected dataset path '$Gdir' does not exist!\n" if $geno && ! -d $Gdir;
die "$0: Expected dataset path '$Adir' does not exist!\n" if $anno && ! -d $Adir;
die "$0: Expected STAR index '$Sdir' does not exist!\n" if $aligner eq 'star' && ! -d $Sdir;

die "$0: fastq1 '$fastq1' does not exist!\n" unless -e $fastq1;  # must exist
die "$0: fastq2 '$fastq2' does not exist!\n" if $fastq2 && ! -e $fastq2;  # must exist if specified
die "$0: cannot make tracks when aligning to transcriptome only!\n" if $txome && ($bw || $tdf);

my ($ribopref, $reptpref);
if ($runmode eq 'ribosomes') {
    $ribopref = "/n/data1/genomes/indexes/$geno/$anno/ribosomes/$geno.$anno.ribosomes";
    die "$0: ribosome build for genome '$geno' / '$anno' is incomplete or missing!\n" unless -e "$ribopref.1.bt2";  # this must exist
} elsif ($runmode eq 'repeats') {
    if ($repeatfa) {
        die "$0: indicated reference fasta '$ref_fa' does not exist!\n" unless -e $repeatfa;
        ($reptpref = $repeatfa) =~ s/\.f[sta]+$//;
        unless (-e "$reptpref.1.ebwt") {
            ## Bowtie1 index missing; create it
            print "Bowtie1 index for repeat fasta '$repeatfa' appears to be missing; creating one now...\n";
            system "$bowtie_build $repeatfa $reptpref > $reptpref.bowtie-build.log";
            print "Bowtie1 indexing complete, see '$reptpref.bowtie-build.log'\n";
        }
    } else {
        my $reptdir = "/n/projects/apa/stuff/bowtie_building/builds/$geno/repeats";
        $reptpref = "$reptdir/$geno.repeats";
        $repeatfa = "$reptpref.fa";
    }
    die "$0: repeat build for genome '$geno' is incomplete or missing!\n" unless -e "$reptpref.1.ebwt";  # testing one of the known bowtie1-index files
} elsif ($runmode eq 'funcRNA') {
} else {
    if ($smallRNA) {
        my $msg = "$0: Use of --smallRNA supersedes";
        if ($stranded) {
            if ($switch) {
                print "$msg --stranded and --switch: ignoring these...\n";
            } else {
                print "$msg --stranded: ignoring it...\n";
            }
        } elsif ($switch) {
            print "$msg --switch: ignoring it...\n";
        }
        $switch = 0;
        $stranded = 1;
    }
}

## Globals
my %files;    # filename control
my $nreads;   # N input reads (reads not pairs)
my $noUna;    # specified aligner does not write unaligned-read files
my $allbam;   # final all-reads bam post-alignment
my $pribam;   # final primary-reads bam (may be same as $allbam)
my $nounq;    # do not run no-multireads quantitation set
my $LOG;
my $PE = $fastq2 ? ' --PE' : '';
my $aligner2 = $aligner;  # initially
$local = ' --local' if $local;


## NEW RUN ##

## Output location
my $label = "alignTo$Runmode";
chomp(my $pwd = `pwd`);
(my $Runmode = $runmode) =~ s/^(.)/\U$1/;
$outdir = "$pwd/$label" unless $outdir;
if ($clobber) {
    $outdir = $clobber;
    system "rm -rf $outdir";
} elsif (-d $outdir) {
    die "$0: output dir '$outdir' already exists!\n";
}
system "mkdir -p $outdir";
die "$0: failed to create output dir '$outdir'!\n" unless -d $outdir;


## Log file
my $outpref = "$outdir/$label";
my $logfile = "$outpref.log";
$LOG = &open2('W', $logfile, 'log file');
&logreport("Initialized: ".`date`, 1, $LOG);


## Quantitation call variable + shell scripts
my $bamquantcall;
my $qsha = "$outpref.quantitate_aln.sh";
my $qshu = "$outpref.quantitate_unq.sh";
system "rm -f $qsha $qshu";


## Alignments
my $alndir = "$outpref.align.$$.tmp";

if ($runmode eq 'ribosomes') {
    
    
    $nounq = 1;
    my $inputs = $PE ? "-1 $fastq1 -2 $fastq2" : "-U $fastq1";
    my $unaflag = $PE ? '' : "--un-gz $outpref.unaligned.fastq.gz";
    my $alnsum1 = "$outdir/bowtie2.align_summary.txt";
    my $alnsum = "$outpref.align_summary.txt";
    &execute("$bowtie2$local -p $cores -x $ribopref $unaflag $inputs -S $outpref.tmp.sam 2> $alnsum1", 1, $LOG, 1);
    &execute("$samtools view -h -bS -F 4 $outpref.tmp.sam > $outpref.tmp.bam", 1, $LOG, 2);
    &execute("$samtools sort -@ $cores -m $mem -o $outpref.bam $outpref.tmp.bam", 1, $LOG, 2);
    
    ## Post-alignment
    &logreport("Running post-alignment:", 1, $LOG);
    &execute("$samtools index $outpref.bam", 1, $LOG, 2);
    &execute("$samtools idxstats $outpref.bam > $outpref.idxstats.txt", 1, $LOG, 2);
    system "rm -f $outpref.align_summary.txt";
    &execute("$mergeAlignSummary -o $alnsum -a bowtie2 $alnsum1", 1, $LOG, 2);  # tabularize
    &execute("perl -i -pe 's!^$alnsum1!alignToRibosomes!' $alnsum", 1, $LOG, 2);
    
    ## Extract final unaligned read sets
    if ($PE) {
        &logreport("Extracting unaligned reads: ".`date`, 1, $LOG);
        &execute("$fastqUnaligned -fq1 $fastq1 -fq2 $fastq2 -p $outpref.unaligned --no-orphans $outpref.bam", 1, $LOG, 3);
    }
    
    ## Get N input reads; set up BQ call; x remove temp stuff
    chomp($nreads = `head -2 $alnsum | tail -1 | cut -f2`);
    $bamquantcall = "$bamQuantitate -b $ribopref.bed --class --transcriptome";
    $aligner2 = 'bowtie2';  # for BQ call finalizing
    &execute("rm -f $outdir/*.tmp.*", 1, $LOG, 2);
    
    
} elsif ($runmode eq 'repeats') {
    
    
    ## Collapse paired fastqs, if paired
    ## Shortstack will only align single-end data
    my $infq = "$outdir/tmp.fastq";
    $infq .= '.gz' if $fastq1 =~ /\.gz$/;
    if ($fastq2) {
        &execute("$collapseFastqPair $fastq1 $fastq2 $infq", 1, $LOG, 2);
    } else {
        &execute("ln -sf \$(readlink -f $fastq1) $infq", 1, $LOG, 2);
    }
    
    ## Alignment
    my $alndir = "$outdir/alignToRepeats.align.tmp";
    my $alnname = 'ShortStack';
    my $alnpref = "$alndir/$alnname";
    my $alnsum = "$outpref.align_summary.txt";
    &execute("rm -rf $alndir", 1, $LOG, 2);  # else ShortStack dies
    &execute("$ShortStack --align_only --nostitch --nohp --bowtie_cores $cores --sort_mem $mem --mismatches 2 --mmap u --bowtie_m $alnmax --ranmax $ranmax --readfile $infq --bowtie_un $outpref.unal.fastq --outprefix $alnname --outdir $alndir --genomefile $repeatfa", 1, $LOG, 3);
    &logreport("Aligned: ".`date`, 1, $LOG);
    
    ## Post-alignment
    &logreport("Running post-alignment:", 1, $LOG);
    &execute("samtools view -h -F 4 $alnpref.bam | samtools view -bS - > $outpref.bam", 1, $LOG, 3);
    &execute("samtools index $outpref.bam", 1, $LOG, 3);
    &execute("samtools idxstats $outpref.bam > $outpref.idxstats.txt", 1, $LOG, 2);
    if ($PE) {
        &execute("$uncollapseFastqPair $outpref.unal.fastq $outpref.unaligned", 1, $LOG, 2);  # outputs gzipped
    } else {
        &execute("gzip -c $outpref.unal.fastq > $outpref.unaligned.fastq.gz", 1, $LOG, 2);
    }
    &execute("$mergeAlignSummary -o $alnsum -a shortstack --nodirname $alnpref.Log.txt", 1, $LOG, 2);
    &execute("mv -f $alndir/* $outdir/", 1, $LOG, 2);
    &execute("rm -rf $alndir", 1, $LOG, 2);
    &execute("rm -f $infq", 1, $LOG, 2);
    &execute("rm -f $outpref.unal.fastq", 1, $LOG, 2);
    
    ## Get N input reads; set up BQ call; x remove temp stuff
    chomp($nreads = `head -2 $alnsum | tail -1 | cut -f2`);
    $bamquantcall = "$bamQuantitate --transcriptome";
    $aligner2 = 'shortstack';  # for BQ call finalizing
    &execute("rm -f $outdir/*.tmp.*", 1, $LOG, 2);
    
    
} elsif ($runmode eq 'funcRNA') {
    
    
    
} else {
    
    
    ## Transcriptome or Genome alignment
    
    $bamquantcall = "$bamQuantitate -g $geno -a $anno";
    
    
    if ($aligner eq 'shortstack') {
        
        
        my $SSname = 'shortstack';
        my $SSpref = "$outdir/$SSname";
        my $SSfqin = "$SSpref-input.fastq.gz";
        my $noUna = 0;
        if ($PE) {
            &execute("$collapseFastqPair $fastq1 $fastq2 $SSfqin", 1, $LOG, 2);
        } else {
            &execute("ln -sf \$(readlink -f $fastq1) $SSfqin", 1, $LOG, 2, 1);  ## REQUIRES BASH
        }
        my $ref_fa = $txome ? "$gtfi.fa" : "$Gpref.fa";
        unless ($params) {
            $params .= " --align_only --nostitch --nohp --mismatches 2 --mmap u --bowtie_m $alnmax --ranmax $ranmax";
        }
        &execute("$shortstack $params --genomefile $ref_fa --bowtie_cores $cores --sort_mem $mem --readfile $SSfqin --bowtie_un $SSpref.unal.fastq --outprefix $SSname --outdir $alndir", 1, $LOG, 2);
        if ($PE) {
            &execute("$uncollapseFastqPair $SSpref.unal.fastq $outpref.unaligned", 1, $LOG, 2);  # outputs gzipped
        } else {
            &execute("mv -f $SSpref.unal.fastq $outpref.unaligned.fastq", 1, $LOG, 2);
            &execute("gzip -f $outpref.unaligned.fastq", 1, $LOG, 2) if -e "$outpref.unaligned.fastq";
        }
        &execute("mv $alndir/* $outdir/", 1, $LOG, 2);
        &execute("rm -rf $alndir", 1, $LOG, 2);
        &execute("rm -f $SSfqin", 1, $LOG, 2);
        
        &execute("samtools view -h -F 4 $SSpref.bam | samtools view -bS - > $outpref.bam", 1, $LOG, 2);
        &execute("rm -f $SSpref.bam", 1, $LOG, 2);
        &execute("$samtools index $outpref.bam", 1, $LOG, 2);
        &execute("$samtools idxstats $outpref.bam > $outpref.idxstats.txt", 1, $LOG, 2);
        &execute("ln -sf \$(readlink -f $outpref.bam) $outpref.primary.bam", 1, $LOG, 2); 
        &execute("ln -sf \$(readlink -f $outpref.bam.bai) $outpref.primary.bam.bai", 1, $LOG, 2); 
        &execute("ln -sf \$(readlink -f $outpref.idxstats.txt) $outpref.primary.idxstats.txt", 1, $LOG, 2); 
        &execute("$mergeAlignSummary -o $outpref.align_summary.txt -a shortstack --nodirname $SSpref.Log.txt", 1, $LOG, 2);
        chomp(my $nreads = `head -2 $outpref.align_summary.txt | tail -1 | cut -f2`);
        
        
    } elsif ($aligner eq 'tophat') {
        
        
        my $THpref = "$outdir/tophat";
        unless ($params) {
            $params = "-g 1 -x 1 -N 2 --bowtie1 --segment-length 10" if $smallRNA;
            $params .= " --library-type fr-firststrand" if $stranded && !$smallRNA;
            $params .= " --no-novel-juncs --no-coverage-search";
        }
        $params .= " --transcriptome-index $gtfi" if $anno && !$txome;
        my $index = $txome ? $gtfi : $Gpref;
        $noUna = 0;
        &execute("$tophat $params -p $cores -o $alndir $index $fastq1 $fastq2", 1, $LOG, 2);  # $fastq2 need not exist
        &execute("mv -f $alndir/accepted_hits.bam $outpref.bam", 1, $LOG, 2);
        &execute("samtools index $outpref.bam", 1, $LOG, 2);
        &execute("samtools idxstats $outpref.bam > $outpref.idxstats.txt", 1, $LOG, 2);
        &execute("$bam2fq $alndir/unmapped.bam $outpref.unaligned", 1, $LOG, 3);
        &execute("mv $alndir/$_ $THpref.$_", 1, $LOG, 2) foreach qw/ deletions.bed insertions.bed junctions.bed prep_reads.info align_summary.txt logs /;  # expected Tophat outputs
        &execute("rm -rf $alndir/", 1, $LOG, 2);    # remove unmapped.bam along with it
        &execute("$samtools view -h -F 256 $outpref.bam | $samtools view -h -bS - > $outpref.primary.bam", 1, $LOG, 2);
        &execute("$samtools index $outpref.primary.bam", 1, $LOG, 2);
        &execute("$samtools idxstats $outpref.primary.bam > $outpref.primary.idxstats.txt", 1, $LOG, 2);
        &execute("$mergeAlignSummary -o $outpref.align_summary.txt -a tophat --nodirname$PE $THpref.align_summary.txt", 1, $LOG, 2);
        chomp(my $nreads = `head -2 $outpref.align_summary.txt | tail -1 | cut -f2`);
        #&execute("$stripMultireads --mode2 --tophat --index $outpref.bam", 1, $LOG, 2);
        #&execute("rename stripMultireads unique $outpref.stripmultireads.*", 1, $LOG, 2);
        #&execute("samtools idxstats $outpref.unique.bam > $outpref.unique.idxstats.txt", 1, $LOG, 2);
        
        
    } elsif ($aligner eq 'star') {
        
        
        my $STpref = "$outdir/star";
        unless ($params) {
            $params = "--outReadsUnmapped Fastx" unless $PE;
        }
        $noUna = $PE ? 1 : 0;
        &execute("STAR --outFileNamePrefix $STpref. --genomeDir $Sdir --readFilesIn $fastq1 $fastq2 --readFilesCommand zcat --runThreadN $cores --outSAMtype BAM Unsorted", 1, $LOG, 2);
        &execute("$samtools sort -@ $cores -m $mem -o $outpref.bam $STpref.Aligned.out.bam", 1, $LOG, 2);
        &execute("rm -f $STpref.Aligned.out.bam", 1, $LOG, 2);
        &execute("$samtools index $outpref.bam", 1, $LOG, 2);
        &execute("$samtools idxstats $outpref.bam > $outpref.idxstats.txt", 1, $LOG, 2);
        &execute("$samtools view -h -F 256 $outpref.bam | $samtools view -h -bS - > $outpref.primary.bam", 1, $LOG, 2);
        &execute("$samtools index $outpref.primary.bam", 1, $LOG, 2);
        &execute("$samtools idxstats $outpref.primary.bam > $outpref.primary.idxstats.txt", 1, $LOG, 2);
        &execute("mv -f $STpref.Unmapped.out.mate1 $outpref.unaligned.fastq", 1, $LOG, 2) unless $PE;
        &execute("gzip -f $outpref.unaligned.fastq", 1, $LOG, 2) if -e "$outpref.unaligned.fastq";
        &execute("$mergeAlignSummary -o $outpref.align_summary.txt -a star --nodirname$PE $STpref.Log.final.txt", 1, $LOG, 2);
        chomp(my $nreads = `head -2 $outpref.align_summary.txt | tail -1 | cut -f2`);
        
        
    } elsif ($aligner eq 'bowtie1') {
        
        
        die "$0: Aligner '$aligner' not available!\n";
        
        
    } elsif ($aligner eq 'bowtie2') {
        
        
        my $BTpref = "$outdir/bowtie2";
        my $BTtmp = "$outdir/bowtie2.tmp";
        my $infq = $PE ? "-1 $fastq1 -2 $fastq2" : "-U $fastq1";
        $noUna = $PE ? 1 : 0;
        unless ($params) {
            $params = "--un-gz $outpref.unaligned.fastq.gz" unless $PE;  # PE has separate unaligned-read extraction approach
        }
        my $idx = $txome ? $gtfi : $Gpref;
        &execute("$bowtie2 $params -p $cores -x $idx $infq -S $BTtmp.sam 2> $BTpref.align_summary.txt", 1, $LOG, 1);
        &execute("samtools view -h -F 4 $BTtmp.sam > $BTtmp.bam", 1, $LOG, 2);
        &execute("$samtools sort -@ $cores -m $mem -o $outpref.bam $BTtmp.bam", 1, $LOG, 2);
        &execute("$samtools index $outpref.bam", 1, $LOG, 2);
        &execute("$samtools idxstats $outpref.bam > $outpref.idxstats.txt", 1, $LOG, 2);
        #### FIXME: TEST $PARAMS FOR MULTIALIGN-INDUCING PARAMS
        &execute("ln -sf $outpref.bam $outpref.primary.bam", 1, $LOG, 2);
        &execute("ln -sf $outpref.bam.bai $outpref.primary.bam.bai", 1, $LOG, 2);
        &execute("ln -sf $outpref.idxstats.txt $outpref.primary.idxstats.txt", 1, $LOG, 2);
        &execute("rm -f $BTtmp.*", 1, $LOG, 2);
        #&execute("$mergeAlignSummary -o $outpref.align_summary.txt -a tophat --nodirname$PE $BTpref.align_summary.txt", 1, $LOG, 2);  ############# NOT: PE bowtie2 not enabled here
        chomp(my $nreads = `head -1 $BTpref.align_summary.txt | cut -f1 -d' '`);
        
        
    } elsif ($aligner eq 'bwa') {
        
        
        $noUna = 1;
        #$aligner2 = 'bwax' if $bwa_version > XXXXX;
        die "$0: Aligner '$aligner' not available!\n";
        
        
    } else {
        
        
        die "$0: Aligner '$aligner' not available!\n";
        
        
    }

    ($allbam, $pribam) = ("$outpref.bam", "$outpref.primary.bam");
    &logreport("Aligment complete: ".`date`, 1, $LOG);


    ## Post-alignment
    &logreport("Running post-alignment:", 1, $LOG);

    #~/lbin/CollectRnaSeqMetricsPlus -g $geno -a $anno -b $pri.bam -m $mem -o $outdir/RnaSeqMetrics.txt --old
    #~/lbin/CollectRnaSeqMetricsPlus -g $geno -a $anno -b $pri.bam -m $mem -o $outdir/RnaSeqMetrics.txt --old -s samtools-sorted

    #cufflinks -p 4 -u -G $annots.cuff.gtf -o $outdir $all.bam --max-bundle-frags 1000000

    if ($bw || $tdf) {
        my $cmd = "$bam2track -t $cores -r $mem -b ";
        $cmd .= $pribam ? $pribam : $allbam;
        $cmd .= " --primary" unless $pribam;
        $cmd .= " -g $geno -n APM";
        $cmd .= " --Mval $M" if $M;
        $cmd .= " --PE" if $PE;
        $cmd .= " --stranded --switch" if $stranded;  # assuming all stranded is TruSeq
        $cmd .= " --BW" if $bw;
        $cmd .= " --TDF" if $tdf;
        &execute($cmd, 1, $LOG, 2);
    }

    &execute("$readCount --uhisto $fastq1 $fastq2 > $outdir/fqReads.txt", 1, $LOG, 2) if $runfqr;  # $fastq2 need not exist


    ## Extract final unaligned read sets
    if ($noUna) {
        &logreport("Extracting unaligned reads: ".`date`, 1, $LOG);
        if ($PE) {
            &execute("$fastqUnaligned -fq1 $fastq1 -fq2 $fastq2 -p $outpref.unaligned $outpref.bam", 1, $LOG, 2);
        } else {
            &execute("$fastqUnaligned -fq1 $fastq1 -p $outpref.unaligned $outpref.bam", 1, $LOG, 2);
        }	
    }
    
    ## Add some bamQuantitate flags as needed
    $bamquantcall .= ' --transcriptome' if $txome;
    
}


## Finalize master-bam quantitation call
$M = $nreads unless defined $M;
$bamquantcall .= " --Mval $M";
$bamquantcall .= ' --stranded' if $stranded;
$bamquantcall .= ' --switch' if $switch;
##$bamquantcall .= ' --3pass' if ...;   # no --3pass option here; must call $bamQuantitate yourself
system "echo \"$bamquantcall -i $outpref.bam -o $outpref.counts_all.txt\" > $qsha"; 
system "echo \"$bamquantcall -i $outpref.bam -o $outpref.counts_unq.txt -n $aligner2\" > $qshu" unless $nounq;  # ribosome reference NOT designed for uniqueness; no-multiread quantitation is pointless


## Split bam on read length
if ($rlsplit) {
    &execute("$readGroupByLength $outpref.bam $outpref.RGL.bam $outpref.RGL", 1, $LOG, 2);
    &execute("mv -f $outpref.RGL.bam $outpref.bam", 1, $LOG, 2);
    &execute("samtools index $outpref.bam", 1, $LOG, 2);
    &execute("samtools idxstats $outpref.bam > $outpref.idxstats.txt", 1, $LOG, 2);
    ## Prepare readlength-bam quantitation
    foreach my $rglbam (glob "$outpref.RGL/*.bam") {
        (my $rglpref = $rglbam) =~ s/\.bam$//;
        system "echo \"$bamquantcall -i $rglbam -o $rglpref.counts_all.txt\" >> $qsha";
        system "echo \"$bamquantcall -i $rglbam -o $rglpref.counts_unq.txt -n $aligner2\" >> $qshu" unless $nounq;    # ribosome reference NOT designed for uniqueness; no-multiread quantitation is pointless
    }
}


## Quantitate regions
if ($quant) {
    &logreport("Quantitating: ".`date`, 1, $LOG);
    system "$qsha";
    system "$qshu";
}


## Exit
&logreport("Complete: ".`date`, 1, $LOG);
close $LOG;
exit;