<commit_before>#!/bin/bash

USAGE=" 

Usage: $0 [-nfh] [-c halorc] projectrc

Positional arguments
  
  projectrc  - project configuration file

Options

  -n         - dry run
  -f         - force execution of all steps
  -c         - halo runtime configuration. Default: ~/.halorc

"

# Getting the path to the current script
# http://stackoverflow.com/questions/59895/can-a-bash-script-tell-what-directory-its-stored-in
DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )"
RESYNCMATES=$DIR/resyncMates.pl

# Force move operations
MV="mv -f"

# Setup environment
DRY_RUN=false
FORCE=false

E_PARAM_ERR=250       # If less than 2 params passed to function.
P_RUN=true            # Run command
P_NORUN=false         # Skip command
STAGING_DIR=tx        # Staging directory in which to generate temporary files
JAVA_OPTS=-Xmx2g
# Initialize stdout and stderr variables; can be overloaded in config
# file. These variables set where *program output* logging goes;
# pipeline status commands will still be printed to console
LOGFILE=/dev/stdout
ERRFILE=/dev/stderr

# Halo configuration file
halorc=~/.halorc
# Check user input
while getopts :nfhc: flag; do
  case $flag in
    n)
      dry_run=true;
      ;;
    f)
      force=true;
      ;;
    c) 
      halorc=$OPTARG;
      ;;
    h)
      echo "$USAGE";
      exit;
      ;;
    \?)
      echo "Invalid option: -$OPTARG" >&2
      echo "$USAGE";
      exit;
      ;;
  esac
done
shift $(( OPTIND - 1 ));

# Exit if no input
if [ $# -ne 1 ]; then
    echo "ERROR: Wrong number of arguments"
    echo "$USAGE"
    exit;
fi

# Functions
# run_command: emulate dry run behaviour
# NOTE: currently not used!
run_command () {
    if [ -z "$1" ]
    then
	echo $(date) "No command passed to run_command"
	exit
    fi
    command=$1
    output=$2
    if [ "$dry_run" ]; then
	echo $(date) "Command: " $command >> $LOGFILE
    else
	up_to_date $output
	retval=$?
	if [ ! "$retval" ]; then
	    echo $(date) "command up to date; skipping"
	    echo -e "\t" $command >> $LOGFILE
	else
	    $command
	fi
    fi
}

# up_to_date: check if job is up to date. Very simple rules: if only
# output file given, we assume that job is done if file exists. If
# pre_output file given, compare time stamps and return 0 if output
# older than pre_output (emulates make).
up_to_date () {
    if [ $force ]; then
	echo "Force flag set; running analysis" >> $LOGFILE
	return 0
    fi
    if [ -z "$1" ]
    then
	echo "No parameter passed to up_to_date" >> $ERRFILE
	exit
    fi
    up_to_date_output=$1
    if [ ! -e $up_to_date_output ]; then
	echo "up_to_date: No such file $up_to_date_output; running analysis" >> $ERRFILE
	return 0
    fi
    if [ ! -z "$2" ]; then
	up_to_date_pre_output=$2
	if [ -e "$up_to_date_pre_output" ]; then
	    up_to_date_output_time=`stat -c %Y $up_to_date_output`
	    up_to_date_pre_output_time=`stat -c %Y $up_to_date_pre_output`
	    if [ $up_to_date_output_time -gt $up_to_date_pre_output_time ];
	    then
		return 1
	    else
		return 0
	    fi
	fi
    fi
    return 1
}

##################################################
# Start processing samples
##################################################
# Source the project runtime configuration. 
projectrc=$1
if [ ! -e $projectrc ]; then
    echo $(date) "no project configuration file found; no such file $projectrc" >> $ERRFILE
    exit
fi
source $projectrc

# Source the halo runtime configuration.
if [ ! -e $halorc ]; then
    echo $(date) "no halo configuration file found; no such file $halorc" >> $ERRFILE
    exit
fi
source $halorc

if [ ! "$samples" ]; then
    echo  $(date) "No samples provided; please provide either sample names or a file name listing sample names in the project configuration file" >> $ERRFILE
    echo "$USAGE" >> $ERRFILE
    exit
fi

if [ ! "$indir" ]; then
    echo  $(date) "No indir provided; please provide indir in the project configuration file" >> $ERRFILE
    echo "$USAGE" >> $ERRFILE
    exit
fi
# Init logging
echo $(date) > $LOGFILE
echo $(date) > $ERRFILE

# Make sure indir is absolute path
if [ "`echo $indir | cut -c1`" != "/" ]; then
    echo  $(date) "$indir is not an absolute path; please use absolute path names for input directory" >> $ERRFILE
    exit
fi

# Set input directories
if [ -f "$samples" ]; then
    echo  $(date) "Reading input directories from file $samples"
    samples=`cat $samples`
else
    echo $(date) "Assuming input directories passed as argument list"
fi
sample_regexp=`echo $samples | sed -e s"/ /\|/"g`

# Find input files based on regular expressions that include sample names
# Here assuming casava-based file names

echo "$(date) Finding files with command 'find $indir -regextype posix-extended -regex \".*(${sample_regexp}).*${READ1_REGEXP}.fastq.gz\"'"
infiles=`find $indir -regextype posix-extended -regex ".*(${sample_regexp}).*${READ1_REGEXP}.fastq.gz"`

echo $(date) "Going to run pipeline on "
for f in $infiles; do 
    echo -e "\t" ${f%.fastq.gz} 
done
read1=$infiles
read2=`echo $infiles | sed -e "s/${READ1_REGEXP}/${READ2_REGEXP}/g"`

# Setup staging directories
for f in $read1 $read2; do
    TMPDIR=`dirname $f`/tx
    mkdir -p $TMPDIR
done

##################################################
# Pipeline code
# The various steps follow a presentation on HaloPlex
# analysis obtained from Agilent
##################################################

##############################
# Primary analysis
##############################
# 1. QC
command=""
for f in $read1 $read2; do
    outdir=`dirname $f`/fastqc
    outfile=$outdir/`basename ${f%.fastq*}`_fastqc/summary.txt
    up_to_date  $outfile $f
    if [ $? = 1 ]; then continue; fi
    echo "running fastqc"  >> $LOGFILE
    mkdir -p $outdir
    cmd="$FASTQC $f -o ${outdir}"
    command="$command\n$cmd"
done
echo -e $(date) 1. QC section 
echo -e $(date) $command
echo -e "$command" | $PARALLEL  >> $LOGFILE 2>> $ERRFILE

# 2a. Trim adapter sequence
command=""
trimfiles=""
for f in $read1; do
    trimfiles="$trimfiles ${f%.fastq.gz}.trimmed.fastq.gz"
    TMP=`dirname $f`/$STAGING_DIR/`basename $f`
    up_to_date ${f%.fastq.gz}.trimmed.fastq.gz $f
    if [ $? = 1 ]; then continue; fi
    cmd="$CUTADAPT $CUTADAPT_OPTS -a $THREEPRIME $f -o ${TMP%.fastq.gz}.trimmed.fastq.gz > ${f%.fastq.gz}.trimmed.fastq.cutadapt_metrics && $MV ${TMP%.fastq.gz}.trimmed.fastq.gz ${f%.fastq.gz}.trimmed.fastq.gz"
    command="$command\n$cmd"
done
for f in $read2; do
    trimfiles="$trimfiles ${f%.fastq.gz}.trimmed.fastq.gz"
    TMP=`dirname $f`/$STAGING_DIR/`basename $f`
    up_to_date ${f%.fastq.gz}.trimmed.fastq.gz $f
    if [ $? = 1 ]; then continue; fi
    cmd="$CUTADAPT $CUTADAPT_OPTS -a $FIVEPRIME $f -o ${TMP%.fastq.gz}.trimmed.fastq.gz > ${f%.fastq.gz}.trimmed.fastq.cutadapt_metrics && $MV ${TMP%.fastq.gz}.trimmed.fastq.gz ${f%.fastq.gz}.trimmed.fastq.gz"
    command="$command\n$cmd"
done
echo -e $(date) 2a. Adapter trimming
echo -e $(date) $command
echo -e "$command" | $PARALLEL  >> $LOGFILE 2>> $ERRFILE

# 2b. Resync mates - sometimes cutadapt cuts reads down to 0, so there
# are some reads without mates
sample_pfx=`for f in $read1; do echo ${f%_${READ1_REGEXP}.fastq.gz}; done`
command=""
syncfiles=""
for f in $sample_pfx; do
    OUTDIR=`dirname $f`
    TMP=$OUTDIR/$STAGING_DIR/`basename $f`
    syncfiles="$syncfiles ${TMP}_${READ1_REGEXP}.trimmed.sync.fastq.gz ${TMP}_${READ2_REGEXP}.trimmed.sync.fastq.gz"
    up_to_date ${f}_${READ1_REGEXP}.trimmed.sync.fastq.gz ${f}_${READ1_REGEXP}.trimmed.fastq.gz
    if [ $? = 1 ]; then continue; fi
    echo $(date) resyncing reads for ${f}_${READ1_REGEXP}.trimmed.fastq.gz, ${f}_${READ2_REGEXP}.trimmed.fastq.gz
    cmd="$RESYNCMATES -i ${f}_${READ1_REGEXP}.trimmed.fastq.gz -j ${f}_${READ2_REGEXP}.trimmed.fastq.gz -o ${TMP}_${READ1_REGEXP}.trimmed.sync.fastq.gz -p ${TMP}_${READ2_REGEXP}.trimmed.sync.fastq.gz && $MV ${TMP}_${READ1_REGEXP}.trimmed.sync.fastq.gz $OUTDIR && $MV ${TMP}_${READ2_REGEXP}.trimmed.sync.fastq.gz $OUTDIR"
    command="$command\n$cmd"
done
echo -e $(date) 2b. Resync mates
echo -e $(date) $command
echo -e "$command" | $PARALLEL  >> $LOGFILE 2>> $ERRFILE

##############################
# Mapping - secondary analysis
##############################

# 3. Align sequences with bwa. Here we run command sequentially since
# bwa takes care of parallelization. From now on we run at sample level.
# Remove staging directory from syncfiles
syncfiles=`echo $syncfiles | sed -e "s/$STAGING_DIR\///g;"`
echo -e $(date) 3. Alignment
for f in $syncfiles; do
    OUTDIR=`dirname $f`
    TMP=$OUTDIR/$STAGING_DIR/`basename $f`
    up_to_date ${f%.fastq.gz}.sai $f
    if [ $? = 1 ]; then continue; fi
    echo $(date) aligning reads $f
    echo $(date) "$BWA aln -t $N_CORES $BWA_REF $f > ${TMP%.fastq.gz}.sai"
    $BWA aln -t $N_CORES $BWA_REF $f > ${TMP%.fastq.gz}.sai 2>> $ERRFILE && $MV ${TMP%.fastq.gz}.sai $OUTDIR
done

# 4. Pair reads
# NB: sam file is kept in staging directory for following step (is this really necessary?)
command=""
for f in $sample_pfx; do
    OUTDIR=`dirname $f`
    TMP=$OUTDIR/$STAGING_DIR/`basename $f`
    label=`basename $f`
    extension="trimmed.sync"
    up_to_date $f.sam ${f}_${READ1_REGEXP}.${extension}.sai
    if [ $? = 1 ]; then continue; fi
    echo $(date) pairing reads for sample $f
    cmd="$BWA sampe -A -P -r \"@RG\tID:${label}\tSM:${label}\tPL:Illumina\tCN:Agilent\" $BWA_REF ${f}_${READ1_REGEXP}.${extension}.sai ${f}_${READ2_REGEXP}.${extension}.sai ${f}_${READ1_REGEXP}.${extension}.fastq.gz ${f}_${READ2_REGEXP}.${extension}.fastq.gz > ${TMP}.sam && $MV ${TMP}.sam $OUTDIR"
    command="$command\n$cmd"
done
echo -e $(date) 4. Pair reads
echo -e $(date) $command
echo -e $command | $PARALLEL >> $LOGFILE 2>> $ERRFILE

# 5. Generate bam file
SORTSAM_OPTS="SO=coordinate"
command=""
for f in $sample_pfx; do
    OUTDIR=`dirname $f`
    TMP=$OUTDIR/$STAGING_DIR/`basename $f`
    up_to_date $f.sort.bam $f.bam
    if [ $? = 1 ]; then continue; fi
    echo $(date) generating bam file for $f
    echo "$SAMTOOLS view -bS $f.sam > ${TMP}.bam && java -jar $PICARD_HOME/SortSam.jar $SORTSAM_OPTS INPUT=${TMP}.bam OUTPUT=${TMP}.sort.bam && echo \"$f.bam removed to save space; see $f.sort.bam\" > ${TMP}.bam && $MV ${TMP}.*.ba[mi] $OUTDIR && echo \"$f.sam removed to save space; see $f.sort.bam\" > $f.sam"
    cmd="$SAMTOOLS view -bS $f.sam > ${TMP}.bam && java -jar $PICARD_HOME/SortSam.jar $SORTSAM_OPTS INPUT=${TMP}.bam OUTPUT=${TMP}.sort.bam && echo \"$f.bam removed to save space; see $f.sort.bam\" > ${TMP}.bam && $MV ${TMP}*.ba[mi] $OUTDIR && echo \"$f.sam removed to save space; see $f.sort.bam\" > $f.sam"
    command="$command\n$cmd"
done<commit_msg>Remove trailing comment<commit_after>#!/bin/bash

USAGE=" 

Usage: $0 [-nfh] [-c halorc] projectrc

Positional arguments
  
  projectrc  - project configuration file

Options

  -n         - dry run
  -f         - force execution of all steps
  -c         - halo runtime configuration. Default: ~/.halorc

"

# Getting the path to the current script
# http://stackoverflow.com/questions/59895/can-a-bash-script-tell-what-directory-its-stored-in
DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )"
RESYNCMATES=$DIR/resyncMates.pl

# Force move operations
MV="mv -f"

# Setup environment
DRY_RUN=false
FORCE=false

E_PARAM_ERR=250       # If less than 2 params passed to function.
P_RUN=true            # Run command
P_NORUN=false         # Skip command
STAGING_DIR=tx        # Staging directory in which to generate temporary files
JAVA_OPTS=-Xmx2g
# Initialize stdout and stderr variables; can be overloaded in config
# file. These variables set where *program output* logging goes;
# pipeline status commands will still be printed to console
LOGFILE=/dev/stdout
ERRFILE=/dev/stderr

# Halo configuration file
halorc=~/.halorc
# Check user input
while getopts :nfhc: flag; do
  case $flag in
    n)
      dry_run=true;
      ;;
    f)
      force=true;
      ;;
    c) 
      halorc=$OPTARG;
      ;;
    h)
      echo "$USAGE";
      exit;
      ;;
    \?)
      echo "Invalid option: -$OPTARG" >&2
      echo "$USAGE";
      exit;
      ;;
  esac
done
shift $(( OPTIND - 1 ));

# Exit if no input
if [ $# -ne 1 ]; then
    echo "ERROR: Wrong number of arguments"
    echo "$USAGE"
    exit;
fi

# Functions
# run_command: emulate dry run behaviour
# NOTE: currently not used!
run_command () {
    if [ -z "$1" ]
    then
	echo $(date) "No command passed to run_command"
	exit
    fi
    command=$1
    output=$2
    if [ "$dry_run" ]; then
	echo $(date) "Command: " $command >> $LOGFILE
    else
	up_to_date $output
	retval=$?
	if [ ! "$retval" ]; then
	    echo $(date) "command up to date; skipping"
	    echo -e "\t" $command >> $LOGFILE
	else
	    $command
	fi
    fi
}

# up_to_date: check if job is up to date. Very simple rules: if only
# output file given, we assume that job is done if file exists. If
# pre_output file given, compare time stamps and return 0 if output
# older than pre_output (emulates make).
up_to_date () {
    if [ $force ]; then
	echo "Force flag set; running analysis" >> $LOGFILE
	return 0
    fi
    if [ -z "$1" ]
    then
	echo "No parameter passed to up_to_date" >> $ERRFILE
	exit
    fi
    up_to_date_output=$1
    if [ ! -e $up_to_date_output ]; then
	echo "up_to_date: No such file $up_to_date_output; running analysis" >> $ERRFILE
	return 0
    fi
    if [ ! -z "$2" ]; then
	up_to_date_pre_output=$2
	if [ -e "$up_to_date_pre_output" ]; then
	    up_to_date_output_time=`stat -c %Y $up_to_date_output`
	    up_to_date_pre_output_time=`stat -c %Y $up_to_date_pre_output`
	    if [ $up_to_date_output_time -gt $up_to_date_pre_output_time ];
	    then
		return 1
	    else
		return 0
	    fi
	fi
    fi
    return 1
}

##################################################
# Start processing samples
##################################################
# Source the project runtime configuration. 
projectrc=$1
if [ ! -e $projectrc ]; then
    echo $(date) "no project configuration file found; no such file $projectrc" >> $ERRFILE
    exit
fi
source $projectrc

# Source the halo runtime configuration.
if [ ! -e $halorc ]; then
    echo $(date) "no halo configuration file found; no such file $halorc" >> $ERRFILE
    exit
fi
source $halorc

if [ ! "$samples" ]; then
    echo  $(date) "No samples provided; please provide either sample names or a file name listing sample names in the project configuration file" >> $ERRFILE
    echo "$USAGE" >> $ERRFILE
    exit
fi

if [ ! "$indir" ]; then
    echo  $(date) "No indir provided; please provide indir in the project configuration file" >> $ERRFILE
    echo "$USAGE" >> $ERRFILE
    exit
fi
# Init logging
echo $(date) > $LOGFILE
echo $(date) > $ERRFILE

# Make sure indir is absolute path
if [ "`echo $indir | cut -c1`" != "/" ]; then
    echo  $(date) "$indir is not an absolute path; please use absolute path names for input directory" >> $ERRFILE
    exit
fi

# Set input directories
if [ -f "$samples" ]; then
    echo  $(date) "Reading input directories from file $samples"
    samples=`cat $samples`
else
    echo $(date) "Assuming input directories passed as argument list"
fi
sample_regexp=`echo $samples | sed -e s"/ /\|/"g`

# Find input files based on regular expressions that include sample names
# Here assuming casava-based file names

echo "$(date) Finding files with command 'find $indir -regextype posix-extended -regex \".*(${sample_regexp}).*${READ1_REGEXP}.fastq.gz\"'"
infiles=`find $indir -regextype posix-extended -regex ".*(${sample_regexp}).*${READ1_REGEXP}.fastq.gz"`

echo $(date) "Going to run pipeline on "
for f in $infiles; do 
    echo -e "\t" ${f%.fastq.gz} 
done
read1=$infiles
read2=`echo $infiles | sed -e "s/${READ1_REGEXP}/${READ2_REGEXP}/g"`

# Setup staging directories
for f in $read1 $read2; do
    TMPDIR=`dirname $f`/tx
    mkdir -p $TMPDIR
done

##################################################
# Pipeline code
# The various steps follow a presentation on HaloPlex
# analysis obtained from Agilent
##################################################

##############################
# Primary analysis
##############################
# 1. QC
command=""
for f in $read1 $read2; do
    outdir=`dirname $f`/fastqc
    outfile=$outdir/`basename ${f%.fastq*}`_fastqc/summary.txt
    up_to_date  $outfile $f
    if [ $? = 1 ]; then continue; fi
    echo "running fastqc"  >> $LOGFILE
    mkdir -p $outdir
    cmd="$FASTQC $f -o ${outdir}"
    command="$command\n$cmd"
done
echo -e $(date) 1. QC section 
echo -e $(date) $command
echo -e "$command" | $PARALLEL  >> $LOGFILE 2>> $ERRFILE

# 2a. Trim adapter sequence
command=""
trimfiles=""
for f in $read1; do
    trimfiles="$trimfiles ${f%.fastq.gz}.trimmed.fastq.gz"
    TMP=`dirname $f`/$STAGING_DIR/`basename $f`
    up_to_date ${f%.fastq.gz}.trimmed.fastq.gz $f
    if [ $? = 1 ]; then continue; fi
    cmd="$CUTADAPT $CUTADAPT_OPTS -a $THREEPRIME $f -o ${TMP%.fastq.gz}.trimmed.fastq.gz > ${f%.fastq.gz}.trimmed.fastq.cutadapt_metrics && $MV ${TMP%.fastq.gz}.trimmed.fastq.gz ${f%.fastq.gz}.trimmed.fastq.gz"
    command="$command\n$cmd"
done
for f in $read2; do
    trimfiles="$trimfiles ${f%.fastq.gz}.trimmed.fastq.gz"
    TMP=`dirname $f`/$STAGING_DIR/`basename $f`
    up_to_date ${f%.fastq.gz}.trimmed.fastq.gz $f
    if [ $? = 1 ]; then continue; fi
    cmd="$CUTADAPT $CUTADAPT_OPTS -a $FIVEPRIME $f -o ${TMP%.fastq.gz}.trimmed.fastq.gz > ${f%.fastq.gz}.trimmed.fastq.cutadapt_metrics && $MV ${TMP%.fastq.gz}.trimmed.fastq.gz ${f%.fastq.gz}.trimmed.fastq.gz"
    command="$command\n$cmd"
done
echo -e $(date) 2a. Adapter trimming
echo -e $(date) $command
echo -e "$command" | $PARALLEL  >> $LOGFILE 2>> $ERRFILE

# 2b. Resync mates - sometimes cutadapt cuts reads down to 0, so there
# are some reads without mates
sample_pfx=`for f in $read1; do echo ${f%_${READ1_REGEXP}.fastq.gz}; done`
command=""
syncfiles=""
for f in $sample_pfx; do
    OUTDIR=`dirname $f`
    TMP=$OUTDIR/$STAGING_DIR/`basename $f`
    syncfiles="$syncfiles ${TMP}_${READ1_REGEXP}.trimmed.sync.fastq.gz ${TMP}_${READ2_REGEXP}.trimmed.sync.fastq.gz"
    up_to_date ${f}_${READ1_REGEXP}.trimmed.sync.fastq.gz ${f}_${READ1_REGEXP}.trimmed.fastq.gz
    if [ $? = 1 ]; then continue; fi
    echo $(date) resyncing reads for ${f}_${READ1_REGEXP}.trimmed.fastq.gz, ${f}_${READ2_REGEXP}.trimmed.fastq.gz
    cmd="$RESYNCMATES -i ${f}_${READ1_REGEXP}.trimmed.fastq.gz -j ${f}_${READ2_REGEXP}.trimmed.fastq.gz -o ${TMP}_${READ1_REGEXP}.trimmed.sync.fastq.gz -p ${TMP}_${READ2_REGEXP}.trimmed.sync.fastq.gz && $MV ${TMP}_${READ1_REGEXP}.trimmed.sync.fastq.gz $OUTDIR && $MV ${TMP}_${READ2_REGEXP}.trimmed.sync.fastq.gz $OUTDIR"
    command="$command\n$cmd"
done
echo -e $(date) 2b. Resync mates
echo -e $(date) $command
echo -e "$command" | $PARALLEL  >> $LOGFILE 2>> $ERRFILE

##############################
# Mapping - secondary analysis
##############################

# 3. Align sequences with bwa. Here we run command sequentially since
# bwa takes care of parallelization. From now on we run at sample level.
# Remove staging directory from syncfiles
syncfiles=`echo $syncfiles | sed -e "s/$STAGING_DIR\///g;"`
echo -e $(date) 3. Alignment
for f in $syncfiles; do
    OUTDIR=`dirname $f`
    TMP=$OUTDIR/$STAGING_DIR/`basename $f`
    up_to_date ${f%.fastq.gz}.sai $f
    if [ $? = 1 ]; then continue; fi
    echo $(date) aligning reads $f
    echo $(date) "$BWA aln -t $N_CORES $BWA_REF $f > ${TMP%.fastq.gz}.sai"
    $BWA aln -t $N_CORES $BWA_REF $f > ${TMP%.fastq.gz}.sai 2>> $ERRFILE && $MV ${TMP%.fastq.gz}.sai $OUTDIR
done

# 4. Pair reads
command=""
for f in $sample_pfx; do
    OUTDIR=`dirname $f`
    TMP=$OUTDIR/$STAGING_DIR/`basename $f`
    label=`basename $f`
    extension="trimmed.sync"
    up_to_date $f.sam ${f}_${READ1_REGEXP}.${extension}.sai
    if [ $? = 1 ]; then continue; fi
    echo $(date) pairing reads for sample $f
    cmd="$BWA sampe -A -P -r \"@RG\tID:${label}\tSM:${label}\tPL:Illumina\tCN:Agilent\" $BWA_REF ${f}_${READ1_REGEXP}.${extension}.sai ${f}_${READ2_REGEXP}.${extension}.sai ${f}_${READ1_REGEXP}.${extension}.fastq.gz ${f}_${READ2_REGEXP}.${extension}.fastq.gz > ${TMP}.sam && $MV ${TMP}.sam $OUTDIR"
    command="$command\n$cmd"
done
echo -e $(date) 4. Pair reads
echo -e $(date) $command
echo -e $command | $PARALLEL >> $LOGFILE 2>> $ERRFILE

# 5. Generate bam file
SORTSAM_OPTS="SO=coordinate"
command=""
for f in $sample_pfx; do
    OUTDIR=`dirname $f`
    TMP=$OUTDIR/$STAGING_DIR/`basename $f`
    up_to_date $f.sort.bam $f.bam
    if [ $? = 1 ]; then continue; fi
    echo $(date) generating bam file for $f
    echo "$SAMTOOLS view -bS $f.sam > ${TMP}.bam && java -jar $PICARD_HOME/SortSam.jar $SORTSAM_OPTS INPUT=${TMP}.bam OUTPUT=${TMP}.sort.bam && echo \"$f.bam removed to save space; see $f.sort.bam\" > ${TMP}.bam && $MV ${TMP}.*.ba[mi] $OUTDIR && echo \"$f.sam removed to save space; see $f.sort.bam\" > $f.sam"
    cmd="$SAMTOOLS view -bS $f.sam > ${TMP}.bam && java -jar $PICARD_HOME/SortSam.jar $SORTSAM_OPTS INPUT=${TMP}.bam OUTPUT=${TMP}.sort.bam && echo \"$f.bam removed to save space; see $f.sort.bam\" > ${TMP}.bam && $MV ${TMP}*.ba[mi] $OUTDIR && echo \"$f.sam removed to save space; see $f.sort.bam\" > $f.sam"
    command="$command\n$cmd"
done
{'fraction_non_alphanumeric': 0.17876489707475623, 'fraction_numerical': 0.008007913702953506, 'mean_word_length': 3.620239390642002, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 119, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7222412', 'n_tokens_mistral': 7901, 'n_tokens_neox': 7646, 'n_words': 2174}