2012-07-16 70 views
1

我有兩個文件:第一個是帶有頭文件和序列的fasta文件,第二個文件僅包含頭文件。在python中匹配fasta文件中的頭文件

File_1:

>DF94KKQ1|265|D0M1LACXX|3|2103|4637|10742|1|N|0|TGACCA 
TTCCAAAGAAACATGGAAGACCCAGGACTTGGAGGCACCAGGCACCAGCACACAGGGGTA 
GGCACATGGCATGGTGTTGGTTGAAGTCTACTTTTCCCACC 
>DF94KKQ1|265|D0M1LACXX|3|2103|4565|10742|1|N|0|TGACCA 
TTAATTTTTTCAGGCAAGTTTTGTGGATTTCAGTGTGTAAGTCTTTCACCTCTTTGGTTA 
AATTTATTCCTATGTATTTTATTCCTTTAGATGCTATTATG 
>DF94KKQ1|265|D0M1LACXX|3|2103|4565|10742|2|N|0|TGACCA 
TTAATTTTTTCAGGCAAGTTTTGTGGATTTCAGTGTGTAAGTCTTTCACCTCTTTGGTTA 
AATTTATTCCTATGTATTTTATTCCTTTAGATGCTATTATG 

File_2:

>DF94KKQ1|265|D0M1LACXX|3|2103|4565|10742|2|N|0|TGACCA 
>DF94KKQ1|265|D0M1LACXX|3|2103|4668|10746|1|N|0|TGACCA 
>DF94KKQ1|265|D0M1LACXX|3|2103|4668|10746|2|N|0|TGACCA 
>DF94KKQ1|265|D0M1LACXX|1|2207|10852|3331|2|N|0|TGACCA 

我想在File_2頭與任何在File_1說,直到7具有完全相同的字符匹配 '|'。

我在File_1中分割項目(標題的每個部分都被編入一個列表)。苯胺是用「>」開始被放置到一個變量:

#!/usr/bin/env python 
import sys 
from Bio import SeqIO 

#Function, split header line into a list 
def getHeaderInfo(blastLine): 
    myFields = blastLine.strip("\n").split("|") 
    HeaderInfo = myFields[:6] 
    return HeaderInfo 

input_file = sys.argv[1] 

#Get input file from the command line 
inFileName = sys.argv[1] 

#open the input file 
inFileHandle = open(inFileName) 

#loop over the input file line by line 
for thisLine in inFileHandle.readlines(): 
    if thisLine [0] == '>': 
     print getHeaderInfo(thisLine) 
     HeaderInfo = getHeaderInfo(thisLine) 

我一直在試圖找到我可以在File_2比較這些相同的指標返回以下輸出的方法:

>DF94KKQ1|265|D0M1LACXX|3|2103|4565|10742|1|N|0|TGACCA 
TTAATTTTTTCAGGCAAGTTTTGTGGATTTCAGTGTGTAAGTCTTTCACCTCTTTGGTTA 
AATTTATTCCTATGTATTTTATTCCTTTAGATGCTATTATG 
>DF94KKQ1|265|D0M1LACXX|3|2103|4565|10742|2|N|0|TGACCA 
TTAATTTTTTCAGGCAAGTTTTGTGGATTTCAGTGTGTAAGTCTTTCACCTCTTTGGTTA 
AATTTATTCCTATGTATTTTATTCCTTTAGATGCTATTATG 

我嘗試過的幾種方法使用索引,但是,我的密鑰不是唯一的。我怎樣才能拿到前六個元素,並讓他們成爲我的鑰匙,還是有比我現在嘗試的更好的方法?謝謝。

+0

標頭線在File_1和File_2開頭「>」 ,他們沒有出現在文章 – user1529819 2012-07-16 19:28:52

回答

1

這是做你想做的嗎?

def make_key(line): 
    return "|".join(line.split("|", 7)[ : 7]) + "|" 

header_set = set() 
with open("file_2.txt") as in_f: 
    for line in in_f: 
     header_set.add(make_key(line)) 

with open("file_1.txt") as in_f, open("file_3.txt", "w") as out_f: 
    accept = False 
    for line in in_f: 
     if line.startswith(">"): 
      key = make_key(line) 
      accept = key in header_set 

     if accept: 
      out_f.write(line) 
-1

此方法涉及存儲器映射file_1.txt,通過它re.finditer運行線加載到由報頭鍵控defaultdict

import re 
import mmap 
from collections import defaultdict 
pat = re.compile('>.+?(?=(>|$))', re.DOTALL) 
with open('file_1.txt') as f: 
    map = mmap.mmap(f.fileno(), 0, access=mmap.ACCESS_READ) 

line_1 = defaultdict(list) 
for line in pat.finditer(map): 
    fields = line.group().split('|') 
    key = '|'.join(fields[:6]) 
    line_1[key].append(line.group()) 

map.close() 
line_2 = [] 
with open('file_2.txt') as f: 
    for line in f: 
     fields = line.split('|') 
     key = '|'.join(fields[:6]) 
     line_2.append(key) 

line_2 = set(line_2) 
for key in line_2.intersection(line_1.keys()): 
    print "".join(line_1[key]) 
+0

謝謝我最終確定瞭如何解析這兩個文件。回覆整個fasta唱片的頭和序列是我的下一個掛斷。將很快發佈代碼。 – user1529819 2012-07-17 21:42:25

+0

@ user1529819 - 上面的正則表達式解決方案返回標題和序列 - 除非我錯過了某些東西 – iruvar 2012-07-17 23:06:07

+0

實際上,grep結束了更快更簡單的工作。問題是文件1是巨大的,並且將它讀入內存是不可行的。使用grep我可以在編輯查詢文件(文件2)後進行模式匹配 – user1529819 2012-07-18 16:54:45