Entering edit mode
6.2 years ago
NicoN64
▴
30
Hello,
I was wondering if there is a method using perl or python to remove columns from MSA that have only gaps - and/or N and/or ? (everything except a nucleotide). It will be to run it on several MSA.
seq1
ATCGNN-??ATCG
seq2
ATCGNN--ATCGCG
seq3
ATCG-?NCGAAAAA
(Remove columns 5,6,7)
I know that there is some tool like TrimAl who has the option to remove gappy positions but I dont know if I can adapt it to look at positions with only characters '-' 'N' '?'. I want to keep ?,N,- in others part of the alignments, so I don't think that the command tr can help?
Thank you
did you try degap from scikit-bio python library? : http://scikit-bio.org/docs/0.2.3/generated/skbio.alignment.Alignment.degap.html ?
Your problem is sophisticated enough to warrant a custom script. You can use a couple of python dicts for this, one dict to hold ids-sequences as key value pairs and the other to hold position-unique_bases_at_position as key value pairs. The second dict is created using the first. By filtering the second dict, you can pick positions to display and write out these positions alone.
A word of caution with the above approach, you'll run into problems with long sequences or when you have a ton of sequences in the MSA, not to mention cases of unequal length sequences that you'll have to handle.
Also, please use the formatting bar (especially the
code
option) to present your post better. I've done it for you this time.You can have a look at this. Alternatively, you may wish to try UGENE for removing gaps. For other characters, it might be possible to convert them into '-' first?
You can use
sed
command if you just want to remove these characterssed 's/[N|?|-]//g' filename
How would
sed
account for the across-all-sequences ("columns") part?sed
can only process one line at a time.