samtools view: writing to standard output failed: Broken pipe

Question

I have the following alignment script which uses BWA:

#!/bin/bash
#usage: sh align_pbs.sh /scratch/waterhouse_team/benth/dbg2olc-40x/good/ output.racon2.fasta bam-soap

mkdir $3

for r1 in $(find $1 -name "*R1*.gz");
do
  base=${r1%_R1*}
  output=$(basename $(echo $r1 | sed 's/R1//g'))
  r2=$(echo $r1 | sed 's/R1/R2/g') 

  #Get read group infomration:
  #Source: https://www.biostars.org/p/280837/#310132
  header=$(zcat $r1 | head -n 1)
  id=$(echo $header | head -n 1 | cut -f 1-4 -d":" | sed 's/@//' | sed 's/:/_/g')
  sm=$(echo $header | head -n 1 | grep -Eo "[ATGCN]+$")
  #echo "Read Group @RG\tID:$id\tSM:$id"_"$sm\tLB:$id"_"$sm\tPL:ILLUMINA"
  echo "@RG\tID:${output}\tSM:${output}\tLB:${output}\tPL:ILLUMINA"

  #cat <<EOF
  qsub <<EOF
#!/bin/bash -l

#PBS -N $output
#PBS -l walltime=48:00:00
#PBS -j oe
#PBS -l mem=70G
#PBS -l ncpus=12
#PBS -M [email protected]
##PBS -m bea

cd \$PBS_O_WORKDIR

conda activate bwa-hic
bwa mem \
$2 $r1 $r2 \
-t 12 \
-R "$(echo "@RG\tID:${output}\tSM:${output}\tLB:${output}\tPL:ILLUMINA")" | samblaster -r | samtools view -@ 12 -h -f 0x2 -F 2316 - | samtools fixmate - - | sed 's/^/LP1-/' | samtools sort -@ 12 - -o ${3}/${output}.sorted.dedup.bam

EOF

done

Unfortunatelye, this pipeline cuased the following erros:

samblaster: Version 0.1.24
samblaster: Inputting from stdin
samblaster: Outputting to stdout
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 1200000 sequences (120000000 bp)...
[M::process] read 1200000 sequences (120000000 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (2, 316829, 6, 4)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (432, 474, 509)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (278, 663)
[M::mem_pestat] mean and std.dev: (471.59, 53.66)
[M::mem_pestat] low and high boundaries for proper pairs: (201, 740)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 1200000 reads in 615.287 CPU sec, 58.490 real sec
samblaster: Loaded 19 header sequence entries.
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
sed: couldn't write 567 items to stdout: Broken pipe
[E::bgzf_flush] File write failed (wrong size)
samtools fixmate: Couldn't write to output file: Broken pipe
[E::bgzf_close] File write failed
[bam_mating] error while closing output file
samtools view: writing to standard output failed: Broken pipe
samtools view: error closing standard output: -1
samblaster: Unable to write to output file.

Thank you in advance,

Answer 1

[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting

SAMtools sort has been unable to parse its input, which it thought was SAM (mostly because it couldn't be recognised as another format e.g. BAM).

This is because sed 's/^/LP1-/' is putting LP1- at the front of every line. If the output of samtools fixmate is SAM, then this LP1 is garbling the SAM header lines. If the output of samtools fixmate is BAM (as it probably is), this LP1 is completely garbling the data.

What were you trying to achieve with this sed command?