如何用Biopython拆分多序列FASTA文件为等长片段并修改序列头

Hey there! I get exactly what you're trying to do—splitting a multi-FASTA file into chunks based on total base pairs, not just number of records, while updating the headers with position info. The Biopython batch_iterator is a solid starting point, but it's built for grouping records, not bp counts. Let's adapt it (and add some extra logic) to fit your needs.

Core Approach #

First, let's outline the key steps we need to make this work:

• Calculate the total number of base pairs across all sequences in your input FASTA.
• Determine the target size for each chunk: most chunks use the ceiling of total_bp / split_num, while the last chunk takes whatever's left.
• Iterate through each sequence, slicing it into fragments to fill up each chunk, tracking start/end positions in the original contig.
• Write each completed chunk to a new FASTA file once it hits the target size.

Full Working Code #

Here's the commented code that implements exactly what you need:

from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

def split_fasta_by_bp(input_fasta, split_num):
    # Step 1: Calculate total base pairs across all records
    total_bp = 0
    records = list(SeqIO.parse(input_fasta, "fasta"))
    for rec in records:
        total_bp += len(rec.seq)
    
    # Step 2: Compute chunk size (ceiling division to get equal first chunks)
    chunk_size = (total_bp + split_num - 1) // split_num  # No need for math.ceil!
    
    # Initialize chunk tracking variables
    current_chunk_bp = 0
    current_chunk_records = []
    chunk_index = 1
    
    # Step 3: Process each record and build chunks
    for rec in records:
        rec_len = len(rec.seq)
        current_pos = 0
        
        while current_pos < rec_len:
            # How much can we take from this record to fill the current chunk?
            remaining_in_chunk = chunk_size - current_chunk_bp
            take = min(remaining_in_chunk, rec_len - current_pos)
            
            # Slice the sequence and create a new record with updated header
            seq_slice = rec.seq[current_pos : current_pos + take]
            new_id = f"{rec.id}:{current_pos}-{current_pos + take - 1}"
            new_rec = SeqRecord(seq_slice, id=new_id, description="")
            
            # Add to current chunk
            current_chunk_records.append(new_rec)
            current_chunk_bp += take
            current_pos += take
            
            # Write chunk if it's full, or if it's the final chunk
            if current_chunk_bp == chunk_size or (chunk_index == split_num and current_pos >= rec_len):
                output_filename = f"chunk{chunk_index}.fas"
                SeqIO.write(current_chunk_records, output_filename, "fasta")
                print(f"Wrote chunk {chunk_index} ({current_chunk_bp} bp) to {output_filename}")
                
                # Reset for next chunk
                current_chunk_bp = 0
                current_chunk_records = []
                chunk_index += 1

# Example usage (match your input and split count)
if __name__ == "__main__":
    split_fasta_by_bp("original.fas", split_num=4)

Key Details Explained #

• Total BP Calculation: We first read all records to get the total length—this ensures we split the entire dataset evenly.
• Chunk Size Trick: (total_bp + split_num - 1) // split_num is a clean way to compute ceiling division without importing math. For your example, 25 + 4 -1 = 28, 28//4 =7, which matches your desired 7bp chunks for the first 3 files.
• Header Updates: The new header uses 0-based start/end positions exactly as you specified, so contig2:0-1 refers to the first two bases of contig2.
• Cross-Contig Chunks: The code seamlessly handles chunks that span multiple contigs (like your chunk1, which includes the end of contig1 and start of contig2).

Test with Your Example #

If you run this with your original.fas and split_num=4, you'll get exactly the output you described:

• chunk1.fas: 7bp total (full contig1 + first 2 bases of contig2)
• chunk2.fas: 7bp (contig2 positions 2-8: GATAGTC)
• chunk3.fas:7bp (contig2's last base + first 6 bases of contig3)
• chunk4.fas:4bp (final 4 bases of contig3)

内容的提问来源于stack exchange，提问作者Gunther