Answers.Body: <p>Assuming you have <strong>paired-end</strong> Illumina reads mapped using BWA or BWA-MEM (with default parameters), these limitations are usually not much of a concern:</p>
  <blockquote>
  <ul>
  <li>Alignments cannot have an unknown read sequence (SEQ=&quot;*&quot;)</li>
  </ul>
  </blockquote>
  <p>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 <code>samtools view -F 256 &lt;bam&gt;</code>, or simply avoid writing out any secondary alignments in the first place.</p>
  <hr />
  <blockquote>
  <ul>
  <li>Alignments cannot contain the &quot;=&quot; character in the SEQ field.</li>
  </ul>
  </blockquote>
  <p>An &quot;=&quot; 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.</p>
  <hr />
  <blockquote>
  <ul>
  <li>Alignments cannot use the sequence match/mismatch (&quot;=&quot;/&quot;X&quot;) CIGAR notation</li>
  </ul>
  </blockquote>
  <p>Please see: <a href="https://bioinformatics.stackexchange.com/questions/404/why-most-aligners-do-not-output-the-x-cigar-operation">Why most aligners do not output the &quot;X&quot; CIGAR operation?</a></p>
  <hr />
  <blockquote>
  <ul>
  <li>RG (read group) tags in the alignment records are ignored -- each file will be treated as representing one sample.</li>
  </ul>
  </blockquote>
  <p>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.</p>
  <hr />
  <blockquote>
  <ul>
  <li>Alignments with basecall quality values greater than 70 are rejected (these are not supported on the assumption that this indicates an offset error)</li>
  </ul>
  </blockquote>
  <p>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).</p>
  <hr />
  <hr />
  <p>Ultimately, I don't think you really need to deal with any of these with some command pipeline if you choose your inputs wisely.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-03-30T16:35:00.650Z
Answers.Text: 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.
Answers.UserDisplayName: Indira
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-03-30T23:52:48.360Z
Answers.Text: @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.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-03-31T00:01:26.210Z
Answers.Text: @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.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-11T13:29:38.117Z
Answers.Text: What do you mean start upstream? I am only provided the bams?
Answers.UserDisplayName: Indira
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-11T14:01:53.027Z
Answers.Text: I tried to edit the script as well.
Answers.UserDisplayName: Indira
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-11T14:05:02.043Z
Answers.Text: @Indira I mean if you can't start with pairs of FASTQ files, documenting the input requirements will go a long way in reducing the number of issues that arise. If you only have BAM files, I would either _test_ to see if the BAM contains suitable alignments _or_ just try to run Manta and include a useful error message in the logs when it fails. People who aren't familiar with your code and don't choose their inputs carefully should still be able to read the error message(s). A notice about the input not being suitable should suffice, and a link to a wiki page could be helpful here.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-11T14:06:59.870Z
Answers.Text: Regrettably the are not the best at reading error message(s)/will process everything which is why I need to create this script.
Answers.UserDisplayName: Indira
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-11T14:26:00.467Z
Answers.Text: @Indira But won't they still get an error when they try to run Manta using a BAM without any read/alignments in it?
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-11T14:29:24.983Z
Answers.Text: What do you mean?
Answers.UserDisplayName: Indira
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-11T14:36:55.383Z
Answers.Text: @Indira Well, don't you want Manta to fail when you feed it an empty BAM? If the input BAM contains ONT or PacBio long reads for example, your code will remove all of these.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-11T15:25:21.963Z
Answers.Text: Considering that there will be batches of samples used I don't want it to fail but to create an out file saying Sample 1 (out of n samples) could not be processed by the main manta script (just not for this limiting factors script which is needed to remove any limiting factors).
Answers.UserDisplayName: Indira
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-12T23:24:54.387Z
Answers.Text: @Indira So I think what you want is for _Manta_ to fail, but not your script? In this case, I would skip the pre-processing and just _try_ the command. For example: `{ >&2 echo "Running Manta for sample x"; ./runWorkflow.py ...; } || { >&2 echo "Could not process sample x"; do_cleanup(); }`
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-13T00:46:49.253Z
Answers.Text: Thank you Steve I will try and use that :)
Answers.UserDisplayName: Indira
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-03-30T14:37:32.573Z
Answers.LastActivityDate: 2023-03-30T14:37:32.573Z
Answers.UserDisplayName: Steve
Body: <p>I am trying to figure out a way to remove the limiting factors that Manta says it cannot handle. To quote from the <a href="https://github.com/Illumina/manta/blob/master/docs/userGuide/README.md" rel="nofollow noreferrer">Manta page</a>:</p>
  <blockquote>
  <p>The following limitations exist on the input BAM or CRAM files provided to Manta:</p>
  <ul>
  <li>Alignments cannot have an unknown read sequence (SEQ=&quot;*&quot;)</li>
  <li>Alignments cannot contain the &quot;=&quot; character in the SEQ field.</li>
  <li>Alignments cannot use the sequence match/mismatch (&quot;=&quot;/&quot;X&quot;) CIGAR notation
  RG (read group) tags in the alignment records are ignored -- each file will be
  treated as representing one sample.</li>
  <li>Alignments with basecall quality values greater than 70 are rejected (these are not
  supported on the assumption that this indicates an offset error)</li>
  </ul>
  </blockquote>
  <p>Would using a <code>find</code> statement in this way and then making <code>.bai</code> files be adequate to remove the limiting factors that Manta has?:</p>
  <pre><code>find <span class="math-container">${input_dir} -name "*.bam" -type f -exec sh -c 'output_file="$</span>{0%.bam}_filtered.bam&quot;; temp_file=<span class="math-container">$(mktemp); samtools view -h -q 70 "$</span>0&quot; | sed &quot;/^\(@\|[^[:blank:]]\+[[:blank:]]\+4\)/!s/\t[=*]/\t/g&quot; | samtools view -Sb - &gt; &quot;<span class="math-container">$temp_file"; mv "$</span>temp_file&quot; &quot;<span class="math-container">$output_file"; samtools index "$</span>output_file&quot;; rm &quot;$0&quot;' {} \;
  # Make the .bai files:
  
  for file in <span class="math-container">${input_dir}/*.bam; do
      samtools index "$</span>{file}&quot; &quot;${file%.*}.bai&quot;
  done
  </code></pre>
  <p>However, I get the error (repeatedly for all my bam files):</p>
  <blockquote>
  <p>[main_samview] fail to read the header from &quot;-&quot;.
  samtools index: &quot;/path/to/MyBatch1And2/SampleNumber_filtered_filtered_filtered.bam&quot; is in a format that cannot be usefully indexed</p>
  </blockquote>
  <p>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?</p>
  
Tags: <sam><bash>
Title: Dealing With Manta Limitations
ThreadText: # Dealing With Manta Limitations
  
  (Asked by: username_0 on 2023-03-29)
  
  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](https://github.com/Illumina/manta/blob/master/docs/userGuide/README.md):
  
  > 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?
  
  **Tags:** sam, bash
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2023-03-30. Score: 2
  
  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?](https://bioinformatics.stackexchange.com/questions/404/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.
  
  > **Comments:**
  > **username_0**: 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.  
  > **username_1**: @username_0 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.  
  > **username_1**: @username_0 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.  
  > **username_0**: What do you mean start upstream? I am only provided the bams?  
  > **username_0**: I tried to edit the script as well.  
  > **username_1**: @username_0 I mean if you can't start with pairs of FASTQ files, documenting the input requirements will go a long way in reducing the number of issues that arise. If you only have BAM files, I would either _test_ to see if the BAM contains suitable alignments _or_ just try to run Manta and include a useful error message in the logs when it fails. People who aren't familiar with your code and don't choose their inputs carefully should still be able to read the error message(s). A notice about the input not being suitable should suffice, and a link to a wiki page could be helpful here.  
  > **username_0**: Regrettably the are not the best at reading error message(s)/will process everything which is why I need to create this script.  
  > **username_1**: @username_0 But won't they still get an error when they try to run Manta using a BAM without any read/alignments in it?  
  > **username_0**: What do you mean?  
  > **username_1**: @username_0 Well, don't you want Manta to fail when you feed it an empty BAM? If the input BAM contains ONT or PacBio long reads for example, your code will remove all of these.  
  > **username_0**: Considering that there will be batches of samples used I don't want it to fail but to create an out file saying Sample 1 (out of n samples) could not be processed by the main manta script (just not for this limiting factors script which is needed to remove any limiting factors).  
  > **username_1**: @username_0 So I think what you want is for _Manta_ to fail, but not your script? In this case, I would skip the pre-processing and just _try_ the command. For example: `{ >&2 echo "Running Manta for sample x"; ./runWorkflow.py ...; } || { >&2 echo "Could not process sample x"; do_cleanup(); }`  
  > **username_0**: Thank you username_1 I will try and use that :)  
  
  ---
  
  