1

I have a fasta file with millions of DNA sequences - for each DNA sequence, there is an ID as well.

>seq1
AATCAG #----> 5mers of this line are AATCA and ATCAG
>seq3
AAATTACTACTCTCTA
>seq19
ATTACG #----> 5mers of this line are ATTAC and TTACG

what I want to do is that if the length of DNA sequence is 6 then I make 5mers of that line (shown in the code and above). So the thing that I have could not solve is that I want to make a dictionary in a way that dictionary shows which 5mers come from which sequence id.

here is my code but it is halfway through:

from Bio import SeqIO
from Bio.Seq import Seq

with open(file, 'r') as f:
    lst = []
    dict = {}
    for record in SeqIO.parse(f, 'fasta'): #reads the fast file 
        if len(record.seq) == 6:  #considers the DNA sequence length
            for i in range(len(record.seq)):
                kmer = str(record.seq[i:i + 5])
                if len(kmer) == 5:
                    lst.append(kmer)
                    dict[record.id] = kmer  #record.id picks the ids of DNA sequences


    #print(lst)
                    print(dict)

the desired dictionary should look like:

dict = {'seq1':'AATCA', 'seq1':'ATCAG', 'seq19':'ATTAC','seq19':'TTACG'}
Apex
  • 1,055
  • 4
  • 22
  • 1
    You can't have duplicate keys in a dictionary e.g. `"seq19"`. – Matthias Mar 23 '21 at 09:48
  • is there an alternative way? – Apex Mar 23 '21 at 09:49
  • Use list instead of dictionary –  Mar 23 '21 at 09:53
  • just to mention that - if I have a list or a dictionary showing which 5mer comes from which sequence id - in the downstream analysis I need to loop over each sequence id and related 5mers – Apex Mar 23 '21 at 09:57
  • See if it helps.. https://stackoverflow.com/questions/10664856/make-a-dictionary-with-duplicate-keys-in-python –  Mar 23 '21 at 10:00

1 Answers1

2

using defaultdict as suggested in Make a dictionary with duplicate keys in Python by @SulemanElahi because You can't have duplicate keys in a dictionary

from Bio import SeqIO
from collections import defaultdict

file = 'test.faa'

with open(file, 'r') as f:
    dicti = defaultdict(list)
    for record in SeqIO.parse(f, 'fasta'): #reads the fast file 
        if len(record.seq) == 6:  #considers the DNA sequence length
                kmer = str(record.seq[:5])
                kmer2 = str(record.seq[1:])
                dicti[record.id].extend((kmer,kmer2))  #record.id picks the ids of DNA sequences

print(dicti)

for i in dicti:
    for ii in dicti[i]:
        print(i, '   : ', ii) 

output:

defaultdict(<class 'list'>, {'seq1': ['AATCA', 'ATCAG'], 'seq19': ['ATTAC', 'TTACG']})
seq1    :  AATCA
seq1    :  ATCAG
seq19    :  ATTAC
seq19    :  TTACG
pippo1980
  • 2,181
  • 3
  • 14
  • 30