This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Problem trimming SAM file to use with assembler

Hi everyone,

I'm trying to perform de novo assembly of the region between two known non-repeating sequences (two "anchors"). I have a fastq files from ONT long read sequencing of the entire genome and I'm trying to implement an iterative workflow with the following steps:

  1. Align reads to one anchor (the "reference anchor")
  2. Trim the side that is not pointing to the other anchor
  3. Assemble a de novo sequence, with the objective of extending the anchor
  4. Check if the reference anchor was extended or if the other anchor appeared in the new de novo sequence
  5. If the assembly was extended but the other anchor did not appear, go back to step 1 with the extended sequence as an anchor
  6. If the assembly could not be extended or the other anchor was found, exit the loop

I'm having trouble with step 2, where I can't seem to make my trimming scripts work.

Is there an already existing functionality to trim a SAM alignment and generate its corresponding fastq file (with trimmed reads)? I have looked for trimming scripts but all I have found is for trimming barcodes or to select reads that align to a specific region (all of my reads align to the small anchor), but nothing about trimming reads, especially not a one-direction trim (all of the programs seem to require specific regions with start and end positions, and I want to extend as much as possible in one direction while trimming in the other).

Alternatively, is there a simpler way to make this work that I'm not seeing? I don't trust re-aligning trimmed reads nor generating the new de novo sequence and then trimming it (what if the alignment changes with trimmed reads? what if the de novo sequence contains more copies of the anchor sequence? this is being done on a repetitive region, after all), but maybe those can be done safely in some way I don't know of.

Lastly, I understand that step 6 may need some other exit conditions (e.g. to avoid looping infinitely if stuck on a repetitive region), but I want to focus on solving the problem at hand first.

Thanks in advance.

fastq sam assembly

The simplest approach to trimming would be to write your own BAM file processing with say pysam or in a different programming language.

Specifically in pysam the get_aligned_pairs() method for each alignment returns a base-by-base alignment, making it straightforward to reconstruct the aligned regions.

Plus, you can run many other steps along the ay

may be but - sorry - I don't understand the question.

This part:

existing functionality to trim a SAM alignment and generate its corresponding fastq file (with trimmed reads)?

0 answers

No answers yet.

Log in to answer this question.