그럼 이제 이 파일에서 무엇인가 뽑아내야 겠죠?
무엇을 뽑아내느냐?
16S rRNA를 뽑아낼겁니다.
어떻게?
다음 스크립트를 작성해 봅시다.
$vi parser_nt.py
import glob,sys,re,gzip
from Bio import SeqIO
try:
input_fa = sys.argv[1]
except:
print "fasta_split.py <in.fa.gz>"
exit(1)
for rec in SeqIO.parse(gzip.open(input_fa),format='fasta'):
desc = rec.description
seq = str(rec.seq)
name = desc.strip().split('\x01')[0]
if name.upper().find(' 16S RIBOSOMAL RNA ') != -1 or name.find(' 16S RIBOSOMAL RNA,') != -1:
if name.upper().find('MITOCHONDRIA')!= -1:
pass
elif name.upper().find('CHLOROPLAST')!= -1:
pass
elif name.upper().find(' 23S ') != -1:
pass
else:
if name.upper().find(' SIMILAR ') != -1:
pass
elif name.upper().find(' INTRON') != -1:
pass
elif name.upper().find(' PLASTID') != -1:
pass
else:
print '>%s\n%s\n'%(desc,seq)
※ 위의 스크립트를 수정해서 입맛에 맞게 교정하시고 사용하시면되겠습니다.
$python parser_nt.py nt.gz > nt.fasta
위에 스크립트 실행시키면 떡 하니 수 gb 짜리 파일이 하나 나올겁니다.
이 파일안에는 nt서열 중에서 16S rRNA 서열 (대신 mitochondria와 chloroplast의 16S rRNA는 제외하고 이것도 16S rRNA지만 저한테는 일단 필요없어서 뺐습니다. 사용하고 싶으시면 사용하셔도 됩니다. :) )이 500만개 정도 들어 있습니다.
이 파일에는 온갖 종의 16S rRNA 서열이 있지만 문제가 있다는 점!
16S rRNA 서열이 품질이...
어떤 녀석은 full 서열이 있지만 어떤녀석은 서열의 일부만 가지고 있는 경우가 있습니다.
그래서 그런 녀석들을 잘 확인해서 제거를 하던지 merge를 하던지...
그건 개인 취향으로 남겨 놓도록 하겠습니다. :)
aka 필터링을 하던 안하던 문제 생기면 그 문제의 책임은 오롯이 당신의 것!!
그리고 즐거운 추석보내시기 바랍니다. :)
 |
| 출처 SM |
