2
$\begingroup$

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?

$\endgroup$
0

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
13
  • 1
    $\begingroup$ As people who don't know how to code and won't choose the inputs as carefully -- will be using my pipeline I cannot afford not to deal with it. $\endgroup$
    – Indira
    Mar 30, 2023 at 16:35
  • $\begingroup$ @Indira For a pipeline, I think your best bet would be to start upstream where you get to stipulate the requirements and choose the aligner. The requirements could be one or more pairs of FASTQ files for example. What you want to do is to make sure to fail early. You really do not want to spend several hours or so filtering out all of the long reads from an ONT BAM, only to then have it fall over when the filtered BAM no longer contains any usable data. If you can't start upstream for whatever reason, documenting the requirements might eliminate the majority of issues that arise. $\endgroup$
    – Steve
    Mar 30, 2023 at 23:52
  • $\begingroup$ @Indira BTW, I think the problem with your find command is that the sed component is failing. Using single quotes instead of double quotes around the sed expression should solve the problem. Also consider moving the pipeline command inside a separate script file and prepend set -euo pipefail to the start of it. The former should make it a bit more readable, and the latter will catch the errors. $\endgroup$
    – Steve
    Mar 31, 2023 at 0:01
  • $\begingroup$ What do you mean start upstream? I am only provided the bams? $\endgroup$
    – Indira
    Apr 11, 2023 at 13:29
  • $\begingroup$ I tried to edit the script as well. $\endgroup$
    – Indira
    Apr 11, 2023 at 14:01

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.