 |
파이썬 마을
우리나라 파이썬 사용자들의 이야기 마을
|
|
| 이전 주제 보기 :: 다음 주제 보기 |
| 글쓴이 |
메시지 |
mgenome
가입:
올린 글: 59
|
 올려짐: 2007 9월 21 6:16 pm 주제: 대용량 화일에서 필요한 부분만 추출 |
 |
|
추석 기념으로 오랜만에 프로그램 좀 해 보려는 프로그래밍 초보 생물학자 입니다.
20M 정도 되는 화일에서 gi number라는 것을 뽑아내어서 이 gi number를 280M 정도 되는 화일에서 검색해서 gi number를 가지고 있는 line만 뽑아내려고 하는데 잘 안되네요.
280M 화일의 형태는
gi_no 'space' taxid_no
gi_no 'space' taxid_no
gi_no 'space' taxid_no
.
.
.
.
와 같이 두 숫자가 공백으로 연결된 아주 단순한 형태입니다.
20M 정도되는 화일에서는 gi number를 다 뽑아내어서 (58000 개정도) 공백으로 연결해서 txt화일을 만들었구요 500KB 정도 크기가 되네요. .
처음에는 각각의 number를 280M 화일을 한 라인씩 읽은 것과 비교하려고 했으나 (re module을 사용해서) 속도 면에서 말도 안되는 짓이라는 생각이 들어서.
일단 280M 화일에서 라인을 하나하나 읽어서 500KB되는 gi number 화일과 비교해 보려고 했으나 이것도 속도 면에서 불가능 한 것 같아서 실행 화일 돌리다가 포기 했습니다.
다행이 280M 화일에서 gi number는 순서대로 써져 있기 때문에 이것을 이용하면 빠르게 검색이 되겠다는 생각이 들어서 프로그램을 다음과 같이 작성해서 실헝을 시켜 봤습니다. 숫자 sorting 하는 방법을 약간 변형해서 line을 검색하는 형식으로 했는데요. 잘 되는 듯 하다가 (한 계속해서 무한 루프를 돌고 있는 것 같습니다. (지금도 프로그램이 돌고 있네요 TT.)
| 코드: |
import re
import string
def extract_taxid(taxidfile, gifile, savefile):
## 첫번째 화일이 taxid화일, 두번째가 gi 화일
tfile = open(taxidfile)
gfile = open(gifile)
sfile = open(savefile, 'w')
gifile_str = gfile.read()
gfile.close()
tfile_len = len(tfile.read())
## 최대한 gi_number를 빨리 찾기 위해서 binary sorting하는 방법을 응용하려고 했는데. 너무 주먹구구식인 것 같다는 생각이 자꾸만 뜨네요.
tup = gifile_str.split()
for gi in tup:
index=frag=tfile_len/2
tfile.seek(index)
tfile.readline()
b=tfile.readline().split()
while int(b[0]) <> int(gi):
frag = frag/2
if int(b[0]) > int(gi):
index = index - frag
if int(b[0]) < int(gi):
index = index + frag
tfile.seek(index)
tfile.readline()
b=tfile.readline().split()
sfile.write(string.join(b, " ") + '\n')
tfile.close()
sfile.close()
|
결과로 나온 화일을 보면
.
.
.
.
30387422 215158
30387423 215158
30387424 215158
30387425 215158
30387426 215158
30387427 215158
30387428 215158
30387454 227720
30387455 227720
30387456 227720
30387457 227720
30387458 227720
30387459 227720
30387460 227720
30387461 227720
30387462 22
이런식으로 나왔네요. 앞쪽은 gi number고 뒤쪽은 taxid number 인데.
한 20000개 정도까지는 잘 나왔는데 그 뒤로 멈추어 버렸네요. 아마도 sorting 방법에 문제가 있어서 그 라인을 못 찾는 것 같다는 생각이 드는데요.
line을 직접적으로 지정해서 검색을 할 수 있으면 좋을 텐데 그렇게 못하고 그냥 화일 사이즈를 가지고 하다보니까 검색 범위를 좁이지 못했거나, 다른 곳을 계속 검색하는 것 같네요.
숫자가 쓰이다 만 것을 보면 다른 문제가 있는 것 같기도 하고요.
어떻게 하면 해결 할 수 있을까요?
아마도 멋진 다른 방법이 있을 듯 하지만... 제가 전공이 생물학이라서 더이상 나아갈 방법이 없어서 질문을 합니다.
고수의 1분 커맨트가 초보의 1주일을 좌우한다고 생각하시고 많은 지도를 부탁드립니다.  |
|
| 위로 |
|
 |
장혜식
사이트 관리자
가입:
올린 글: 1374
|
| 올려짐: 2007 9월 21 11:43 pm 주제: |
|
|
말씀하신 큰 파일이 항상 변경되는 경우가 아니라면, 일단 인덱스가 된 형태로
만들어서 변환한 다음 읽는 방법이 실험할 때 가장 편합니다.
보통 자료처리할 때 아주 편하게 쓰이는 모듈이 shelve가 있는데요.
import shelve
db = shelve.open('파일이름')
db[gi_number] = 데이터
이런 식으로 한 번만 쭉 읽어놓으면 다음부터 같은 파일을 shelve.open으로
열어서 쓰면 그 내용에 계속 접근할 수 있습니다.
실제 사용 예로 EMBL CDS fasta파일을 읽어서 E. coli에서 나온 것만 골라서
뽑아주는 소스는 이렇게 만들 수 있습니다.
| 코드: |
from Bio import SeqIO
import shelve
seqs = SeqIO.parse(open('/data/embl/cds_nr.fasta'), format='fasta')
db = shelve.open('ecoli.db')
for rec in seqs:
if ('Escherichia coli' in rec.description and
'hypothetical' not in rec.description):
name = rec.description.split()[0].split(':')[1]
seq = rec.seq.tostring()
db[name] = seq
|
나중에 EMBL ID로 검색할 때는 다음과 같이 open만 하면 검색할 수 있습니다.
| 코드: |
import shelve
db = shelve.open('ecoli.db')
genes = ['100746833', '1008171518']
for gene in genes:
print db[gene]
|
GenBank에서 받은 32기가짜리 파일을 대상으로 shelve에 넣어보니까
대략 20분이면 되는 걸로 봐서는, 400메가면 금방 될 것 같습니다. 
_________________
http://me2day.net/hyeshik |
|
| 위로 |
|
|
forest
가입:
올린 글: 66
소속: 연구소
|
| 올려짐: 2007 10월 02 10:50 am 주제: 저에게도 유용한 정보네요. |
|
|
저에게도 유용한 정보네요. 감사합니다. 사전과 비슷하면서도 리소스는 훨씬 적게 잡아먹을것 같습니다. 대용량 자료를 사전으로 피클링 할 때는 여러 파일로 쪼개 놓곤 했는데, 좋은 모듈이 있었네요. 사전과 비교했을 때 단점이나 다른 장점도 있을듯 한데 좀 더 공부해봐야겠습니다.
_________________
:: R U Ready? :: |
|
| 위로 |
|
|
|
|
새로운 주제를 올릴 수 없습니다
답글을 올릴 수 없습니다
주제를 수정할 수 없습니다
올린 글을 삭제할 수 없습니다
투표를 할 수 없습니다
|
|
FAQ
검색
멤버리스트
사용자 그룹
사용자 등록하기
개인 정보
비공개 메시지를 확인하려면 로그인하십시오
로그인
