I adapted the script and it doesn't work. Sometimes, if one of the aligned sequences has a relatively big gap in the middle, everything downstream is eliminated.
I have used muscle to align a newly assembled genome to a reference genome. The seqs.afa file looks like this:
>reference
--------------------------ACTGAC
ACTGACTGACTGACTGACTGACTGACTGACTG
... Lots of bases here .........
ACTGACTGACTGACTGACTGACTGACTG----
-----------------
>my_assembly
AAAAAAAAAAAAAAAAAAAAAAAAAAACTGAC
ACTGACTGACTGACTGACTGACTGACTGACTG
... Lots of bases here .........
ACTGACTGACTGACTGACTGACTGACTGAAAA
AAAAAAAAAAAAAAAAA
As you can see, my program has a tendency to leave dangling bases downstream and upstream the reference genome. and I need to get rid of them in post processing. Is there a program or simple Python script I can use to trim the bases that overhang from the reference. How to do it efficiently? (I have a considerable amount of data sets)
3 answers
Is your input data in the same format as OP's of that post? No linebreaks in sequences?
trimal (http://trimal.cgenomics.org/publications) can trim multiple sequence alignment results
Jalview (www.jalview.org) is a good visualisation workbench for multiple sequence alignments. Along with several useful editing features (check out the youtube video), it has pad gaps feature that can be toggled on and off. When selected, the alignment will be kept at a minimal width (so there are no empty columns before or after the first or last aligned residue) and all sequences will be padded with gap characters before and after their terminating residues. The pad gaps feature is demonstrated at 3.20min in this video. The sequence can then be exported in a variety of file formats.
Log in to answer this question.