This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Any ideas on converting IUPAC ambiguity codes in a fasta file

Hi, I have fasta files that contain ambiguity codes. And I would like to convert the sequences so that it's just ACGT and create a new entry for each possibility. i.e. Y = C or T

>sampleID
AGCTAGYACG

into 

>sampleID.1
AGCTAGCACG
>sampleID.2
AGCTAGTACG

I was thinking about writing a Perl script that does that using the transliterate function (tr). Before I embark on this, has anyone already written or used some code/software that does this?

fasta iupac regex perl bash
• 6,044 views • link
updated by schagas • schagas • written by SaltedPork SaltedPork

3 answers

see all posible sequences from consensus

$ gcc biostars282490.c
$ ./a.out AGCTAGYACG
AGCTAGCACG
AGCTAGTACG
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

This python code should work if you FASTA is single-line. If not you can linearize it. 

#!/usr/bin/env python
import sys
with open(sys.argv[1], 'r') as f:
    with open(sys.argv[2], 'w') as out:
        for line in f:
            if line.startswith(">"):
                header = line.strip()
                seq = next(f).strip()
                if 'Y' in seq:
                    out.write(header + '.1' + '\n' + 'C'.join(seq.split('Y')))
                    out.write(header + '.2' + '\n' + 'T'.join(seq.split('Y')))
                else:
                    out.write(header + '\n' + seq)

Save as convert.py, run as python convert.py input.fasta output.fasta

EDIT

IUPAC_ambiguity_codes.txt

M       A       C
R       A       G
W       A       T
S       C       G
Y       C       T
K       G       T
V       A       C       G
H       A       C       T
D       A       G       T
B       C       G       T
N       G       A       T       C
  
#!/usr/bin/env python
import sys
from collections import defaultdict

with open('IUPAC_ambiguity_codes.txt', 'r') as amb:
        codes = defaultdict(list)
        for line in amb:
                for line in f:
                        codes[line.strip().split('\t')[0]] = line.strip().split('\t')

with open(sys.argv[1], 'r') as f:
    with open(sys.argv[2], 'w') as out:
        for line in f:
            if line.startswith(">"):
                header = line.strip()
                seq = next(f).strip()
                for c in codes:
                        if c in seq:
                                for n in range(len(codes[c])):
                                        out.write(header + '.' + str(n+1) + '\n' + codes[c][n].join(seq.split(c)))

Untested, but should output a new numbered sequence header for each replacement required to make for each ambiguity code encountered.

• 0 views • link
written by st.ph.n st.ph.n

I am reasonably sure OP wants this to work for all valid IUPAC codes.

• 0 views • link
written by GenoMax GenoMax

Yep, that would be the idea, I suppose I could take this and have many if statements, but then how do I account for sequences which have multiple IUPAC codes in them?

• 0 views • link
written by SaltedPork SaltedPork

You sure can. If you want something that works now then use one of the stackoverflow or @Pierre's solution below.

• 0 views • link
written by GenoMax GenoMax

It appears as though you want a new sequence for each ambiguity code encountered. See edited post for reference. Otherwise if you want a new sequence and ALL codes changed at once, it will be slightly different.

• 0 views • link
written by st.ph.n st.ph.n

As my lab has the limitation to work only with the Sanger sequencer, I needed to get the two versions (ignoring the infinity of combinations/possibilities) from sequences with ambiguity code just to visualize the nucleotide diversity in heterozygous samples. I tested the code above and made some changes.

#!/usr/bin/env python

import subprocess
import sys
import os
from collections import defaultdict

def install(package):
    subprocess.check_call([sys.executable, "-m", "pip", "install", package])

def read_fasta(fasta):
    with open(fasta, 'r') as fast:
        headers, sequences = [], []
        for line in fast:
            if line.startswith('>'):
                head = line.replace('>','').strip()
                headers.append(head)
                sequences.append('')
            else :
                seq = line.strip()
                if len(seq) > 0:
                    sequences[-1] += seq
    return (headers, sequences, fast)


def write_fasta(headers, sequences, fasta):
    with open(fasta, 'w') as fast:
        for i in range(len(headers)):
            fast.write('>' + headers[i] + '\n' + sequences[i] + '\n')

install('biopython')
print (sys.path)
from Bio.SeqIO import FastaIO







with open('IUPAC_ambiguity_codes.txt', 'r') as amb:
        codes = defaultdict(list)
        for line in amb:
            codes[line.strip().split()[0]] = line.strip().split()

try:
    os.makedirs('./msa-wo', exist_ok=True)
    print("Directory ./msa-wo created successfully" )
except OSError as error:
    print("Directory ./msa-wo can not be created")

fasta = read_fasta("big.fasta")
new_fasta = write_fasta(fasta[0], fasta[1],'new_fasta_nobreaklines.fasta')

with open('new_fasta_nobreaklines.fasta', 'r') as f:
    with open("./msa-wo/msa_wo_amb.fasta", 'w') as out:
        for line in f:
            if line.startswith(">"):
                        header = line.strip()
                        seq = next(f).strip()
                        seq1 = seq
                        seq2 = seq
            for c in codes:
                        if c in seq1:
                            seq1 = codes[c][1].join(seq1.split(c))
                        if c in seq2:
                            seq2 = codes[c][2].join(seq2.split(c))
            out.write(header + '.1' + '\n' + seq1 + '\n\n')
            out.write(header + '.2' + '\n' + seq2 + '\n\n') 

print("Output: msa_wo_amb.fasta, check msa-wo folder at convert.py directory")

Input:

>1f
ACATTA

>2f
WCCCRA

Result:

>1f.1
ACATTA

>1f.2
ACATTA

>2f.1
ACCCAA

>2f.2
TCCCGA
• 0 views • link
written by schagas • schagas

Log in to answer this question.

More posts like this