Question: Software for Unaligning Sequences
0
gravatar for printingdots
20 months ago by
printingdots0 wrote:

The dataset I have consists of aligned sequences in fasta format. However, I am developing a certain multiple sequence alignment method, for which I need to run my program on unaligned sequences. Is there any software for converting from aligned fasta sequences to unaligned sequences in fasta format?

alignment • 749 views
ADD COMMENTlink modified 20 months ago by Joe17k • written 20 months ago by printingdots0

please post an example of what you mean "aligned sequences". FASTA as a format is not an alignment format, hence it is not clear what you are looking for

ADD REPLYlink written 20 months ago by Istvan Albert ♦♦ 84k

This is part of what my aligned fasta file contains:

>>
> SDEA
ATT-------AG-------------------------GAGGGT-----------------
-----ACCATTGT--------------ACCGGTTGGT-T---------------------
---------ATCAAA---------------------TCG------CCG------------
-GA-----------------------------------------AACCTA----T-----
------------------GTGC----------CCG-------CAGCG-------------
----TGATTTC----A-----------------------CGATGGT------C-------
---------------------CC---T--------------------ATGTGGGC-----
--CA---------------------------A----TCCCGA---------------CAA
AGGCCGAAAA------A-ATT-GAC-G---AC---------TATGT--------------
---------------GAG-----GTGC--------------------GTC-----AA---
-----------------C------------TGT------------------A--------
-----TG------CCTCAAAGTAGCC-----A---------------TGAGTT-------
------------------------CC----------TT----CGCAGGC-----------
--------------------------TTCCCCCCATAGGGGTTT--GTA--TTAC-----
-----TGAA------GATCT----CTAT----AC-------AA-----------------
--GGGAAGTTTTC-----------------------------------------------
G--GT-------------AAGG---CG----------TC------------------CTG
C------------------------------------------G-A-T----------G-
--G------G--------------------CG--------------------------CC
GATTGA--------TGCT------------------------------------------
-------CATG---------------------------------T---------------
---T--------------------------------------AGCAGCGG---------A
--------AATTG--TCTCCG-----------------CGC--------------GTC--
---------A------GAGTGAAAACTCCTC------------CA---------------
------CTATT-----------------------------AGCCCCA-------------
--------------------------AGGCGGT-------AG---A--------------
------TGT-------------------------C-----GATG-------TGG------
---------------------------------------------------------GT-
------------TAAAT-------------------------------------------
--------------TGA---C-CG-----AA-----------------------------
----------GC
ADD REPLYlink modified 20 months ago by Istvan Albert ♦♦ 84k • written 20 months ago by printingdots0
1

Run it through awk and delete - gap characters from the sequence string.

ADD REPLYlink modified 20 months ago • written 20 months ago by Alex Reynolds30k
2
gravatar for Alex Reynolds
20 months ago by
Alex Reynolds30k
Seattle, WA USA
Alex Reynolds30k wrote:
$ awk '{ if ($0 !~ /^>/) { gsub(/-/, "", $0); } print $0; }' in.fa > out.fa
ADD COMMENTlink modified 20 months ago • written 20 months ago by Alex Reynolds30k

Error:

$ awk '{ if ($0 !~ /^>/) { gsub(//, '-', $0); } print $0; }' test.fa
awk: cmd. line:1: { if ($0 !~ /^>/) { gsub(//, -, $0); } print $0; }
awk: cmd. line:1:                               ^ syntax error

Modified:

 $ awk '{ if ($0 !~ /^>/) { gsub("-", "",$0); } print $0; }' test.fa 
>seq-1-1
atgc
>seq2
gc
>seq_3
gatgatg

awk version:

 $ awk --version
    GNU Awk 4.1.4, API: 1.1 (GNU MPFR 4.0.1, GNU MP 6.1.2)

input example fasta:

$ cat test.fa 
>seq-1-1
a-t-g----c
>seq2
--g--c---
>seq_3
-gatga-tg
ADD REPLYlink modified 20 months ago • written 20 months ago by cpad011213k

Thanks, fixed the order of substitution.

ADD REPLYlink written 20 months ago by Alex Reynolds30k
4
gravatar for cpad0112
20 months ago by
cpad011213k
India
cpad011213k wrote:
$ cat test.fa 
>seq-1-1
a-t-g----c
>seq2
--g--c---
>seq_3
-gatga-tg

$ sed '/>/! s/-//g' test.fa 
>seq-1-1
atgc
>seq2
gc
>seq_3
gatgatg
ADD COMMENTlink written 20 months ago by cpad011213k
2
gravatar for Joe
20 months ago by
Joe17k
United Kingdom
Joe17k wrote:

In practice I’d use a lot of the solutions above, but a ‘safe’ option would be to use a dedicated parser like BioPython. For instance, Istvan’s solution (I think) would edit your headers if they had - in:

from Bio import AlignIO
import sys

aln = AlignIO.read(sys.argv[1], 'fasta')
for s in aln:
    print('>{}\n{}'.formats.id, str(s.seq).replace('-',''))
ADD COMMENTlink modified 20 months ago • written 20 months ago by Joe17k
1

To be further safe, onr can use either degap from scikitbio or ungap from Bio.seq module instead of python for replacing the gaps. In addition, there is an extra ) at the end of the code.

from Bio import SeqIO
import sys

for record in SeqIO.parse(sys.argv[1], "fasta"):
    print(">", record.id, "\n", record.seq.ungap("-"), sep="")

output:

$ python3 test.py test.fa 
>seq-1-1
atgc
>seq2
gcatgc
>seq_3
gatgatgatgcctg

input:

$ cat test.fa
>seq-1-1
a-t-g----c
>seq2
--g--c---
atgc
>seq_3
-gatga-tg
----atgc---
-ctg
ADD REPLYlink modified 20 months ago • written 20 months ago by cpad011213k

Ah that's good. I was sure I'd seen a method for gap removal somewhere before but couldn't find it again earlier.

Incidentally, casting to str() and then running a replace is exactly what ungap does internally ;) but always nice to use the designed methods.

From the source:

1245          if hasattr(self.alphabet, "gap_char"): 
1246              if not gap: 
1247                  gap = self.alphabet.gap_char 
1248              elif gap != self.alphabet.gap_char: 
1249                  raise ValueError( 
1250                      "Gap {0!r} does not match {1!r} from alphabet".format( 
1251                          gap, self.alphabet.gap_char)) 
1252              alpha = Alphabet._ungap(self.alphabet) 
1253          elif not gap: 
1254              raise ValueError("Gap character not given and not defined in " 
1255                               "alphabet") 
1256          else: 
1257              alpha = self.alphabet  # modify! 
1258          if len(gap) != 1 or not isinstance(gap, str): 
1259              raise ValueError("Unexpected gap character, {0!r}".format(gap)) 
1260          return Seq(str(self).replace(gap, ""), alpha)
ADD REPLYlink modified 20 months ago • written 20 months ago by Joe17k
1
gravatar for Istvan Albert
20 months ago by
Istvan Albert ♦♦ 84k
University Park, USA
Istvan Albert ♦♦ 84k wrote:

Delete the - characters then reformat to "nice" fasta with:

cat seq.fa | tr -d -  | seqtk seq -A
ADD COMMENTlink written 20 months ago by Istvan Albert ♦♦ 84k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 773 users visited in the last hour