for loop over bwa aln

Question

[phylofun]

My for loop around bwa aln and bwa sampe doesn't seem to work. First, I created paths to the read data:

ls -1r  reads/*/cleaned/*READ1.fastq.gz > READ1_list
ls -1r  reads/*/cleaned/*READ2.fastq.gz > READ2_list

Next, I made a for loop using READ1_list and READ2_list around the bwa aln function:

for i in READ1_list 
do
  bwa aln reference/reference.unaligned.fasta $i > ${i%READ1.fastq.gz}read1.sai 
done

for i in READ2_list
do
  bwa aln reference/reference.unaligned.fasta $i > ${i%READ2.fastq.gz}read1.sai 
done

Then, I use a for loop around bwa sampe:

 for i in $(cat READ1_list | sed s'/\-READ1.fastq.gz//' | sort) 
do
  bwa sampe reference/reference.unaligned.fasta $i\-READ1.sai $i\-READ2.sai $i\-READ1.fastq.gz $i\-READ2.fastq.gz \
  | samtools view -bS - > $i\-pe.bam
done

Can anyone spot the issue?

[h.mon]

I can't spot the issue. What happens? Are there error messages?

Unless you reads are shorter than 70bp, you should be using bwa mem instead of bwa aln.

[phylofun]

There is an error message when I run bwa sampe, but I think the issue begins with bwa aln. No error message occurs when I run bwa aln; however, the .sai output file goes to the current working directory instead of following the path in READ1_list, and the sample name is not written for the .sai file (e.g. READ1_listread1.sai instead of sample_name-read1.sai)

[WouterDeCoster]

Put an echo in your loop to see how it looks like:

 for i in $(cat READ1_list | sed s'/\-READ1.fastq.gz//' | sort) 
do
  echo "bwa sampe reference/reference.unaligned.fasta $i\-READ1.sai $i\-READ2.sai $i\-READ1.fastq.gz $i\-READ2.fastq.gz \
  | samtools view -bS - > $i\-pe.bam"
done

Please note that you could also do the following and pipe directly to samtools sort, without samtools view:

bwa sample ref.fa R1 R2 | samtools sort -o alignment.bam

[JC]

I see two problems in your code:

1. The second alignment outputs a read1.sai file, it should be read2.sai
2. The final BAM file requires files named as READ but your alignments are read (file names are generally case-sensitive)

Also it is hard to find the errors without any detail, are you running this on Linux/Mac or Windows?, what are the names of the input files?, what is inside the READ1_list and READ2_list files?

[5heikki]

for i in READ1_list

That shouldn't even work. You should instead e.g.

for i in $(cat READ1_list); do ...

Just try e.g.

for i in READ1_list; do echo "$i"; done

vs.

for i in $(cat READ1_list); do echo "$i"; done

Only the latter works in bash (as intended)..

But anyway, I wouldn't even use a for loop for this but instead just define a function and use it with find and GNU parallel. Also, why are you aligning the reads separately? Surely it would make more sense to align read pairs properly in the same runs?

Answer 1

[Joe]

This is overkill, but you might be interested in my bash code for aligning/remapping: (You'll mainly want to look at Line 95 onward, everything before that is just picking up commandline args etc.)

With a little more bash magic, you can loop your files and pass the relevant paths as arguments:

paste -- "READS1.txt" "READS2.txt" |
while IFS=$'\t' read -r reads1 reads2; do
 bash align_and_remap.sh -r /path/to/reference.fasta -1 "$reads1" -2 "$reads2" -o /path/to/outdir
done

*Disclaimer, this script is not well tested at all. Do some testing of you own before letting it loose on all your data. This is also based on an old setup with bwa aln where you had to provide each read file separately. That doesn't need to happen with bwa mem I believe.

	#!/bin/bash
	# Script to align fastq reads to a given reference, and output the necessary files sorted etc.
	# Script requires bwa and samtools in path - check for existence:
	command -v bwa >/dev/null 2>&1 || { echo >&2 "BWA doesn't appear to be installed. Aborting."; exit 1; }
	command -v samtools >/dev/null 2>&1 || { echo >&2 "samtools doesn't appear to be installed. Aborting."; exit 1; }
	# Capture inputs
	usage()
	{
	cat << EOF
	usage: $0 options
	This script aligns fastq reads to a given
	reference, and outputs alignment maps etc.
	OPTIONS:
	-h | --help Show this message
	-r | --ref Reference fasta sequence
	-1 | --read1 First set of reads ("R1")
	-2 | --read2 Second set of reads ("R2")
	-o | --outdir Directory to output all the files to
	EOF
	}
	# Tolerate long arguments
	for arg in "$@"; do
	shift
	case "$arg" in
	"--help")
	set -- "$@" "-h"
	;;
	"--ref")
	set -- "$@" "-r"
	;;
	"--read1")
	set -- "$@" "-1"
	;;
	"--read2")
	set -- "$@" "-2"
	;;
	"--outdir")
	set == "$@" "-o"
	;;
	*)
	set -- "$@" "$arg"
	esac
	done
	# getopts assigns the arguments to variables
	while getopts "hr:1:2:o:" OPTION
	do
	case $OPTION in
	r)
	reference=$OPTARG
	;;
	1)
	read1=$OPTARG
	;;
	2)
	read2=$OPTARG
	;;
	o)
	outdir=$OPTARG
	;;
	h)
	usage
	exit
	;;
	esac
	done
	if [[ -z $reference ]]
	then
	usage
	echo "No Reference sequence provided. Exiting." ; exit 1
	fi
	if [[ -z $read1 ]]
	then
	usage
	echo "Forward reads not supplied. Exiting." ; exit 1
	fi
	if [[ -z $read2 ]]
	then
	usage
	echo "Reverse reads not supplied. Exiting." ; exit 1
	fi
	if [[ -z $outdir ]]
	then
	outdir=$(pwd)
	echo "No output directory was specified. Defaulting to the current working directory:"
	pwd
	fi
	#####
	reads=("$read1" "$read2")
	# Index the reference sequence
	bwa index "$reference"
	echo "BWA FINISHED INDEXING THE REFERENCE"
	# Align reads
	alignedreads=()
	for readelement in ${reads[@]} ; do
	bwa mem "$reference" "$readelement" > ${readelement%.*.*}.sai
	echo "BWA FINISHED ALIGNING $readelement "
	alignedreads+=("${readelement%.*.*}.sai")
	done
	echo "BWA FINISHED ALIGNING EACH READ."
	bwa sampe "$reference" "${alignedreads[0]}" "${alignedreads[1]}" "$read1" "$read2" > "${reference%.*}"_aligned.sam
	echo "SAM FILE CREATED"
	samtools view -S "${reference%.*}"_aligned.sam -b -o ./"${reference%.*}"_aligned.bam
	echo "SAM FILE CONVERTED TO BAM"
	samtools sort "${reference%.*}"_aligned.bam "${reference%.*}"_aligned_sorted
	echo "BAM FILE SORTED"
	samtools index "${reference%.*}"_aligned_sorted.bam
	echo "BAM FILE INDEXED"

view raw align_and_remap.sh hosted with ❤ by GitHub