간만에 denovo RNA-Seq 해보기 -유전자구조예측편-

denovo RNA-Seq를 사용해서 조립을 했다면
조립된 서열들은 어떤 유전자인지 궁금해 지겠쥬?

안 궁금하면 그냥 끝! 하고
서열을 NCBI에 fastq파일 디파짓하고 Bye さようなら하면
그냥 남 좋은일 하는 겁니다.
(나는 대인배다 나는 시퀀싱 비용이 아깝지 않다 하시는 분이라면
친하게 지내요!! 제발~ )

여튼 RNA-Seq을 했고, 생산된 RNA-Seq을 가지고 조립까지 했다면
조립된 서열들은 어떤 유전자들일까 궁금한게 인지사정!

그러면 그런 tool에서는 어떤 것들이 있을까?

바로 이런게 있습니다.
TransDecoder

TransDecoder Wiki

조립된 RNA-Seq서열 에서 coding 서열을 찾아주는 프로그램 입니다.
(현재 버전이 5.5.0이네요.. 다행히 어떤 업데이트도 일어나지 않았네요..)

풋 아마추어같이 RNA-Seq서열이니깐 ATG로 시작하는 것 찾으면 되지 무슨 프로그램이야 프로그램은 아마추어 같으니라고!!
라고 하신다면 당신은 느응력자!

다들 알고계시다 싶이 ATG로 시작하는 것들 major긴 하지만 RNA-Seq을 해서 조립하게되면 ATG로 시작하지 않은 partial로 어딘가가 짤려진 gene 서열들이 존재하기 때문에 그런것들도 잘 알아서(모 대략, 못찾는것도다는) 찾아주는 녀석이 바로 이녀석 되겠습니다.

-사실 이거 말고 다른것도 많이 있을겁니다. 제가 이것밖에 안써서 이거 소개합니다. ㅎㅎ


그냥 위에 파일 다운 받아서 압축 풀고 trinity로 조립한 fasta파일을 넣고 돌리면

$ ~/TransDecoder-TransDecoder-v5.5.0/TransDecoder.LongOrfs -t Trinity.fasta --gene_trans_map Trinity.fasta.gene_trans_map

(근데 저 --gene_trans_map이 무슨 옵션이었는지 까먹었네요...)

여튼 이렇게 돌리면 대략적인 결과 나오고 그 결과가지고 연구하면됩니다.
이거가지고 부족해!! 하시면 genome project 진행하시면되겠습니다!!

ps. 위의 글은 유전자 예측이 아닌 유전자 구조 예측이 맞는 표현입니다. Orz

출처: SM

qiime2 파헤치기

간만에 qiime2도 간만에 공부하려고 docs를 한번 읽어보는 겸해서 정리도 같이? :)

그럼 qiime2가 무엇이냐 우리 롭 훃님께서 만든 microbiome분석 파이프라인,
아니 파이프라인을 뛰어넘어 토탈 분석 패키지(?)이라고 할 수 있는 툴입죠.
*토탈패키지라고까지 얘기하는 이유는 qiime1을 사용해보신 분이라면 감동의 눈물이.. ㅋ

Core Concepts는...

qiime2는 qza, qzv와 같은 독자적인 파일 형식을 지원하고 있네요
qiime1에서는 듣도보지도 못한.. 개념이라서 처음에 artifacts라는 단어를 보고
artifacts는 제거해야 제맛아니겠슴까 했는데..  (   ' ')


일단 간단히 qiime2는 qiime1과 달리 지저분하게 파일을 만들지 않고 qza, qzv 두개의 파일로 관리를 시도 하고 있습니다.

qza는 텍스트 파일을 조금더 관리가 용이하도록 json이나 table 형식으로 만든것 같고
qzv는 시각화 파일로 plot이나 그래프를 나름에 규격에 맞게 qza파일과 같이 구조화된 형식으로 만들어놓고 시각화하기 편리하도록 구현해 놓았다고 보면 될것 같습니다.
https://view.qiime2.org라는 사이트를 운영해서 qzv파일만 있으면 언제든지 확인 할 수 있도록 만든것 보면 예상 가능하죠. 어느 멍청한 놈이 만들어도 우리 사이트에 qzv파일 올리면 내용 볼수 있게 우리가 만들었어 라고 하고 있죠.

그래서 qiime2 처음 나왔을때 상용화하려고 하려나 하는 생각을 했음죠
예전에는 그냥 순수 bio업자들에게 친숙하지 않은 plain txt파일이나 csv파일만
디립다 생성하고 있는데 파일 내용을 굳이 보여주지 않지만
웹페이지에 파일 가져다 놓으면 이쁜 그림이 뙇!!

두번째로 공식적으로 plugin을 사용하여 확장성을... (이것 또한 상용화의 냄새가...)
잡았다고 얘기하고 싶겠지만... 일단 제가 qiime2를 한번도 안써봐서
모 내머리에서 번뜩 생각나는 것들인 이미 plugin으로 만들어놨을듯한..
어차피 루틴하게 쓸만한 plugin은 이미 다 있을테니 가져다가 사용하면 될듯 합니다.
plugin개발하고 싶으시면 >> 여기로 (전 너님이 개발하면 유용하게 쓰겠습니다. ㅎㅎ )

모 처음으로 qiime2 첫 페이지를 읽어 보았습니다.

다음에는 install qiime2를 해보고 싶은데..
사용할수 있는 서버가 닫혀있어서...
테스트를 해볼 수가 없네요 ㅎㅎ :)

여튼 사적으로 접근 가능한 서버가 열리면 한번 테스트해서 자세한 리뷰한번 해보는걸로 :)

ps. 이제 고대 유물이 된 qiime1 base로 만든 파이프라인은 어디에 팔아먹나...

출처: SM

간만에 denovo RNA-Seq 해보기 -설치편-

최근 간만에 해보기가 올라가고 있는데...
진짜 2년만에 RNA-seq 분석을 해봐서..

걍 분석하는 단계나 프로그램 사용법 정리 차원에서 글을 올리고 있습니다.

4짜 산업 시대에 발맞춰 유전체 데이터 전문 설거지팀 하나 꾸리는것도 나쁘지 않을듯.... (대신 건당 비용때문에 수주가 안들어올 것 같다는게 함정 ㅎㅎ )

여튼 오늘은 de novo RNA-Seq 분석입니다.

일단 de novo RNAseq 시장을 석권했던.. 지금도 지배하고 있는 것으로 보이는데..
제가 사용했던 버전은 2.0.6이었는데.. ㄷㄷㄷ 벌써 2.8.4네요..
다들 아시는 삼위일체 Trinity 입니다.

지금 사용하는 서버에서는 cmake버전이 2.x라서 2.8.4대신 낮은 버전인 2.6.6버전으로 테스트를 수행하고 있습니다.
같은 input에 옵션이 비슷한데 2.6과 2.8의 결과가 많이 달라질지는 잘 모르겠습니다.
버전별 output 비교는 나중에 한번 기회되면 도전해보는것으로!!

$ wget https://github.com/trinityrnaseq/trinityrnaseq/archive/Trinity-v2.6.6.tar.gz
$ tar zxf Trinity-v2.6.6.tar.gz
$ cd trinityrnaseq-Trinity-v2.6.6/
$ make && make install

참고로 make했을때 어쩌구 저쩌구 /usr/local/bin 권한없다라는 메세지를 보여주고 에러를 밷어낸다면 trinityrnaseq-Trinity-v2.6.6/util/support_scripts/ 밑에 있는 trinity_installer.py 파일의 destination_package_dir 변수명의 내용을 수정해주시면됩니다.
(제 경우 make할때 DESTDIR 설정을 해주어도 계속 /usr/local/bin을 요구해서... trinity_install.py 파일을 직접 수정했습니다. ㅎㅎ 다른 방법이 분명 있을거 같은데.. )

여튼 에러가 발생한다면 해당 에러를 잡고 설치하면(당연한 소리를..) 문제 없을것이라고 말씀드릴 수 있습니다!!

출처: SM 

간만에 RNAseq 분석 해보기 -Reference편-

Alignment를 수행하기 위해서는 reference가 필요합니다.
모 어떤 alignment 툴에서는 그냥 genome 서열만 있어도 되지만
하이 쓰루풋 시퀀싱 데이터를 다룰 때는 대부분 genome 서열을
나름의 index를 새로 생성하게 됩니다.

앞에서 설치했던 aligner들의 index를 만드는 작업의 로그를 남겨보도록 하겠습니다.

BWA

$ ~/bwa/bwa index -p index_name genome.fa

hisat2

$ ~/hisat2/hisat2_extract_exons.py genome.gtf > genome.exon
$ ~/hisat2/hisat2_extract_splice_sites.py genome.gtf > genome.ss
$ ~/hisat2/hisat2-build -f genome.fa --ss genome.ss --exon genome.exon genome_index_base

STAR

$ ~/STAR/STAR --runThreadN 16 --runMode genomeGenerate --genomeDir genomeOutFolder --genomeFastaFiles genome.fa --genomeSAindexNbases index_base --sjdbGTFfile genome.gtf --sjdbOverhang 99

Kallisto

$ ~/tophat/gtf_to_fasta genome.gtf genome.fa genome.gtf2fa.fa
$ ~/kallisto/kallisto index --index=index_name genome.gtf2fa.fa

Salmon

$ ~/salmon/bin/salmon index -t genome.gtf2fa.fa -i genome_idx --type quasi -k 31

이렇게 하면 각 align tool을 사용하기 위한 reference는 준비되었습니다.

출처: SM

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

간만에 Microbiome DB 관련 글을 올립니다!
그래봤자 또 다운로드하자 되겠습니다. :)

오늘은 국내 microbiom 기업중 가장 기술력을 가지고 있는
업체인 chunlab의 DB 되겠습니다.

이 DB의 경우 당연히 fee를 내셔야합니다.

라이센스를 정확히 안읽어봤는데..
비영리면 무료였나? 정확하진 않지만
무엇인가 영리 서비스 하고자 EzBioCloud 16S database를 사용하려면
fee를 내야합니다.

다운로드는 다음 링크에 가셔서 다운로드 받으시면 됩니다.

>이곳<

천랩의 DB의 경우 제가 사용하겠다 말겠다할수 있는건 아니라서..
다운로드 받고싶으신 분들만 받아서 잘 사용하시면 되겠습니다.
qiime으로 분석하시는 분들께서는 별 무리없이 사용하시는 파이프라인에 붙여서 사용하실 수 있게 제작되었습니다. 사실 그냥 작동합니다. :)

출처: SM

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를 만들어볼까? -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

Microbiome Database를 만들어 볼까?

Microbiome분석을 위해서 여기저기 기웃거려봤다면
여러가지 16S rRNA 데이터베이스가 있다는것을 아실겁니다.

보통 microbiome분석에 입문해서 사용하는 것이라면
대게 처음 분석하는 tool에따라 결정되는데
우리 롭횽님의 qiime를 접한다면 greengene을, mothur을 접하게된다면 silva를 database로 만나게됩니다.

사실 대부분의 연구 결과들이 greengene과 silva로 나오기 떄문에 이 두 database를 사용하면 당연히 그 누구도 갠세이 놓지 않습니다.
-아니면 에디터나 리뷰어에게 외쳐보자. Drop the DB, yo!

근데 매번 분석하다보면 family수준 밑에만 내려가면 unknown은 왜이리 많을걸까..

그렇다면 그냥 우리가 손수 microbiome분석을 위한 db를 만들어보면 어떨까?

당근 이렇게 만들경우 실제 연구에 사용하기는 마뜩치 않다는걸 미리 말씀드립니다.

일단 tree를 만들기가;;; 녹녹치 않습니다. (물론 서열수를 줄이면 수십G 메모리를 가지는 워크스테이션이 있으면 가능합니다.)
그리고 제가 분류학자도 아니고 제가 서열보고도 얘가 몬지 알지도 못하고
서열에 taxonomy 붙여도 제대로 연결시킨건지 확인이 되지 않는다는 큰 문제가 있습죠

#물론_이사진을보면_누군지_압니다, 출처:SM Town

그럼에도 불구하고 왜하냐? 그냥 재미삼아, 경험삼아 만들어보는것입니다.

내가 사용하는 DB를 만드는데 얼마나 많은 고민이 녹아있고
얼마나 많은 생각들이 들어가 있는지 이해를 해보는것도 나쁘지 않을듯하고요 ㅎㅎ ;)

일단 이런걸 하겠다고 블로그에 띄워놨으니 언젠가는 후속글을 올리지 않을까요?
따라가능하도록 소스같은것은 각 글이나 github에 업로드하는 걸로 :)
jupyter notebook으로 올리면 더더욱 좋겠지만 제가 아직 notebook이 익숙치가 않아서..

우선 다음 글에서는 Custom Microbiome Database에 필요한 기초 자료 수집에 관련된 내용들을 올리도록 하겠습니다. :)