After surveying the existing tools for de-multiplexing barcoded paired-end sequencing reads, I grew frustrated and rolled my own solution. The issue with most of the tools for de-multiplexing reads from massively parallel sequencing is that they operate on fastq files. Since paired-end sequencing typically generates two fastq files (one for forward reads and one for reverse reads), it becomes more difficult to apply existing single-end read tools without doing some extra matching and filtering. Most people seem to either:
- Merge the forward and reverse reads into one “megaread”, demultiplex these based on barcode sequence, and then split the resulting reads back into forward and reverse before mapping.
- Filter the forward reads based on barcodes, and take advantage of same sorting order of forward and reverse reads to match pairs from the reverse reads.
Both of these routes involve creating many individual fastq files that will be individually mapped. Depending on the aligner and amount of sequence, mapping something like 75 de-multiplexed sets of reads could be inefficient since you would be initializing the aligner and indexed reference genome 75 times. This does not seem like an optimal solution.
Since my immediate purpose is amplicon resequencing on a
MiSeq, the number of reads I will be dealing with is fairly low, so I think I can design a more logical workflow. Ideally, I would like to map all the reads at once, moving the barcode from each read into the SAM
BC tag. Then I can split the resulting mapped SAM file into de-multiplexed files for analysis.
First off, mapping my reads with BWA. I’m using
Bpipe for pipeline management.
@Transform("sam")
bwa_aln_bc = {
exec "bwa aln -t $threads -B $length $bwa_reference $input1 > ${input1}.sai"
exec "bwa aln -t $threads $bwa_reference $input2 > ${input2}.sai"
exec "bwa sampe $bwa_reference ${input1}.sai ${input2}.sai $input1 $input2 > $output"
} The above function truncates $length number of bases from each forward read, and assigns that as a BC tag in the resulting SAM alignment e.g.: “BC:Z:TTAATGC”.
Next, we will split the barcoded SAM file into multiple files, based on the barcodes, while preserving the SAM header. All reads that do not match a barcode in the supplied table will be written to a separate file.
#!/bin/sh
if [ $# -eq 0 ] ; then
echo 'Usage: splitSam.sh input.sam barcodes.txt'
echo ''
echo 'barcodes.txt must be two columns with tab delimeter'
echo 'column1 = barcode name, column2 = barcode sequence'
fi
SAM=$1
BC_FILE=$2
#Capture the SAM header
SAM_H=`samtools view -SH $SAM`
#Read barcode table into array line by line
#then grep barcoded SAM reads to file
while IFS=$'\t' read -r -a array
do
sampleName=${array[0]} #barcode name
BC=${array[1]} #barcode sequence
BC=${BC%"${BC##*[![:space:]]}"} #remove trailing whitespace
if grep -q -m 1 -e "BC:Z:$BC" $SAM; then
printf "$SAM_H\n" > ${sampleName}.sam #write SAM header to file
grep -e "BC:Z:$BC" $SAM >> ${sampleName}.sam #write barcoded reads
fi
printf "BC:Z:$BC\n" >> /tmp/bc #patterns for unmatched reads
done < $BC_FILE
#Write unmatched reads to file
grep -v -f /tmp/bc $SAM > unmatched.sam #write unmatched reads
rm /tmp/bcThis approach, while efficient for a small number of reads, may be inappropriate for larger projects.
