Trinity

논문 소개는 여기서 했고 실질적으로 setting 하면서 기억해야 할거 command line에서 주의 해야 할것 등을 기록하려 한다. 워낙 붕어 머리라. 


일단 trinity download 및 documentation은 여기.


<install>
그냥 tar 랑 gzip 풀고 기본 폴더에서 make 한번 눌러주면 된다. 그럼 inchworm 이랑 chrysalis 가 컴파일 된다(c++로 되어 있음). butterfly 는 java인데 precompiled software라 할 필요 없단다.


<running>
기본 명령어는 sample_data 폴더 안에 runMe.sh 에 있다. run 하는 batch 파일은 기본 폴더에 있는 Trinity.pl 파일이다. 
위 링크에 걸려 있는 documentation을 요약하자면..

• strand-specific sequencing 일때는 -SS_lib_type 옵션을 사용하란다(default 는 strand 안가린다). 
• 또 fungal genome 처럼 gene이 compact 한 genome에서는 UTR 부위에서 transcript의 overlap이 생기기 때문에 -jaccard_clip 옵션을 사용하란다. 요 옵션은 Bowtie가 깔려 있어야 한단다(정확하게는 jaccard_clip의 옵션의 기능을 모르겠다).
• computing grid 일 때는 --run_butterfly 옵션을 사용하지 말고, grid 가 아닐때 단순시 cpu가 많을땐 --run_butterfly 옵션이랑 --cpu 옵션 넣어서 돌려라. --run_butterfly 옵션 없을때 결과 폴더 안에 Chrysalis 폴더 안에 butterfly_commnads 가 있는데 이걸 이용해서 grid를 사용하라.

<hardware requirements>
1G of RAM per 1M pairs of Illumina reads. Chrysalis 단계가 deep recursion 땜시 stack 메모리가 많이 필요하단다. 이럴 때 세팅은 여기처럼 하란다(Q&A에 나온다). 
전체 시간은 1 hour per million pairs of reads(이 표현이 정확이 1Mbp 당인지 아니면 1M paired-end count를 이야기 하는건지 정확히 모르겠네). 


<output file>
--output 옵션에 지정한 대로 아니면 default 인 Trinity_dir 폴더안에 Trinity.fasta 파일이 생기는데 이는 full-length transcript의 fasta 파일.




좀더 자세하게 알고 싶으면 여기,  advanced guide


<Inchworm>
read를 k-mer로 쪼개서 hash table에 저장한다. 이때 base는 2-bit encoding을 사용한다(A:00,G:01 이런식). k-mer는 25 로 고정(package로 돌릴때), 단 따로 Inchworm만 돌리면 32까지는 허용.
monitor.out 이란 파일이 inchworm의 log 파일. 
inchworm.K25.L48.fa 가 inchworm의 결과 파일. K는 k-mer를 뜻하고 L은 contig의 minimum Length를 의미. 그 안에 보면 각 contig의 이름에 >a1;142 K: 25 length: 3697로 되어 있는데 a1뒤의 숫자는 k-mer의 average frequency. 이는 나중에 expression quantity에 사용된다.


<Chrysalis>
Chrysalis는 inchworm에서 생긴 contig들을 alternative spliced transcript 나 paralogous transcript들끼리 clustering 하는 역할. 각 cluster의 contig를 사용해서 de bruijn graph 생성하고 read를 mapping한다. 각 cluster의 de bruijn graph는 butterfly 단계를 위해 각각 다른 파일로 저장됨
이때 생기는 파일은 chrysalis/RawComps.0 폴더 안에 있는데 파일들은 이렇단다.



*chrysalis intermediate outputs*
chrysalis/RawComps.0/comp9.raw.graph  :de Bruijn graph based on Inchworm contigs only
chrysalis/RawComps.0/comp9.raw.fasta :raw RNA-Seq reads that map to the Inchworm contigs based on k-mer composition
*chrysalis outputs to be used by Butterfly as inputs*
chrysalis/RawComps.0/comp9.out :the de Bruijn graph with edge weights incorporating the mapped reads
chrysalis/RawComps.0/comp9.reads :the read sequences and anchor points within the above graph

이것과 함께 butterfly가 돌수 있는 명령문이 있는 butterfly_commands와 butterfly_commands.adj

가 생성되는데 butterfly_commands.adj가 heap size랑 parameter등이 있는 파일로 butterfly 실행시 이 파일을 실행해야 한다.

<Butterfly>
butterfly는 이제 chrysalis에서 생긴 *.out과 *.reads 파일을 보고 가능성 있어 보이는 trascript를 fasta 파일로 report 한다.
Trinity.pl 에다가 --run_butterfly를 줬으면 이건 /util/cmd_process_forker.pl을 돌리는 거다(forker란걸 보니 프로세스 늘려서 실행 시키는 프로그램인 갑다).
결과로 comp*_allProPaths.fasta 라는 파일이 생기는데 이게 possible transcrip의 fasta 파일. accession에 대해선 그냥 메뉴얼보자(쓰기 귀찮다). 
butterfly 돌릴때 그러니까 butterfly_commands.adj에 있는 명령문에서 -V 5 라고 하면 돌아가는 상황이 투명하게 보여지는데, 뭔소리냐면 compacted graph structure를 report 하고 graph 순회하면서 들린 node를 report한다. 그 파일이 comp*_justBeforeFidingPaths.dot. 이 파일은 GraphViz 로 볼 수있단다. 이 파일의 accession도 설명하는데 것도 manual 보자

De Bruijn Graph

genome이건 transcriptome이건 de novo assembly에 사용되는 알고리즘인데..
wiki에 나온건 간단하게 그냥 그 graph의 특징만 나왔다.
그런데 다음 링크는 코드까지 있다. 내 오늘은 그냥 링크만 기록하지만 코드 분석 후에 다시 posting을 하도록 하겠다.


음 코드를 봤는데 생각보다 별게 없다. De bruijn graph를 그리는 것도 아니고.. 링크에 댓글처럼 나도 댓글 달아서 graph까지 그리는 코드를 달라고 해볼까나..


그래서 다른 코드도 찾아봤다.
새로 찾은 블로그는 왠지..


이 코드가 그나마 de bruijn graph를 그리긴 한다. 단.. 거기까지만.. 
visualization 하는 코드도 있는데 xml-rpc인가로 원격으로 하는것 같다.

A Primer on Scientific Programming with Python

사실은 python으로 applet 어떻게 web에 심는건지 알고 싶어서 검색하다가 나온책. applet 심기랑 전혀 상관없다. 과학적 데이터를 어떻게 프로그래밍으로 다루나 뭐 이런 거 나온 책. 수치로 식 세우고 plot 그리고 그런 책. 생각보다 python으로 숫자, 식 관련해서 그래프 그리고 뭐 이런 프로그래밍을 거의 안해봤는데(이런건 R써야 할거 같아서).. 알아두면 괜찮을거 같아서 띄엄 띄엄 보고 있는책. 영건씨가 알려준 어둠의 경로로 받아서 보고 있는데 괜찮은거 같다. 언넝 ipad와서 넣어가지고 다니면서 봤음 하네. 


배우게 된거. 

• 젤 중요한것 하나 : vectorization instead of scalar code. 명시적을 source에 for가 있으면 건 scalar code.
• arrays : list 랑 비슷하긴 한데(list의 variant 라고 한다).. element가 같은 type이고 갯수가 고정일 때 사용. list에 비해 계산이 빠르고 memory도 적게 먹음. 아.. 이게 원래 python의 data type이 아닌데.. array 이란 객체가 있긴 있다. python에. 근데 이거 별로 효율적이지 않단다. 그럼 어딨는거냐? numpy 패키지에 있는 array 쓰란다. from matplotlib.pylab import * 하면 from numpy import * 까지 실행되고 이건 뭐하는 거냐면 matlab 스타일의 syntax를 사용하는 패키지. 
• numpy.array는 b = a[1:-1]하면 copy되는것이 아니라 reference가 달리는것. 곧 b = a[1:-1]로 변수 b를 a 데이터가 가르키는 곳을 가르키므로 b를 변형하면 a도 변형된다.
• 원인을 찾지 못했는데 scitools.std 를 import 하려면 from maplotlib.pylab import *을 먼저 해야 한다. 나만 그런건지 모르겠는데.. 왜 그러지.. 

genome sequencing 정리

genome assembly를 회사에서 시작할려는지 박사님께서 데이터를 하나 보내왔다(구글에서 찾아진다.. 이 말인즉 누구 박사가 어쩌다 하나 찾아서 참고하라고 보냈나싶다).
사실 예전 부터 genome sequencing 부분에서 조금씩 이해가 안가는 부분이 있었는데 이를 확실하게 정리해본다. 첨 human genome sequencing 했을 때 사용됐던 방법이 shotgun sequencing 이랑 bac-to-bac 방식 거기에 대한 소개는 다음에 잘 나와있다. shotgun이야 익숙한데 bac-to-bac은 익숙하지 않다. 특히나 physical mapping 부분은 이해가 되지 않는다(genetics 수업을 제대로 안들었다). 해서 여기를 참조한다(아 요 physical mapping pdf 어렵다.. ).


처음 데이터에 있는 내용은 음.. 각종 시퀀싱 platform을 섞어서 한거 같은데 Bac-by-Bac도 진행해서 굉장히 사용가능한 정보가 많은거 같은데.. ppt 형식만 봐서는 정확히 이해는 가지 않는다. 논문 method를 보지 않는 이상. 

DNAnexus

Galaxy에 이어 DNAnexus도 알아본다. DNAnexus는 사실 회사라 이용할려면 돈을 내야 한다. 듣기로 데이터 10기가에 얼마.. 뭐 이런식으로 돈을 받는다고 하는거 같은데 정확한 정보는 아니다. 

DNAnexus는 Amazon의 cloud computing 을 이용하여 기가 혹은 테라 단위의 large-scale whole-genome sequencing 데이터를 저장 및 분석할 수 있는 web-based platform을 제공하고 있다. 그러니까 사용자는 NGS 데이터와 돈만 있으면 데이터의 저장 및 각종 분석이 web browser에서 클릭으로 수행되어 진다. 편하다. 그럼 어떤 분석들을 제공하냐.. 아래와 같다(아래 내용은 DNAnexus에서 제공하는 whitepaper를 해석한 정도다, white paper는 계정 만들고 나면 conduct new analysis에서 experiment type 저하면 옆에 생긴 link에 있다).

-분석 모듈-

1. RNA-Seq / Transcriptome-based quantification
Annotated transcript의 발현양을 profiling 하기 위한 분석 타입. Annotated trascript의 발현 양을 분석하는 것이기 때문에 novel한 splice site나 novel 3’ end를 찾지는 않는다(이를 위해서는 3번 과 9번 타입의 분석을 이용할 것). Splice junction을 포함한 transcript로 부터의 데이터가 있을 수 있기 때문에 genomic mapping을 사용하지 않고 Refseq이나Ensembl과 같은 데이터베이스에서 얻은 reference transcriptome에 re-align 한다. Reference의 transcript에다가 read를mapping 한 뒤에 RPKM 값을 구함으로서 normalization을 한다. 하나 이상의 trascript를 갖는 유전자의 경우 가장 높은 RPKM 값을 같는 transcript가 report 된다.

2. 3SEQ / Transciptome-based quantification

3SEQ protocol 에 의해 만들어진 library를 분석 할 때 이용하는 분석 타입. 1번 분석 타입과의 차이라면 3SEQ library의 데이터를 이용한다는 것. 일반적인 RNA-seq 은 trascript의 전반에 걸쳐 read가 나오는 반면 3SEQ library는transcript의 3’ end(대부분 3’ UTR) 부분에서 read가 나온다. read 하나가 transcript 하나를 의미. 그렇기에 RPKM을 이용하지 않고 posterior probability를 이용한 weight를 준 read count 값을 이용해서 transcript 양을 측정한다. 샘플 간의 비교를 할 때는 read수로 normalize 하거나 아니면 Z score를 계산해서 비교한다. 

3. 3SEQ / Expressed regions in genome

3SEQ protocol에 의해 만들어진 library를 이용하여 reference annotation에는 없는 novel한 유전자나 alternative 3’UTR을 찾을 때 사용하는 분석 타입. Kernel density estimator방법으로 read의 start site의 smoothed profile을 계산한다. Strand 별로 read의 enrichement peak를 찾고 minimum region enrichment, minimum region length, minimum reads의 3가지 filtering 조건으로 peak를 filtering 한다(4번 타입의 분석 방법과 유사하다).  

4. ChIP-Seq / Peaks or regions

ChIP 데이터를 분석하기 위한 분석 타입. 이 분석에 있어서 Peak랑 region의 개념 차이가 중요한데, peak는 single genomic site에 의한 것이고 region은 single site에 의한 것이 아닌 좀 더 broad 한 영역에 걸쳐진 것. 그래서 promoter binding protein의 ChIP-seq 데이터의 경우 peak-based analysis를 해야 하고 histone binding protein의 경우 regional analysis를 선택해야 한다.

Kernel Density Smoothing : 일단은 read를 groupping해서 read density mapping을 만든다. 이때 이용하는 방식이KDEs(kernel density estimators). Posterior probability가 90% 가 넘는 mapped read들의 mapping position의 중점을 그read의 위치로 생각해서 genome 상의 position마다의 read 개수를 합한다. 그리고 이 density map을 smoothing 하게 하기 위해 KDE를 사용하는데 K, 즉 kernel은 정규분포를 이용하고 h, breadth는 분석 샘플에 따라 다르게 한다(transcription factor의 경우 30, RNA polymerase의 경우 60, histon biding protein의 경우 100, 숫자가 클수록density가 좀 더 부드러워 진다). Regional analysis의 경우 expected extent of the region을 위한 parameter를 제공해야 하는데 이는 나중에 q-value를 구할 때 이용된다.

Combinding forward and reverse strands : ChIP-seq 자체가 read의 끝을 읽는 것이므로 각 strand의 한쪽 끝에서read가 나오기 때문에 원래 하나의 peak가 생겨야 할 곳에 library의 길이의 반 만큼 shift 되서 양쪽에 두개의 peak가 생긴다. 그래서 이 두 개의 peak를 library 길이 반 만큼 shift 해서 합쳐서 원래의 peak를 추론하게 된다.

Experiment versus Background/Input Samples : peak의 KDE값을 normalize 하기 위해 uniformCoverage를 구한다. 이는peak가 없을 때, 곧 총 mappable read가 고루 genome 상에 분포 하였을때의 depth를 구하고 이를 peak의 KDE 값에 나누어 주어 KDE 값을 normalize한다. 만약 input 데이터(immunoprecipitation을 하지 않은 데이터)가 이용 가능 하다면 sample의 normalized KDE/ input의 normalized KDE 값을 이용한다.

False Discovery Rate :

5. Enriched regions in genome

DNase1 hypersensitivity experiment나 open chromatin의 region을 찾는 실험과 같은 특정 실험에서 생성된 read가enrichment된 region을 찾는데 이용되는 분석 타입. 분석 방법은 strand에 상관없이 하나의 density profile을 만든다는 것을 제외하곤 4번 타입의 ChIP-seq 분석방법과 동일하다. 또한 중요 parameter 3가지는 3번째 타입인 3SEQ / region 과 동일하다.

6. Nucleotide-level variation

Sample의 genome과 reference 간에 차이를 확인하기 위한 분석 타입. 차이라고 하는 것은 SNP, MNP, insertion, deletion을 뜻한다.

7. Cancer nucleotide variation

Cancer 관련 variation을 찾기 위한 분석 방법. 6번 분석 방법과 동일한데 두 개의 sample 그러니까 cancer sample과normal sample의 reference에 대한 variation 분석이 동시에 진행된다.

8.Restriction enzyme quantification

Methyl-seq (methylation-sensitive restriction enzymes을 이용하여 methylated 된 DNA를 잘라서 시퀀싱 한 것)과 같이 reference genome에 대한 restriction site를 찾기 위한 분석 타입

9. Alternative splicing, Exome

알려진 exon의 novel splice junction을 찾기 위한 분석 타입. Novel 한 splicing junction을 찾음과 동시에 expression quantification이 가능하다.

10.Population allele frequency

Multiple population에 걸쳐서 일어나는 variation을 찾기 위한 분석 방법. 6번 타입의 분석 방법을 기초로population의 그룹 간의 비교를 가능하게 해준다

Galaxy

삼성의 핸드폰 galaxy말고 large-scale genome analysis의 분석을 위한 platform인 galaxy. 논문은 여기. 사이트는 요기. genome 관련 데이터 분석을 총망라한 사이트. 몇달전에 함 보고 vod 몇개 보다가 그냥 넘겼는데.. 회사 박사님께서 기분이 상하시는 바람에 요약 보고서를 만들어야 한다. 뭐.. 맘에 안드는 점도 있지만.. 회사기 때문에 시키는거 제대로 하는게 내 일이기에.. 글고 뭐 알아봐서 남주는 것도 아니고 이미 많은 논문에서 citation 하기에 알아야 하는 사항인지라.. 즐거운 맘으로 정리 해보자. 대략적으로다가..


사실 홍창범씨(?.. 성함이 맞나 모르겠네.. 이 분 블로그 보면 참 깊다는 생각이 든다. 같이 일하면서 배워보고 싶은 분) 블로그에 잘 나와있는데.. 이것만 봐도 대략 보고서 나온다. 아마도 홍창범씨 블로그 내용이 galaxy screencast의 내용이 아닌가 싶다. 이 VOD에 이미 설명이 충분히 잘 되어 있기에.. 원래는  구구절절하게 내용을 쓸려 했는데 그냥 여기서 접을련다. 나중에 metagenome 분석 부분만 좀 깊게 review를 해볼련다. 끝.


아 하나 더 Galaxy 모듈 중에 큰 항목이 RGenetics가 있는데 이게 뭔지 아직 정확히 모르겠으나 그 tutorial은 여기.

SMASHCommunity

저번 포스팅에 리뷰 했던 논문(enterotypes of human gut microbiome) 쓴 사람들이 만든 metagenome pipeline인 SMASHCommunity. 네이쳐에 나온 논문을 보면.. 아.. 정말 잘했구나라는 생각밖에 안든다(사실 아직 다 읽진 않았다.. 이런 게으름.. 항상 이런식이다.. 상흠 선배 말투가 떠오르는군). 우린 갈길이 멀었구나..




일단 bioinformatics에 나온 논문부터 요약하자면
SMASH는 sanger나 454 read를 분석하기 위한 standalone 프로그램. metagenome pipeline이 없는것은 아니지만(CAMERA, IMG/M, MG-RAST 등) web-server이거나 아니면 단순 read count를 세는 정도라 한다(참고로 네이쳐 논문의 메소드를 보게되면 read count를 그냥 쓰는 것이 아니라 gDNA의 경우 genome size로 나누고 rRNA의 경우 genome 상의 16s rRNA의 copy수로 나눠서 계산한다). 그 밖에 취약점들을 있는데 생략한다.