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,