gistfile1

#!/bin/bash
#
set -e
#set -x

dist=`pwd`
#fq="$dist/reads/ecoli.reads.fastq"
#fq="$dist/reads/reads.problem_zlib.fastq"
fq="$dist/reads/hybrid.fastq"
ref_h="/stornext/snfs4/next-gen/solid/bf.references/h/hsap.36.1.hg18/hsap_36.1_hg18.fa"
ref_e="$dist/bfast.indexes/ecoli/ecoli.drio.test.fa"
ref_e_s="$dist/bfast.indexes/one.ecoli/ecoli.drio.test.fa"
ref=$ref_e
ref_bwa="/Users/drio/projects/synthetic.pipe/ecoli/bwa.indexes/ecoli.drio.test.fa"
dwg="/Users/drio/projects/bfast_related/dnaa/dwgsim/dwgsim"

dist="/Users/drio/projects/bfast_hybrid/bfast.git"
bbin="$dist/bfast/bfast"
#bbin="/Users/drio/projects/bfast_related/versions/bfast-0.6.4d/bfast/bfast"

# Recompile if necessary
cd $dist
make
#make distclean; sh ./autogen.sh; ./configure ; make ; make check
cd -

# Change to output dir 
o_dir=output
mkdir -p $o_dir
cd $o_dir
echo "$bbin"

# Generate reads
$dwg -c -e 0.08 -N 20000 -1 50 -2 25 $ref ./bwa.read1.fastq ./bwa.read2.fastq ./bf.fastq 2> /dev/null > /dev/null

# Separate R1/R2 in separate fastq files
cat bf.fastq | ruby -ne 'BEGIN{i=1}; puts $_ if i <= 4; i = i == 8 ? 1 : i + 1' > ./50.fastq
cat bf.fastq | ruby -ne 'BEGIN{i=1}; puts $_ if i  > 4; i = i == 8 ? 1 : i + 1' > ./25.fastq

: << COMMENT_2
# Generate raw solid data from fastqs
(
cat<<END_RUBY
#!/usr/bin/env ruby1.9
#

if ARGV[0].nil?
  read = 1
else
  raise "Invalid argument (1||2)" unless ARGV[0] =~ /(1|2)/
  read = ARGV[0].to_i
end

i=0
r_name = seq = qual = ""
STDIN.each_line do |l|
  i = i + 1
  case i
    when 1,5
      r_name = l.gsub(/^@/, '').chomp
    when 2,6
      seq    = l.chomp
    when 4,8
      qual   = l.chomp
      if read == 1 and i == 4 # read 1
        \$stdout.printf ">#{r_name}\n#{seq}\n" 
        \$stderr.printf ">#{r_name}\n#{qual}\n"
      elsif read == 2 and i == 8 # read 2
        \$stdout.printf ">#{r_name}\n#{seq}\n" 
        \$stderr.printf ">#{r_name}\n#{qual}\n"
      end
      i = 0 if i == 8
  end
end
END_RUBY
) > ./fastq2raw.rb
chmod 755 ./fastq2raw.rb

cat bf.fastq | ./fastq2raw.rb 1 >  50.csfasta 2> 50.qual
cat bf.fastq | ./fastq2raw.rb 2 >  25.csfasta 2> 25.qual
COMMENT_2

# Fix the read names so they match (bwa is -1 coordinate)
: << COMMENT
# BF:  @Chr1_1124762_1124213_1_1_6:0:0_3:0:0_0
# BWA: @Chr1_1124763_1124214_1_1_6:0:0_3:0:0_0/1
cat bwa.read2.fastq | \
ruby -ne '$_ = $_.split("_").inject([]){|t, c| t << (t.size == 2 || t.size == 1 ? c.to_i - 1 : c) }.join("_") if $_ =~ /^@/; puts $_' > tmp
mv tmp bwa.read2.fastq
COMMENT

# Match different modes
echo "Match 50 (BF)"
$bbin match -T/tmp/ -n8 -i1 -A1 -l -f $ref -r 50.fastq > 50.bmf 2> match.50.log 
echo "to bmf (50)"; $bbin bmfconvert -O1 50.bmf &>/dev/null 
echo "Match 25 (BF)"
$bbin match -T/tmp/ -n8 -i1 -A1 -l -f $ref -r 25.fastq > 25.bmf 2> match.25.log
echo "to bmf"; $bbin bmfconvert -O1 25.bmf &>/dev/null
echo "Match 50+25 (BF)"
$bbin match -T/tmp/ -n8 -i1 -A1 -l -f $ref -r bf.fastq > 50_25.bmf 2> match.50_25.log
echo "Match 25 (BWA)"
$bbin bwaaln -c -n3 -t8 $ref_bwa 25.fastq > bwa.25.bmf 2> match.25.bf-bwa.log
echo "to bmf (BWA 25)"; $bbin bmfconvert -O1 bwa.25.bmf &>/dev/null

# localalign different modes
echo "Local (25)"
$bbin localalign -A1 -t -f $ref -n 8 -m 25.bmf > 25.baf 2> la.25.log
echo "Local (BWA 50)"
$bbin localalign -U -A1 -t -f $ref -n 8 -m bwa.25.bmf > la.bwa.25.baf 2> la.bwa.25.log
echo "Local (hybrid)"
$bbin localalign -U -A1 -t -f $ref -n 8 -O 50.bmf -T bwa.25.bmf > la.hybrid.baf 2> la.hybrid.log
echo "Local (regular MP)"
$bbin localalign -A1 -t -f $ref -n 8 -m 50_25.bmf > la.50_25.baf 2> la.50_25.log

# PP different modes
echo "pp (hybrid)"
$bbin postprocess -f $ref -i la.hybrid.baf -a 3 -O 1 > out.hybrid.sam  2> post.hybrid.log
echo "pp (regular BF)"
$bbin postprocess -f $ref -i la.50_25.baf -a 3 -O 1 > out.50_25.sam  2> post.50_25.log

# stats
echo "stats (hybrid)"
echo "--------------"
samtools view -b -S out.hybrid.sam > out.hybrid.bam 2>/dev/null
/Users/drio/projects/drd.bio.toolbox/bam.stats.sh -i out.hybrid.bam -e 2 > stats.hybrid.txt 2>/dev/null
cat stats.hybrid.txt

echo "stats (Regular)"
echo "--------------"
samtools view -b -S out.50_25.sam > out.50_25.bam 2>/dev/null
/Users/drio/projects/drd.bio.toolbox/bam.stats.sh -i out.50_25.bam -e 2 > stats.50_25.txt 2>/dev/null
cat stats.50_25.txt

#samtools flagstat out.bam