2
$\begingroup$

I want to fix a chromosome of around 50-100Mb (chr.fasta) for a region that is missing a few kilobases of data. I have the missing sequence with a window upstream and downstream that matches perfectly to the reference (contig.fasta).

What's the best way to produce a new chromosome, fixed.fasta that contains the new data from contig.fasta?

See depiction below:

chr.fasta    |---------------------------------------------------   -------------|
contig.fasta                                                   |-===-|

Final intended result:

fixed.fasta  |---------------------------------------------------===-------------|

Any ideas?

$\endgroup$
1
$\begingroup$

I have not received any answer to the question, so I created a small tool to fix the file myself.

I am assuming as per the question that there are perfect matches upstream and downstream to the new sequence, and that anything in between can be deleted/updated by the new sequence.

Given the first_n and last_n bases of the diagram in the question, and the entire patch sequence, we can create a single-line string of the reference with seqtk seq -l9999999 - and then edit it with the line below (in Perl regex):

$chromosome_str =~ s/$first_n.+$last_n/$patch/i

This creates the new string that can be used as the new fasta file.

$\endgroup$
0
$\begingroup$

You can give a try with megamerger from the emboss suite if you have (almost) perfect sequence match.

In case the gap within your chr.fasta complicates the sequence alignment and hence the resulting contig, you can generate a "tiling path" with: upstream of the gap from chr.fasta + contig.fasta + downstream of the gap from chr.fasta (there should be sequence overlaps of course)

In a (more realistic) scenario where there is no perfect sequence match between the overlapping regions, you can fall back to emboss merger, however, it is not designed to work on large sequences.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.