Wrapper script for paired-end segemehl

map_em_segemehl_pe.sh

#!/bin/bash

segemehl="/scratch/mtw/src/segemehl-svn/segemehl/segemehl.x"
results="./alignments"
reffa="../../genomes/TAIR10/TAIR10.fa"
refidx="../../genomes/TAIR10/TAIR10.idx"
samples=2
rep=3

if [ -d "$results" ];
then
        echo "$results available"
else
        mkdir $results
fi

if ! [ -f "$reffa" ]
then
        echo "Reference FASTA not found"
        exit 2
fi

if ! [ -f "$refidx" ]
then
        echo "Reference Index not found"
        exit 2
fi

for s in $(seq 1 $samples)
do
        for r in $(seq 1 $rep)
        do
        sample=C${s}_R${r}
        echo processing $sample
        rd1="../raw/${sample}_read_1.trim.fq.gz"
        rd2="../raw/${sample}_read_2.trim.fq.gz"
        set -x
        $segemehl -q $rd1 -p $rd2 -u $results/${sample}.unmatched -d $reffa -i $refidx -s -S -t 30 2> $results/${sample}.segemehl.err | samtools view -Sb - | samtools sort -o - ${sample} - > $results/${sample}.bam
        set +x
        done
done