BWA alignment
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
• 2,478 views
•
link
updated
by
a.alnawfal.1992
•
written
by
Vahid •
•
written
by
Vahid •
1 answer
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
• 0 views
•
link
written
by
a.alnawfal.1992
it is not!!
• 0 views
•
link
written
by
a.alnawfal.1992
Log in to answer this question.
Thank you very much for your assistance. The problem has been solved; it was an index issue.
That is not what the logs above tell, but good you solved it.
that's great, but i believe the issue is related to stdin rather than index