PageRenderTime 20ms CodeModel.GetById 11ms app.highlight 6ms RepoModel.GetById 2ms app.codeStats 0ms

/scripts/metagenomics/convert_title.py

https://bitbucket.org/cistrome/cistrome-harvard/
Python | 39 lines | 35 code | 1 blank | 3 comment | 1 complexity | 5140bdd51b315c55e0e435c1d7adff06 MD5 | raw file
 1#!/usr/bin/env python
 2
 3"""
 4convert nt and wgs data (fasta format) to giNumber_seqLen
 5run formatdb in the command line: gunzip -c nt.gz |formatdb -i stdin -p F -n "nt.chunk" -v 2000
 6"""
 7
 8import os, sys, math
 9
10if __name__ == '__main__':
11    
12    seq = [] 
13    len_seq = 0
14    invalid_lines = 0
15    
16    for i, line in enumerate(sys.stdin):
17        line = line.rstrip('\r\n')
18        if line.startswith('>'):
19            if len_seq > 0:
20                print ">%s_%d" %(gi, len_seq)
21                print "\n".join(seq)
22            title = line
23            fields = title.split('|')
24            if len(fields) >= 2 and fields[0] == '>gi':
25                gi = fields[1]
26            else:
27                gi = 'giunknown'
28                invalid_lines += 1
29            len_seq = 0
30            seq = []
31        else:
32            seq.append(line)
33            len_seq += len(line)
34    if len_seq > 0:
35        print ">%s_%d" %(gi, len_seq)
36        print "\n".join(seq)
37        
38    print >> sys.stderr, "Unable to find gi number for %d sequences, the title is replaced as giunknown" %(invalid_lines)
39