<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>

CC BY-SA 4.0
2023-09-29T22:54:41.397Z
2023-09-29T22:54:41.397Z
Maximilian Press
<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>

CC BY-SA 4.0
2023-09-29T10:13:01.120Z
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.
terdon
CC BY-SA 4.0
2023-09-29T10:22:40.370Z
well, I mean it's not an aligner in the sense that's not aligning reads to a reference but to each other.
user1467422
CC BY-SA 4.0
2023-09-29T10:25:05.733Z
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)?
terdon
CC BY-SA 4.0
2023-09-29T10:33:55.740Z
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.
haci
CC BY-SA 4.0
2023-10-06T06:32:08.620Z
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.
ATpoint
<pysam>
Merge fragment reads
# 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.


---

