간만에 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 분석 해보기 -Align-

이번 글은 Align하는 명령어 편! 되겠습니다.

HISAT2

$ ~/hisat2/hisat2 -p num_thread -x /path/to/reference_index --known-splicesite-infile /path/to/splicesite_file --tmo -1 /path/to/pair_1.fastq.gz -2 /path/to/pair_2.fastq.gz | samtools view -Sb - | samtools sort -m memory -@ num_thread - output.sortByCoord
혹은
$~/hisat2/hisat2 -p num_thread -x /path/to/reference_index --known-splicesite-infile /path/to/splicesite_file --dta
-1 /path/to/pair_1.fastq.gz -2 /path/to/pair_2.fastq.gz | samtools view -Sb - | samtools sort -m memory -@ num_thread - output.sortByCoord
혹은
$~/hisat2/hisat2 -p num_thread -x /path/to/reference_index --known-splicesite-infile /path/to/splicesite_file --dta-cufflinks -1 /path/to/pair_1.fastq.gz -2 /path/to/pair_2.fastq.gz | samtools view -Sb - | samtools sort -m memory -@ num_thread - output.sortByCoord

--tmo : Known transcript에 대해서만...

--dta : StringTie를 사용할 경우...

--dta-cufflinks : cufflinks를 사용할 경우...
인것같은 설명이 되어 있는 옵션입니다.

자세한 설명은 메뉴얼 사이트로... (일단 3가지 다 돌리는 걸로)

STAR

$ ~/STAR/STAR --runMode alignReads --genomeDir /path/to/reference_idx --runThreadN num_thread --outSAMtype BAM SortedByCoordinate --outWigType wiggle read1_5p --sjdbGTFfile /path/to/gtf --outBAMsortingThreadN num_thread --readFilesCommand gunzip -c --readFilesIn /path/to/pair_1.fq.gz /path/to/pair_2.fq.gz --outFileNamePrefix /path/to/output/prefix_ --quantMode TranscriptomeSAM --limitBAMsortRAM num_memory

자세한 설명을 원하시면 이곳으로!!

Kallisto

~/kallisto/kallisto quant -i reference_idx -o /path/to/output_dir --bias -t num_thread /path/to/pair_1.fq.gz /path/to/pair_2.fq.gz

Kallisto 메뉴얼은 여기로...


Salmon

~/salmon/salmon quant -i reference_idx -p num_thread -l A -o /path/to/output_dir --validateMappings --seqBias -1 /path/to/pair_1.fq.gz -2 /path/to/pair_2.fq.gz

보다 자세한 설명은 그곳으로..

Align에 필요한 명령어는 이정도면....  현재 명령어에는 각 tools에서 strand specific reads에 대한 옵션은 생략되어 있으니 library를 strand specific으로 했으면 해당 옵션을 활성화 하면됩니다. :) 저는 strand specific한 데이터가 아니라서... ㅎㅎ

출처: JYP

간만에 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

간만에 RNAseq 분석 해보기 -설치편-

근 몇년 만인가요..
간만에 RNAseq을 손에 뭍혀봅니다.!! :)

여러가지 Tool을 조합해서 일단 해보고 좋은 놈 골라 쓸 생각이라서
여러 aligner와 abundance tool을 사용할 예정입니다.
그래서 오늘은 서버에 설치된 프로그램 유무 확인 및 업데이트를 해봐야겠습니다. Orz

그래서 오늘 설치편되겠습니다.

오랜만에 작업하는것이라 잘 될지 모르겠습니다. ㅋ

사실 이건 그냥 스킵하셔도 됩니다.
그냥 제 작업 로그 남기는 거거든요 ㅎㅎ :)

bwa (https://github.com/lh3/bwa, 2019.02.27 기준 0.7.17)

$ git clone https://github.com/lh3/bwa.git
$ cd bwa
$ make

hisat2 (https://ccb.jhu.edu/software/hisat2/index.shtml, 2019.02.27 기준 2.1.0)

$ wget http://ccb.jhu.edu/software/hisat2/dl/hisat2-2.1.0-Linux_x86_64.zip
$ unzip isat2-2.1.0-Linux_x86_64.zip

STAR (https://github.com/alexdobin/STAR, 때마침 2019.02.25일에 업데이트 ㄷㄷ)

$ wget https://github.com/alexdobin/STAR/archive/2.7.0e.tar.gz
$ tar -xzf 2.7.0e.tar.gz
$ cd STAR-2.7.0e

StringTie (https://ccb.jhu.edu/software/stringtie/ 2019.02.27 기준 1.3.5)

$ wget http://ccb.jhu.edu/software/stringtie/dl/stringtie-1.3.5.Linux_x86_64.tar.gz
$ tar zxf stringtie-1.3.5.Linux_x86_64.tar.gz
$ cd stringtie-1.3.5.Linux_x86_64

HT-Seq (https://github.com/simon-anders/htseq, 2019.02.27 기준 0.11.0)

$ pip install HTSeq

혹은

$ pip install 'HTSeq==0.11.0'

혹은

$ easy_install 'HTSeq==0.11.0'

잘 안되시면 주위에 리눅스에 HTSeq를 설치할 줄 아는 친구를 찾으십시요! :)

Kallisto (https://pachterlab.github.io/kallisto/, 2019.02.27 기준 0.45.0)

$ wget https://github.com/pachterlab/kallisto/releases/download/v0.45.0/kallisto_linux-v0.45.0.tar.gz
$ tar zxf kallisto_linux-v0.45.0.tar.gz
$ cd kallisto_linux-v0.45.0

Salmon (https://combine-lab.github.io/salmon/, 2019.03.04 기준 0.12.0)

$ wget https://github.com/COMBINE-lab/salmon/releases/download/v0.12.0/salmon-0.12.0_linux_x86_64.tar.gz
$ tar zxf salmon-0.12.0_linux_x86_64.tar.gz
$ cd salmon-0.12.0_linux_x86_64/bin

일단 설치는 여기까지,
옛날에는 컴파일하느라 설치가 반이었는데
요즘은 걍 컴파일이 잘 되어 있어서 ㅎㅎ :)

다음 포스팅은 reference 생성 내용입니다. ㅋ
-별거 없습니다. :)

출처: JYP

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

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

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

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

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

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

>이곳<

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

출처: SM

Circos를 그까이꺼 한번 해보자 -설치편-

몇년 전부터 Genome을 그린다면 염색체 모양이 아닌
죄다 동그랗게 그리게 한 장본인!!
바로 그녀석! Circos (발음은 나도몰라 그냥 부르고 싶은데로!!)

공식 사이트: http://circos.ca/

사실 circos는 NY Times와 같은 미언론의 인포그래픽의 한 축을 담당했었는데
지금은 나도 모르겠다능.. 여튼 그래서 circos가 기능은 막강하지만
그만큼 사용하기 까다로워서 동그안 genome은 그려주지만 비슷한 기능을 하는
간단한 프로그램들도 많이 나왔다능....
(그래서 결론은 꼭 circos를 사용해야 그릴 수 있는건 아니라는걸 먼저 밝혀둡니다.)

현재 버전은 0.69-6 랍니다.
제가 마지막에 사용했던 circos버전은 0.67-4라는건 비밀...

그래서 일단 다운받아서 실행해 보았습니다. 업데이트 되면서 몬가 더 필요로 하고 있는지..

$ wget http://circos.ca/distribution/circos-0.69-6.tgz
$ tar zxf circos-0.69-6.tgz
$ cd circos-0.69-6/bin
$ ./circos
하면 (정상적으로 perl들이 설치되어 있어도)에러 납니다.
$ ./circos -modules
하면 설치가 안되어 있는 모듈들을 보여줍니다.
시스템 perl이나 circos가 사용하려는 perl에 circos가 필요로 하는 모듈들을 설치해 주면됩니다. :)
circos 설치 방법을 참고하셔서 설치하시면됩니다.

출처: JYP

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