5
$\begingroup$

I've been running bwa mem -a for alignment, using the -a flag---this will

output all alignments for SE or unpaired PE

I've noticed in the SAM that there are several alignments with * in the SEQ and QUAL fields. Based on the documentation:

  1. SEQ: segment SEQuence. This field can be a ‘*’ when the sequence is not stored. If not a ‘*’, the length of the sequence must equal the sum of lengths of M/I/S/=/X operations in CIGAR. An ‘=’ denotes the base is identical to the reference base. No assumptions can be made on the letter cases.
  2. QUAL: ASCII of base QUALity plus 33 (same as the quality string in the Sanger FASTQ format). A base quality is the phred-scaled base error probability which equals −10 log10 Pr{base is wrong}. This field can be a ‘*’ when quality is not stored. If not a ‘*’, SEQ must not be a ‘*’ and the length of the quality string ought to equal the length of SEQ.

it appears the sequence isn't stored.

I would strongly prefer to have the sequence in this case. Is there any to direct bwa to output these sequences?

Improve this question
$\endgroup$
2
  • $\begingroup$ Are you aware that many of these records are secondary/supplementary alignments? Does your workfow/model require/benefit from these alignments and handle them correctly? $\endgroup$
    – Daniel Standage
    Commented Dec 5, 2018 at 15:23
  • 2
    $\begingroup$ @DanielStandage "Are you aware that many of these records are secondary/supplementary alignments?" Yes. " Does your workfow/model require/benefit from these alignments and handle them correctly?". Benefit, yes. $\endgroup$
    – EB2127
    Commented Dec 5, 2018 at 15:59

1 Answer 1

Reset to default
5
$\begingroup$

I don't believe that is an option. Internally we created a patch and submitted it to BWA however that was 4 years ago and it has not been accepted.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks. Did @user172818 ever give a reason why this feature isn't included? $\endgroup$
    – EB2127
    Commented Dec 5, 2018 at 3:29
  • 1
    $\begingroup$ see also biostars.org/p/352930 $\endgroup$
    – Pierre
    Commented Dec 5, 2018 at 9:26
  • $\begingroup$ @Pierre This is a great tool Pierre! I think it runs successfully for me. The problem is, I get errors downstream when using samtools view, i.e.. [E::sam_parse1] CIGAR and query sequence are of different length Are there work arounds? $\endgroup$
    – EB2127
    Commented Dec 5, 2018 at 19:13
  • $\begingroup$ I didn't test it with all the possible cases. Can you please try to isolate a minimal sam file containing the faulty reads ( all the reads with the same name) and send a bug report ? github.com/lindenb/jvarkit/issues $\endgroup$
    – Pierre
    Commented Dec 5, 2018 at 19:54
  • $\begingroup$ @Pierre I may send you an e-mail instead. Thanks Pierre! $\endgroup$
    – EB2127
    Commented Dec 7, 2018 at 5:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.