BWA alignment

Question

[Vahid]

Hello everybody,

I wanted to align some files against the reference genome using the following script:

files="Chionobathyscus_dewitti_12
Chionobathyscus_dewitti_14
Chionobathyscus_dewitti_15
Chionobathyscus_dewitti_16
Chionobathyscus_dewitti_4
Chionobathyscus_dewitti_5
Chionobathyscus_dewitti_6
Chionobathyscus_dewitti_8
Chionobathyscus_dewitti_1
Chionobathyscus_dewitti_2
Chionobathyscus_dewitti_3
Chionobathyscus_dewitti_4
Chionobathyscus_dewitti_5
Chionobathyscus_dewitti_7
Chionobathyscus_dewitti_9
Chionobathyscus_dewitti_10
Chionobathyscus_dewitti_11
Chionobathyscus_dewitti_13"

bwa_db=GCA_943594065.1_fChiDew1_genomic.fna.gz (# reference genome)
for sample in $files
do echo $sample
    bwa mem -t 2 $bwa_db ${sample}.1.fq.gz ${sample}.2.fq.gz  |   samtools view -b | samtools sort --threads 4 > ${sample}.bam
done

However, I got the following error:

Chionobathyscus_dewitti_12
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 296298 sequences (20000115 bp)...
[main_samview] fail to read the header from "-".
[W::hts_set_opt] Cannot change block size for this format
samtools sort: failed to read header from "-"

Your insights on resolving this issue would be greatly appreciated. Cheers, Vahid

[Vahid]

Thank you very much for your assistance. The problem has been solved; it was an index issue.

[ATpoint]

That is not what the logs above tell, but good you solved it.

[a.alnawfal.1992]

that's great, but i believe the issue is related to stdin rather than index

Answer 1

[a.alnawfal.1992]

The error [main_samview] suggests that Samtools is expecting a SAM file header but isn't getting one. This often occurs when the pipe (|) between bwa mem and samtools view isn't correctly passing output due to an upstream error or empty output.

try this one:

bwa_db="GCA_943594065.1_fChiDew1_genomic.fna.gz" 
files=("Chionobathyscus_dewitti_12" "Chionobathyscus_dewitti_14" "Chionobathyscus_dewitti_15" ...) # Add you samples
for sample in "${files[@]}"
do
    echo "Processing $sample"
    bwa mem -t 2 $bwa_db ${sample}.1.fq.gz ${sample}.2.fq.gz | samtools view -b | samtools sort --threads 4 -o ${sample}.bam -
    if [ $? -eq 0 ]; then
        echo "$sample processed successfully"
    else
        echo "Error processing $sample"
    fi
done

[ATpoint]

The relevant code line is precisely the same as in the question.

[a.alnawfal.1992]

it is not!!