1
$\begingroup$

I am trying to extract a specific sequence from a multifasta file, from each sequence in the aligned file.

The sequences look like this, and there are 32 sequences within the multiFASTA:

>NC_045512.2_Severe_acute_respiratory_syndrome_coronavirus_2_isolate_Wuhan-Hu-1_complete_genome 
ATT AAA GGT TTA TAC CTT CCC AGG TAA CAA ACC AAC CAA CTT TCG ATC TCT TGT AGA TCT GTT CTC TAA ACG AAC TTT 
AAA ATC TGT GTG GCT GTC ACT CGG CTG CAT GCT TAG TGC ACT CAC GCA GTA TAA TTA ATA ACT AAT TAC TGT CGT TGA 
CAG GAC ACG AGT AAC TCG TCT ATC TTC TGC AGG CTG CTT ACG GTT TCG TCC GTG TTG CAG CCG ATC ATC AGC ACA TCT 
AGG TTT CGT CCG GGT GTG ACC GAA AGG TAA GAT GGA GAG CCT TGT CCC TGG TTT CAA CGA GAA AAC ACA CGT CCA ACT 

Using the following code I can pull the sequence, using coordinates, from a single FASTA file, but how would one iterate this to pull the section (by coordinates) from each sequence within the multiFASTA?

from Bio import SeqIO
from Bio.SeqIO.FastaIO import FastaWriter
import os

os.chdir("C:/Users/")

record = SeqIO.read("Whole genome.fas", "fasta")

with open("Spike protein.fasta", "wt") as out:
    writer = FastaWriter(out, wrap = [])
    writer.write_header()
    writer.write_record(record[21563:25384])
    writer.write_footer()

Thanks in advance!

$\endgroup$
2
$\begingroup$

You can use SeqIO for the writing in addition to reading the sequences. This has the advantage that it will always write the correct fasta header, which is not always just the .id if there is whitespace, e.g. (untested)

from Bio import SeqIO

with open("Spike protein.fasta", "w") as out:
    for seq_record in SeqIO.parse("Whole genome.fas", "fasta"):
        seq_record.seq = seq_record.seq[21563:25384]
        SeqIO.write(seq_record, out, "fasta")

BTW don't put spaces in your filenames, it can lead to issues

$\endgroup$
1
$\begingroup$

I developed an answer as follows:

from Bio import SeqIO
from Bio.SeqIO.FastaIO import FastaWriter
import os

os.chdir("C:/Users/")

with open("Spike protein.fasta", "wt") as out:
    for seq_record in SeqIO.parse("Whole genome.fas", "fasta"):
        out.write(str(seq_record.id) + "\n")
        out.write(str(seq_record.seq[21563:25384]) + "\n")
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.