Microbiome Database를 만들어볼까? -NCBI편 4-

지난번 글에서 받기시작했던 nt.gz파일은 잘 받아졌나요?

그럼 이제 이 파일에서 무엇인가 뽑아내야 겠죠?

무엇을 뽑아내느냐?
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

Microbiome Database를 만들어볼까? -NCBI편 3-

지난 시간에 이곳에서 귀한 자료를 받아 봤을 겁니다.
그게 무엇이냐!!!

지금까지 공개된 bacteria의 서열들이죠 정확하게 말하면
현재까지 수많은 연구자들이 자발적(자의타의)으로
공개해준 것을 NCBI가 아름답게 정리한 RNA서열들만 다운받았습니다.

근데 우리는 RNA 서열이 아니라 rRNA 서열이 필요하죠

이전 글에서 우리는 rna라는 폴더 안에 12만개에 달하는 파일들을 다운받았습니다.
아마 그냥 ls하시면 어쩌구 저쩌구 long 할겁니다.
ls로 불러오기에 item이 너무 많다 이거죠
그럼 어쩌지?? @.@

우리에겐 귀도 훃님의 파이썬이 있지 않겠습니꽈

간단하게 다음과 같은 스크립트를 뚝딱 뚝딱 만들어보죠

import os,glob
from Bio import SeqIO
for files in glob.glob('rna/*gz`):
    for rec in SeqIO.parse(gzip.open(files), format='fasta'):
        name =rec.description
        seq = rec.seq
        if name.find('[product=16S ribosomal RNA]') != -1:
            print '>%s\n%s\n'.format(name,seq)

(python2.7에 Biopython이 설치되어 있어야 하고 *_rna_form_genomic.fna.gz 파일이 rna폴더 밑에 위치하고 있어야 합니다.)

$python script.py > ncbi.fa

라고 해주면 헤더에 "[product=16S ribosomal RNA]"가 포함된 aka 우리가 원하는 바로 그것! 16S rRNA 서열을 각 종에서 샤샤샥 ncbi.fa라는 파일에 저장할 수가 있습니다.

-12만개 파일 읽어 오는거라 순차적으로 하면 3-4시간 걸릴것이고
subprocess로 잘 해주시면 적어도 십수분? 30분이내면 충분히 끝날각 되겠습니다.

출처: SM Town

Microbiome Database를 만들어 볼까? -RDP편-

RDP
이름에서 똭 느껴지지 않으십니까?
Ribosomal Database Project

일반 fasta/gb 자료가 있고
trainset 데이터가 있습니다.

Bacteria unalign seq (fa,gb)
Archaea unalign seq (fa,gb)

Train Set 은 여기서 받으시면됩니다.

근데 RDP 같은 경우 사실 그냥 별 이유 없이 받아보는 겁니다. ㅋ
-저도 이게 어떻게 쓰일지 잘 모르겠어요 :)

출처 JYP

Microbiome Database를 만들어 볼까? -Silva편-

SILVA는 de.NBI 가 뒷배인 DB로 지속적으로
업데이트되고 있는 DB중 하나입니다.

서열 개수도 아마 가장 많을 것으로 생각됩니다.
대신에 그만큼 정리가 가장 안되어 있기도 합니다. ㅋ
-사실 이정도 양의 서열에 저정도 수준을 갖춰놓은것도 용하긴합니다.

그리고 작년? 제작년쯤에 SILVAngs라는 서비스
MG-RAST의 대항마쯤으로 런칭을 했는데 사견으로는
그때 이미 MG-RAST는 심심하면 서비스 다운되서 MG-RAST는
이미 논외 대상이었을것 같다능...(다년간 쌓아놓은 데이터가 문제였지..)

여튼 silva에서 제가 사용하고자 하는 파일은
silva에서 제공하는 qiime용 파일입니다.
>다운로드링크<

직접 찾아가려는 경우 [Download] > [Archive]를 클릭하시면
페이지에 온갖 리스트들이 나오는데 그중에서 [qiime]를 클릭하시면
제가 사용하고자 하는 파일이 보이실 겁니다. ;)

ㅋㅋ 당분간은 계속 다운로드만 주구장창 받는 글들입니다. :)

출처: JYP

Microbiome Database를 만들어볼까? -NCBI편 2-

지난 시간에 이곳에서 자료를 받아 봤습니다.

오늘은 또 다른 ncbi 자료를 받아 볼겁니다.

※지겹다고요? 제가 그랬잖아요? 당분간은 맨날 다운 받을 거라고 ㅋ
빨라야 추석 이후에야 무엇인가 하지 않을까 기대합니다.

지난번에는 ncbi에서 제공하는 모 그런걸 받았습니다.
-자세한 내용은 다음 글에... ㅎㅎ :)

그렇다면 이번에는 우리가 ncbi라고 하면
맨날 다운받았던 blast의 db를 만드는 source fasta파일인
nt.gz 을 받아보도록 하겠습니다.

자 웹브라우저에서 ftp://ftp.ncbi.nlm.nih.gov/blast/db/fasta/
들어가시면 nt.gz이 보입니다.

물론

$wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/fasta/nr.gz 

하셔서 다운받아도 무방합니다.

그럼 대략 40G 정도의 gzip파일을 받아봅시다. :)

출처: SM Town