Find length of a contig in one fasta, using the header of another fasta as query in python

user_7614836 0
Viewed 638

I'm trying to find a python solution to extract the length of a specific sequence within a fasta file using the full header of the sequence as the query. The full header is stored as a variable earlier in the pipeline (i.e. "CONTIG"). I would like to save the output of this script as a variable to then use later on in the same pipeline.

Below is an updated version of the script using code provided by Lucía Balestrazzi.

Additional information: The following with-statement is nested inside a larger for-loop that cycles through subsamples of an original genome. The first subsample fasta in my directory has a single sequence ">chr1:0-40129801" with a length of 40129801. I'm trying to write out a text file "OUTPUT" that has some basic information about each subsample fasta. This text file will be used as an input for another program downstream.

Header names in the original fasta file are chr1, chr2, etc... while the header names in the subsample fastas are something along the lines of:

batch1.fa >chr1:0-40k

batch2.fa >chr1:40k-80k

...etc...

import Bio.SeqIO as IO

record_dict = IO.to_dict(IO.parse(ORIGINAL_GENOME, "fasta")) #not the subsample
with open(GENOME_SUBSAMPLE, 'r') as FIN: 
    for LINE in FIN:
        if LINE.startswith('>'):
                #Example of "LINE"... >chr1:0-40129801
                HEADER = re.sub('>','',LINE)
                #HEADER = chr1:0-40129801
                HEADER2 = re.sub('\n','',HEADER)
                #HEADER2 = chr1:0-40129801 (no return character on the end)
                CONTIG = HEADER2.split(":")[0]
                #CONTIG = chr1 
                PART2_HEADER = HEADER2.split(":")[1]
                #PART2_HEADER = 0-40129801 
                START = int(PART2_HEADER.split("-")[0])
                #START = 0
                END = int(PART2_HEADER.split("-")[1])
                #END = 40129801
                LENGTH = END-START
                #LENGTH = 40129801 minus 0 = 40129801
            
                #This is where I'm stuck...
                ORIGINAL_CONTIG_LENGTH = len(record_dict[CONTIG]) #This returns "KeyError: 1"
                #ORIGINAL_CONTIG_LENGTH = 223705999 (this is from the full genome, not the subsample). 

                OUTPUT.write(str(START) + '\t' + str(HEADER2) + '\t' + str(LENGTH) + '\t' + str(CONTIG) + '\t' + str(ORIGINAL_CONTIG_LENGTH) + '\n')
                #OUTPUT = 0    chr1:0-40129801    40129801    chr1    223705999
OUTPUT.close()

I'm relatively new to bioinformatics. I know I'm messing up on how I'm using the dictionary, but I'm not quite sure how to fix it.

Any advice would be greatly appreciated. Thanks!

biopython fasta bioinformatics python
2 Answers

You can do it this way:

import Bio.SeqIO as IO
record_dict = IO.to_dict(IO.parse("genome.fa", "fasta"))
print(len(record_dict["chr1"]))

or

import Bio.SeqIO as IO
record_dict = IO.to_dict(IO.parse("genome.fa", "fasta"))
seq = record_dict["chr1"]
print(len(seq))

EDIT: Alternative code

import Bio.SeqIO as IO

record_dict = IO.to_dict(IO.parse("genome.fa", "fasta")
names = record_dict.keys()
for HEADER in names:
    #HEADER = chr1:0-40129801
    ORIGINAL_CONTIG_LENGTH = len(record_dict[HEADER])
    CONTIG = HEADER.split(":")[0]
    #CONTIG = chr1 
    PART2_HEADER = HEADER.split(":")[1]
    #PART2_HEADER = 0-40129801 
    START = int(PART2_HEADER.split("-")[0])
    END = int(PART2_HEADER.split("-")[1])
    LENGTH = END-START

The idea is that you define the dict once, get the value of its keys (all the contigs headers) and store them as a variable, and then loop through the headers extracting the info you need. No need to loop through the file.

Cheers

user_95007790

This works, just changed the "CONTIG" variable to a string. Thanks Lucía for all your help the last couple of days!

import Bio.SeqIO as IO

record_dict = IO.to_dict(IO.parse(ORIGINAL_GENOME, "fasta")) #not the subsample
with open(GENOME_SUBSAMPLE, 'r') as FIN: 
    for LINE in FIN:
        if LINE.startswith('>'):
                #Example of "LINE"... >chr1:0-40129801
                HEADER = re.sub('>','',LINE)
                #HEADER = chr1:0-40129801
                HEADER2 = re.sub('\n','',HEADER)
                #HEADER2 = chr1:0-40129801 (no return character on the end)
                CONTIG = HEADER2.split(":")[0]
                #CONTIG = chr1 
                PART2_HEADER = HEADER2.split(":")[1]
                #PART2_HEADER = 0-40129801 
                START = int(PART2_HEADER.split("-")[0])
                #START = 0
                END = int(PART2_HEADER.split("-")[1])
                #END = 40129801
                LENGTH = END-START
                #LENGTH = 40129801 minus 0 = 40129801
            
                #This is where I'm stuck...
                ORIGINAL_CONTIG_LENGTH = len(record_dict[str(CONTIG)]) 
                #ORIGINAL_CONTIG_LENGTH = 223705999 (this is from the full genome, not the subsample). 

                OUTPUT.write(str(START) + '\t' + str(HEADER2) + '\t' + str(LENGTH) + '\t' + str(CONTIG) + '\t' + str(ORIGINAL_CONTIG_LENGTH) + '\n')
                #OUTPUT = 0    chr1:0-40129801    40129801    chr1    223705999
OUTPUT.close()
user_76148360
Related