Python：fastaファイルからkmersをカウントする

Question

fastaファイルからkmersを数えたいと思う。

import operator 

seq = open('file', 'r') 
kmers = {} 
k = 5 
for i in range(len(seq) - k + 1): 
    kmer = seq[i:i+k] 
    if kmer in kmers: 
     kmers[kmer] += 1 
    else: 
     kmers[kmer] = 1 

for kmer, count in kmers.items(): 
    print (kmer + "\t" + str(count)) 

sortedKmer = sorted(kmers.items(), key=itemgetter(1), reverse=True) 
for item in sortedKmer: 
    print (item[0] + "\t" + str(item[1])) 

これは1つの配列のみを持つファイルを正常に動作しますが、今、私はいくつかのコンティグとFASTAファイルを持っている：私は、次のスクリプトを持っています。 マイFASTAファイルは次のようになります。

>1 
GTCTTCCGGCGAGCGGGCTTTTCACCCGCTTTATCGTTACTTATGTCAGCATTCGCACTT 
CTGATACCTCCAGCAACCCTCACAGGCCACCTTCGCAGGCTTACAGAACGCTCCCCTACC 
CAACAACGCATAAACGTCGCTGCCGCAGCTTCGGTGCATGGTTTAGCCCCGTTACATCTT 
CCGCGCAGGCCGACTCGACCAGTGAGCTATTACGCTTTCTTTAAATGATGGCTGCTTCTA 
AGCCAACATCCTGGCTGTCTGG 
>2 
AAAGAAAGCGTAATAGCTCACTGGTCGAGTCGGCCTGCGCGGAAGATGTAACGGGGCTAA 
ACCATGCACCGAAGCTGCGGCAGCGACACTCAGGTGTTGTTGGGTAGGGGAGCGTTCTGT 
AAGCCTGTGAAGGTGGCCTGTGAGGGTTGCTGGAGGTATCAGAAGTGCGAATGCTGACAT 
AAGTAACGATAAAGCGGGTGAAAAGCCCGCTCGCCGGAAGACCAAGGGTTCCTGTCCAAC 
GTTAATCGGGGCAGG 

は、どのように私はそれが「> 1」の後の最初のシーケンスを取り、出力は、「> 2」に行くことを印刷する、などその出力を印刷することを、スクリプトを変更できますか？

Answer 1

私はキラーやファストについて聞いたことがありませんが、あなたがしようとしていることを理解していると思います。

'>'を含む正規表現を分割しようとすることはできますが、ファイルを1行ずつ処理し、kmersを蓄積してから '> 1'行まで適切に印刷することをお勧めします。コメント

import operator 

def printSeq(name, seq): 
    # Extract your code into a function and print header for current kmer 
    print("%s\n################################" %name) 
    kmers = {} 
    k = 5 
    for i in range(len(seq) - k + 1): 
     kmer = seq[i:i+k] 
     if kmer in kmers: 
      kmers[kmer] += 1 
     else: 
      kmers[kmer] = 1 

    for kmer, count in kmers.items(): 
     print (kmer + "\t" + str(count)) 

    sortedKmer = sorted(kmers.items(), reverse=True) 

    for item in sortedKmer: 
     print (item[0] + "\t" + str(item[1])) 

with open('file', 'r') as f: 
    seq = "" 
    key = "" 
    for line in f.readlines(): 
     # Loop over lines in file 
     if line.startswith(">"): 
      # if we get '>' it is time for a new sequence 
      if key and seq: 
       # if it wasn't the first we should print it before overwriting the variables 
       printSeq(key, seq) 
      # store name after '>' and reset sequence 
      key = line[1:].strip() 
      seq = "" 
     else: 
      # accumulate kmer until we hit another '>' 
      seq += line.strip() 
    # when we are done with all the lines, print the last sequence 
    printSeq(key, seq) 

Answer 2

とのコード私はあなたの例のFASTAファイルに次のことを試してみました、それが動作するはずです以下を参照してください。

def count_kmers(seq, k, kmers): 
    for i in range(len(seq) - k + 1): 
     kmr = seq[i:i + k] 

     if kmr in kmers: 
      kmers[kmr] += 1 

     else: 
      kmers[kmr] = 1 

filename = raw_input('File name/path: ') 
k = input('Value for k: ') 
kmers = {} 

# Put each line of the file into a list (avoid empty lines) 
with open(filename) as f: 
    lines = [l.strip() for l in f.readlines() if l.strip() != ''] 

# Find the line indices where a new sequence starts 
idx = [i for (i, l) in enumerate(lines) if l[0] == '>'] 

idx += [len(lines)] 

for i in xrange(len(idx) - 1): 
    start = idx[i] + 1 
    stop = idx[i + 1] 
    sequence = ''.join(lines[start:stop]) 

    count_kmers(sequence, k, kmers) 

print kmers 

はそれが役に立てば幸い:)