I am writing my own parser for fasta format. I can't use Biopython or anything else, because it's a part of an assigment and our teacher wants us to try to do it manually. For now, I have done this:
def read_file(fasta_file): "Parse fasta file" count = 0 headers =  sequences =  aux =  with open("yeast.fna", 'r') as infile: for line in infile: record = line.rstrip() if record and record == '>': headers.append(record[1:]) if count > 0: sequences.append(''.join(aux)) aux =  else: aux.append(record) count += 1 sequences.append(''.join(aux)) return headers, sequences headers, sequences = read_file("yeast.fna") lengths = 0 countN = 0 acum = 0 lengths =  n = 0 lengths = [len(entry) for entry in sequences] count = sum(lengths) seqlengths = [entry for entry in sequences] countN = [entry.count("N") + entry.count("n") for entry in sequences] print("\nTotal length of the assembly : ", count) print("\n Total number of unsequenced nucleotides (Ns) in the assembly : ", countN) #arrange the secq sizes in decrossing order all_len = sorted(lengths, reverse=True) print(all_len) #Search for size accumulation until the size is <= half the size for y in range (len(lengths)): if acum <= count/2: acum = all_len[y] + acum n = y #L50 n = n + 1 # because the vector begin by a 0 print("The L50 is", n) print("The N50 is", all_len[n-1]) #Seq is n-1 because vector begin in 0
As you can see this parser does some basic stats like calculating L50, ... But the problem with this parser is that it is not very universal since I designed it for the assembled genome of S. Cerivisae which has 17 chromosomes. I know that the best way to do a universal effective parser is to use a dictionary (and make 2 dictionaries).
Can you help me with that?
Here is a sample of input file
yeast.fna. It is basically a txt file with headers which begin by a
> followed by the dna sequences. I just copy part of it because it's a very big txt file with 17 chromosomes sequences (each one has it own header).
>NC_001133.9 Saccharomyces cerevisiae S288C chromosome I, complete sequence ccacaccacacccacacacccacacaccacaccacacaccacaccacacccacacacacacatCCTAACACTACCCTAAC ACAGCCCTAATCTAACCCTGGCCAACCTGTCTCTCAACTTACCCTCCATTACCCTGCCTCCACTCGTTACCCTGTCCCAT TCAACCATACCACTCCGAACCACCATCCATCCCTCTACTTACTACCACTCACCCACCGTTACCCTCCAATTACCCATATC >NC_001134.8 Saccharomyces cerevisiae S288C chromosome II, complete sequence AAATAGCCCTCATGTACGTCTCCTCCAAGCCCTGTTGTCTCTTACCCGGATGTTCAACCAAAAGCTACTTACtaccttta ttttatgtttactttttatagGTTGTCTTTTTATCCCACTTCTTCGCACTTGTCTCTCGCTACTGCCGTGCAACAAACAC TAAATcaaaacaatgaaataCTACTACATCAAAACGCATTTTccctagaaaaaaaattttcttacAATATACTATACTAC ACAATACATAATCACTGACTTTCgtaacaacaatttccttcacTCTCCAACTTCTCTGCTCGAATCTCTACatagtaata >NC_001148.4 Saccharomyces cerevisiae S288C chromosome XVI, complete sequence AAATAGCCCTCATGTACGTCTCCTCCAAGCCCTGTTGTCTCTTACCCGGATGTTCAACCAAAAGCTACTTACtaccttta ttttatgtttactttttatagaTTGTCTTTTTATCCTACTCTTTCCCACTTGTCTCTCGCTACTGCCGTGCAACAAACAC TAAATCAAAACAGTGAAATACTACTACATCAAAACGCATATTccctagaaaaaaaaatttcttacaATATACTATACTAC
Here is my output when I run my code :
Total length of the assembly : 12157105 Total number of unsequenced nucleotides (Ns) in the assembly : [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] [1531933, 1091291, 1090940, 1078177, 948066, 924431, 813184, 784333, 745751, 666816, 576874, 562643, 439888, 316620, 270161, 230218, 85779] The L50 is 6 The N50 is 924431
Everything is okay but I need to make a more universal fasta parser which doesn't work only for my specific input file with 17 chromosome. In order to do that I have to use dictionary but I really don't know how to do that. You can also download the full input file. At the top of the page in line "Download sequences in FASTA format for genome, transcript, protein" you click on "genome" and you will download the file which is zipped.