Dealing With Manta Limitations

Question

I am trying to figure out a way to remove the limiting factors that Manta says it cannot handle. To quote from the Manta page:

The following limitations exist on the input BAM or CRAM files provided to Manta:

• Alignments cannot have an unknown read sequence (SEQ="*")
• Alignments cannot contain the "=" character in the SEQ field.
• Alignments cannot use the sequence match/mismatch ("="/"X") CIGAR notation RG (read group) tags in the alignment records are ignored -- each file will be treated as representing one sample.
• Alignments with basecall quality values greater than 70 are rejected (these are not supported on the assumption that this indicates an offset error)

Would using a find statement in this way and then making .bai files be adequate to remove the limiting factors that Manta has?:

find ${input_dir} -name "*.bam" -type f -exec sh -c 'output_file="${0%.bam}_filtered.bam"; temp_file=$(mktemp); samtools view -h -q 70 "$0" | sed "/^\(@\|[^[:blank:]]\+[[:blank:]]\+4\)/!s/\t[=*]/\t/g" | samtools view -Sb - > "$temp_file"; mv "$temp_file" "$output_file"; samtools index "$output_file"; rm "$0"' {} \;
# Make the .bai files:

for file in ${input_dir}/*.bam; do
    samtools index "${file}" "${file%.*}.bai"
done

However, I get the error (repeatedly for all my bam files):

[main_samview] fail to read the header from "-". samtools index: "/path/to/MyBatch1And2/SampleNumber_filtered_filtered_filtered.bam" is in a format that cannot be usefully indexed

I am not sure why I got this error once and then it did not occur again. Why? Or would this not be enough remove the limiting factors (as there are more to be considered or the code is not written correctly) of Manta?

Answer 1

Assuming you have paired-end Illumina reads mapped using BWA or BWA-MEM (with default parameters), these limitations are usually not much of a concern:

• Alignments cannot have an unknown read sequence (SEQ="*")

Best practice for representing alignments in SAM format, suggests that the SEQ and QUAL of secondary alignments should be set to an asterisk to reduce file size. You can either filter these out after the fact, using something like samtools view -F 256 <bam>, or simply avoid writing out any secondary alignments in the first place.

---

• Alignments cannot contain the "=" character in the SEQ field.

An "=" character can be used to denote that the base is identical to the reference base. I've not come across an aligner or downstream process that produces such output. If you know of such a tool or process, please comment below to let me know.

---

• Alignments cannot use the sequence match/mismatch ("="/"X") CIGAR notation

Please see: Why most aligners do not output the "X" CIGAR operation?

---

• RG (read group) tags in the alignment records are ignored -- each file will be treated as representing one sample.

If you have multiple libraries for a single sample, and they differ in some not so small way (perhaps you have a paired-end library and an old mate-pair library, or maybe the distribution of insert sizes/template lengths of each library varies wildly), you may want to consider running each library through separately.

---

• Alignments with basecall quality values greater than 70 are rejected (these are not supported on the assumption that this indicates an offset error)

Short-read sequencing technologies do not produce such high base quality scores. Only long-read technologies, like PacBio and Oxford Nanopore, produce base quality scores above 70 (up to 93).

---

---

Ultimately, I don't think you really need to deal with any of these with some command pipeline if you choose your inputs wisely.