#! /bin/bash
### qsub -q normal gsnapset.sh
#PBS -N gsnapdebug
#PBS -l mem=168gb,nodes=2:ppn=32,walltime=9:55:00
#PBS -o gsnapdebug.$$.out
#PBS -e gsnapdebug.$$.err
#PBS -V

# ncpu=32
ncpu=64
npart=64
jstart=0
jend=$(( $jstart + $ncpu ))

datad=/N/dc/scratch/$USER
workd=$datad/chrs/cacao
bindir=$HOME/bio/bin
gmapd=$HOME/bio/gmap118
rund=$workd/rnas
gdb=gmap118
gdbpath=$workd/genome/

# snpd=cacao11snp10
# dgenome=cacao11allasm ; gtag=mars11

####dgenosize=$dgenome.chr_size.txt

# .. test run thru gene.cds.fa only
snpd=snp10pub3hcds
dgenome=genes3h_goodcds ; gtag=cds3h

#1 snapopt1="--use-snps=$snpd -N 1 --quality-protocol=illumina "
#2 snapopt1="--use-snps=$snpd --show-refdiff --print-snps -N 1 --quality-protocol=illumina "
# this works: --print-snps  --show-refdiff  << print-snps lists snps: nn@idsnp,...
## stranded??  --antistranded-penalty=1
# cgb2 reads are diff format.. not same qual, -stranded, lacking pair /1 /2 ..
# also has Ill chastity string : --filter-chastity= none,either,both << this works
snapopt1="--use-snps=$snpd --show-refdiff --print-snps -N 1 --antistranded-penalty=1 --filter-chastity=both "

#...
## tie outdir to fastdir or not?
outdir=$workd/rnas/snp4$gtag
#...

notef=$workd/rnas/gsnapset.$$.RUNNING
donef=`echo $notef | sed 's/RUNNING/DONE/'`

cd $workd/rnas/
if [ "X$fastdir" = "X" ]; then
  fastdir="fastq"
fi

touch $notef
echo "START " >> $notef
echo `date`   >> $notef

mkdir $outdir
ls -l $fastdir >> $notef

#... START LOOP rnain : limit to suffix .fq, .fastq, sequence.txt ..

for rnain in $fastdir/*{.1.fastq,_1_sequence.txt}.gz ;  do 
{ 
  cd $rund/
  # mkdir $rund/sams
  if ! test -f "$rnain" ; then continue; fi
  
  drna=`basename $rnain .gz | sed 's/.gz//; s/.fastq//; s/.fq//; s/_sequence.txt//;'`
  outna=$drna-$gtag
  echo "START gsnap : $rnain to $outna.bam" >> $notef
  echo `date`  >> $notef

  snapopts="--gunzip $snapopt1"
  
  rnain2=`echo $rnain | sed 's/.1.fastq/.2.fastq/; s/_1_sequence/_2_sequence/;'`
  if test -f "$rnain2" ; then
    inset="$rnain $rnain2"
  else 
    inset=$rnain
  fi

  # to narrow down break point, use i, j counters, 
  # i=1,ncpu, j=jstart,jend  subset of npart
  i=0;  
  j=$jstart;
  while [ $i -lt $ncpu ]; do { 
  
   $gmapd/bin/gsnap $snapopts --part=$j/$npart \
    -D $gdbpath/$gdb -d $dgenome $inset > $outdir/$drna.gsnap$i.out &
  
   i=$(( $i + 1 ))
   j=$(( $j + 1 ))
  }
  done
  echo "wait gsnap : $drna" >> $notef
  
  wait

  echo "DONE gsnap : $drna" >> $notef
  du -h $outdir >> $notef
  ls -l $outdir >> $notef
  #..................
  ##...... skip sam sort/merge........
} 
done

#... END LOOP rnain
echo "DONE " >> $notef
echo `date`  >> $notef
mv $notef $donef