-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathkmerfreq.py
70 lines (35 loc) · 1.09 KB
/
kmerfreq.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#!/usr/bin/env python3
#r17 finding if any common motifs before numpts using kmer freqs
import sys
import os
import re
import glob
import subprocess as sp
import concurrent.futures
from Bio import SeqIO
import gzip
import math
import collections
def tokenize(filename):
digits = re.compile(r'(\d+)')
return tuple(int(token) if match else token for token, match in ((fragment, digits.search(fragment)) for fragment in digits.split(filename)))
def kmerfreq(dna_sequence,wordsize):
wordsizeint=int(wordsize)
kmers=[]
seqlen=float(len(dna_sequence))
wordsize=float(wordsize)
#posscount=float((4**wordsize)*(1-math.exp(-(1/(4**wordsize))*seqlen)))
#print(posscount)
for x in range(len(dna_sequence)-wordsizeint):
kbit=dna_sequence[x:x+wordsizeint]
kmers.append(kbit) #also try if to test before set see which is faster:
kmerset=set(kmers)
for i in kmerset:
print(i,kmers.count(i))
infile=sys.argv[1]
wordn=int(sys.argv[2])
if __name__ == '__main__':
seq1=""
for x in SeqIO.parse(infile,'fasta'):
seq1=seq1+str(x.seq)
kmerfreq(seq1,wordn)