Answers.Body: <p>I'd suggest looking at similar tools that already exist, e.g. <a href="https://github.com/ponnhide/flashpy" rel="nofollow noreferrer">here is one that I hadn't heard of that lives on github</a>. Just glancing at the code, it seems to use a scoring heuristic to decide how to derive a merge point, so it's using that to get around the absence of unambiguous points, probably at some rate of error.</p>
  <p><a href="https://www.biostars.org/p/284637/" rel="nofollow noreferrer">BioPython has some functions</a> that will create consensus sequences for you as well based on a pairwise alignment. (<code>Bio.Align.AlignInfo.SummaryInfo.dumb_consensus</code>, linked in that post.)</p>
  <p>As a command line tool <a href="https://pandaseq-tutorial.readthedocs.io/en/latest/tutorial/" rel="nofollow noreferrer">pandaseq</a> does the same thing, I'd suggest looking at how they handle edge cases.</p>
  <p>There are similar utilities in <a href="http://qiime.org/scripts/join_paired_ends.html" rel="nofollow noreferrer">QIIME</a>, <a href="https://dash.harvard.edu/bitstream/handle/1/42661733/MS_9.pdf?sequence=10" rel="nofollow noreferrer">NGMerge</a>, <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5657622/" rel="nofollow noreferrer">BBMerge</a>, [illumina-utils]- this is a pretty well-trodden path, I think all of those will be illuminating to look at in terms of reusing or writing original code.</p>
  <p>I'm guessing it will be easier to merge FASTQs and then align, rather than align and then merge. (Note that these tools don't seem to operate on alignments, because that is trickier).</p>
  <p>Alternately, if your real goal is to create a consensus genome sequence from your reference based on your alignments, there are a large number of tools that do that, notably <code>samtools</code> (see e.g. <a href="https://www.biostars.org/p/419927/" rel="nofollow noreferrer">https://www.biostars.org/p/419927/</a>). No new code needed, you just input your BAM file and reference sequence.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-09-29T22:54:41.397Z
Answers.LastActivityDate: 2023-09-29T22:54:41.397Z
Answers.UserDisplayName: Maximilian Press
Body: <p>I'm writing a function (python) that merges overlapping paired reads of a fragment, and a corresponding reference sequence. There's quite a lot of edge cases (e.g. when the reads disagree on indels) so I wonder if there's some standard library that does that?</p>
  <p>Context - the input is pysam AlignedSegment DNA reads. The output is fragment, reference - both strings corresponding to the merged fragment and its corresponding reference.</p>
  <p>e.g.</p>
  <pre><code>read1:           ACGTACGTACGT
  read2:                CGTACGTACGT
  expected output: ACGTACGTACGTACGT
  
  read1:           ACGTANGTACGT
  read2:               ACGTACGTACGT
  expected output: ACGTACGTACGTACGT
  </code></pre>
  <p>this is just an illustration. The idea is to merge the overlapping DNA read sequences into the fragment and its corresponding reference using the information in the AlignedSegment (i.e. query_sequence, reference_start, query_qualities, cigartuples, and get_reference_sequence())</p>
  <p>Thanks.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-09-29T10:13:01.120Z
Comments.Text: Thanks. But why is this not an assembler? This is pretty much exactly what an assembler does, it assembles reads into contigs. Or, if you already have the reference genome, what a read aligner does.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-09-29T10:22:40.370Z
Comments.Text: well, I mean it's not an aligner in the sense that's not aligning reads to a reference but to each other.
Comments.UserDisplayName: user1467422
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-09-29T10:25:05.733Z
Comments.Text: OK, so it's an assembler. But then what do you mean by the word "reference" in your question? Are you referring to the consensus sequence you can build by assembling the reads, the inferred genomic source sequence? Also, do you absolutely have to write this from scratch and in a scripting language? Can't you use one of the many [existing genome assemblers](https://en.wikipedia.org/wiki/De_novo_sequence_assemblers#Commonly_used_programs)?
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-09-29T10:33:55.740Z
Comments.Text: I believe the OP wants to work at the "read-pair level" to obtain contigs of mate pairs and not all of the available reads. If that is the case, I would check [`make.contigs`](https://mothur.org/wiki/make.contigs/) from the `mother` suite.
Comments.UserDisplayName: haci
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-06T06:32:08.620Z
Comments.Text: I would use existing high-performance software for this, for example BBmerge from BBmap suite, or [FLASH](http://ccb.jhu.edu/software/FLASH/), probably many others exist as well. I see no point in re-writing new software when it already exists, unless this is a learning exercise, or you think you can do better algorithm-wise, but then, with respect, you would probably not be asking here.
Comments.UserDisplayName: ATpoint
Tags: <pysam>
Title: Merge fragment reads
ThreadText: # Merge fragment reads
  
  (Asked by: username_0 on 2023-09-29)
  
  I'm writing a function (python) that merges overlapping paired reads of a fragment, and a corresponding reference sequence. There's quite a lot of edge cases (e.g. when the reads disagree on indels) so I wonder if there's some standard library that does that?
  
  Context - the input is pysam AlignedSegment DNA reads. The output is fragment, reference - both strings corresponding to the merged fragment and its corresponding reference.
  
  e.g.
  
  ```
  read1:           ACGTACGTACGT
  read2:                CGTACGTACGT
  expected output: ACGTACGTACGTACGT
  
  read1:           ACGTANGTACGT
  read2:               ACGTACGTACGT
  expected output: ACGTACGTACGTACGT
  
  ```
  
  this is just an illustration. The idea is to merge the overlapping DNA read sequences into the fragment and its corresponding reference using the information in the AlignedSegment (i.e. query\_sequence, reference\_start, query\_qualities, cigartuples, and get\_reference\_sequence())
  
  Thanks.
  
  **Tags:** pysam
  
  > **Comments:**
  > **username_1**: Thanks. But why is this not an assembler? This is pretty much exactly what an assembler does, it assembles reads into contigs. Or, if you already have the reference genome, what a read aligner does.  
  > **username_0**: well, I mean it's not an aligner in the sense that's not aligning reads to a reference but to each other.  
  > **username_1**: OK, so it's an assembler. But then what do you mean by the word "reference" in your question? Are you referring to the consensus sequence you can build by assembling the reads, the inferred genomic source sequence? Also, do you absolutely have to write this from scratch and in a scripting language? Can't you use one of the many [existing genome assemblers](https://en.wikipedia.org/wiki/De_novo_sequence_assemblers#Commonly_used_programs)?  
  > **username_2**: I believe the OP wants to work at the "read-pair level" to obtain contigs of mate pairs and not all of the available reads. If that is the case, I would check [`make.contigs`](https://mothur.org/wiki/make.contigs/) from the `mother` suite.  
  > **username_3**: I would use existing high-performance software for this, for example BBmerge from BBmap suite, or [FLASH](http://ccb.jhu.edu/software/FLASH/), probably many others exist as well. I see no point in re-writing new software when it already exists, unless this is a learning exercise, or you think you can do better algorithm-wise, but then, with respect, you would probably not be asking here.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_4 on 2023-09-29. Score: 1
  
  I'd suggest looking at similar tools that already exist, e.g. [here is one that I hadn't heard of that lives on github](https://github.com/ponnhide/flashpy). Just glancing at the code, it seems to use a scoring heuristic to decide how to derive a merge point, so it's using that to get around the absence of unambiguous points, probably at some rate of error.
  
  [BioPython has some functions](https://www.biostars.org/p/284637/) that will create consensus sequences for you as well based on a pairwise alignment. (`Bio.Align.AlignInfo.SummaryInfo.dumb_consensus`, linked in that post.)
  
  As a command line tool [pandaseq](https://pandaseq-tutorial.readthedocs.io/en/latest/tutorial/) does the same thing, I'd suggest looking at how they handle edge cases.
  
  There are similar utilities in [QIIME](http://qiime.org/scripts/join_paired_ends.html), [NGMerge](https://dash.harvard.edu/bitstream/handle/1/42661733/MS_9.pdf?sequence=10), [BBMerge](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5657622/), [illumina-utils]- this is a pretty well-trodden path, I think all of those will be illuminating to look at in terms of reusing or writing original code.
  
  I'm guessing it will be easier to merge FASTQs and then align, rather than align and then merge. (Note that these tools don't seem to operate on alignments, because that is trickier).
  
  Alternately, if your real goal is to create a consensus genome sequence from your reference based on your alignments, there are a large number of tools that do that, notably `samtools` (see e.g. <https://www.biostars.org/p/419927/>). No new code needed, you just input your BAM file and reference sequence.
  
  
  ---
  
  