Question: change specific bases in a fasta file
0
gravatar for samuel.lipworth
2.5 years ago by
University of Oxford
samuel.lipworth30 wrote:

Hi,

I have a list of positions I would like to change in a fasta file with the base I would like to change to tab separated as follows:

1 T 2 C 10 T 50 G

etc

so that eg: ATGTGTC...

becomes TCGTGTC...

There must be a way to do this but I can't do it - have tried to code it in biopython but can't manage.. any ideas?

biopython fasta • 2.3k views
ADD COMMENTlink modified 2.5 years ago • written 2.5 years ago by samuel.lipworth30

post what have you tried in biopython.

ADD REPLYlink written 2.5 years ago by geek_y10k

see below, thanks...

ADD REPLYlink written 2.5 years ago by samuel.lipworth30
import csv
from Bio import SeqIO
from Bio.Alphabet import generic_dna
from Bio.Seq import MutableSeq
from sys import argv
import sys
script, reference, changes, out = argv
changes_=open(changes, 'rb')
change = csv.reader(changes_, delimiter='\t')
for line in change:
        pos = line[0]
        base = line[1]
        with open(out, 'w') as output:
                with open(reference, 'rb') as fasta_file:
                        for seq_record in SeqIO.parse(fasta_file, "fasta"):
                                seq_record[pos]=base
                                SeqIO.write(fasta_file, output, 'fasta')

I know this is not quite right but not sure how to make it work

ADD REPLYlink written 2.5 years ago by samuel.lipworth30

you don't need to import so many modules. Also, you need to subtract 1 from your positions, because python starts at 0.

ADD REPLYlink modified 2.5 years ago • written 2.5 years ago by st.ph.n2.5k

This is the solution you need but you need to work a little bit to prepare input files

Introducing Known Mutations (From A Vcf) Into A Fasta File

ADD REPLYlink written 2.5 years ago by venu6.3k
3
gravatar for st.ph.n
2.5 years ago by
st.ph.n2.5k
Philadelphia, PA
st.ph.n2.5k wrote:

You don't specify how many sequences are in your file, or how to correlate the positions/changes to each sequence. I suggest adding header information, if there is more than one sequence to your position table. The code below assumes multiple sequences, but will work with only one as well.

First, change this 1 T 2 C 10 T 50 G

to this (tab-delimited, pos.txt):

head1 1 T
head1 2 C 
head1 10 T 
head1 50 G
....
head2 ....

Linearize fasta (a preference I have. You can still read in the sequences with Biopython).

#!/usr/bin/env python
from collections import defaultdict

with open('pos.txt', 'r') as f:
        pos = defaultdict(list)
        for line in f:
                pos[line.strip().split('\t')[0]].append((int(line.strip().split('\t')[1]), line.strip().split('\t')[2]))

with open('input.fasta', 'r') as fasta:
        with open('input_corr.fasta', 'w') as out:
                for line in fasta:
                        if line.startswith(">"):
                                h = line.strip().split('>')[1]
                                s = list(next(fasta).strip())
                                if h in pos:
                                        for n in pos[h]:
                                                s[n[0]-1] = n[1]
                                        out.write('>' + h + '\n' + ''.join(s))
ADD COMMENTlink modified 2.5 years ago • written 2.5 years ago by st.ph.n2.5k
2
gravatar for Matt Shirley
2.5 years ago by
Matt Shirley9.2k
Cambridge, MA
Matt Shirley9.2k wrote:

You can use the MutableFastaRecord in pyfaidx for this exact purpose.

The mutable Fasta instance will modify your file in-place (commits changes directly back to disk) so be careful which file you're editing.

ADD COMMENTlink written 2.5 years ago by Matt Shirley9.2k
0
gravatar for samuel.lipworth
2.5 years ago by
University of Oxford
samuel.lipworth30 wrote:

Thanks for the help!

I also managed to write my own version which seems to work to!

https://github.com/samlipworth/base_changer

ADD COMMENTlink written 2.5 years ago by samuel.lipworth30

and for the same position in different chromossomes?

ADD REPLYlink written 6 weeks ago by mgdias.jose10
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: 2043 users visited in the last hour