seqTools/TESTcalcN50N90.py at master · kdiverson/seqTools · GitHub

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
#!/usr/bin/python
import sys

infile = open(sys.argv[1], 'r')

def get_next_fasta (fileObject):
    '''usage: for header, seq in get_next_fasta(fileObject):
    '''
    header = ''
    seq = ''
    for line in fileObject:
        if line.startswith('>'):
            header = line.strip()
            break

    for line in fileObject:
        if line.startswith('>'):
            yield header, seq
            header = line.strip()
            seq = ''
        else:
            seq += line.strip()
    #yield the last entry
    if header:
        yield header, seq

contigs = []
total = 0

for header, seq in get_next_fasta(infile):
   length = len(seq)
   contigs.append(length)
   total += length

halfTot = 0
contigs.sort()

for contig in contigs:
    halfTot += contig
    if halfTot>total/2:
        print "N50:", contig
        break