programming python

요즘 보고 있는책. 누구지.. 검색하다가 찾은 블로그 하나 있었는데.. 그 분은 완전 전산 일 하는거 같던데. 그 분이 단계별 추천한 책 중 중급(?)에 속하는 책이였던 걸로 기억한다. 예전에 빠르게 활용하는 python3 프로그래밍인가 하는 책을 한번 훑은적이 있었는데.. 그래서 이 책도 그정도 겠지란 생각에 안볼려가다 결국 본다. 사실 원서라 좀 꺼려지기도 했고.. 워낙 내용도 방대해서 너무 양도 많고.. 근데 결국은 봐야 겠단 생각에 하루에 한 sub-chapter만 읽어보자는 생각으로 읽기 시작했는데. 뭐 몇일 되지도 않았다. 이제 겨우 part 1봤나. sneak preview 라고 하는것 까지 봤다. sneak preview에서는 이 책의 전체적인 내용을 특정 mission을 가지고 간략하게 설명한다. 무슨 내용있는 책인지 알려준다. 음.. 꼭 다 읽어봐야 겠다는 생각이 든다.


책 좋은 거 같음.




<기억해야 할 것>


1.script 상단에 #!를 사용하여 direct하게 executable 하게 만들 때 python path를 적는 것보다 unix의 env를 이용하여 #!/usr/bin/env python이라고 하는게 좀 더 general 하다.
2.os.folk를 쓰지 않고 os.system()에서 실행문 뒤에 &를 붙이게 되면 parallel하게 작업을 수행할 수 있다.
3.spawned process는 모 process에서 환경변수를 상속 받기 때문에 os.environ의 환경변수를  process 간에 communication의 변수로 사용 가능하다. 이 변경된 환경 변수는 top-level process의 상위 레벨에서는 영향을 미치지 않는다. 곧 system shell에는영향이 없다. 
4.for line in open('something') 의 방식을 추천. 원래는 난 readline()을 많이 사용했는데. 같은 의미란다. 차이는 파일을 다읽었을때 return 하는 값. readline은 ''을 iterator는 exception을. 
5. 왜 built-in function인 open의 두번째 argument가 b이냐 아니냐가 unix 계열에서 다르지 않는지 설명한다.
6.built-in open 함수도 있지만 os.open 함수도 있다. 왜 있냐? low-level의 처리가 가능(다양한 mode; O_EXCL을 사용하면 synchronization도 가능).
7.parallel processing 방법으로 process forks 와 spawned thread가 있다. 근데 thread가 platform에 영향을 안받는다.
8.os.execlp를 써서 성공적으로 프로그램을 실행시키면 return이 없다. 거의 program replace라고 생각하면 됨.

package in python

프로그램 짜면서 package화를 잘 안시켜봤다.. 사실.. 내가 짜는 프로그램이라고 해봤자.. 조각난 script 정도라서.. 여튼 package의 개념이 좀 없었는데.. 이번 De bruijn graph package보면서 함 봐봤다. __init__.py의 파일의 개념이 좀 이해가 안갔는데 첨에볼때.. 어떤 책에서는 마치 import * 할 때 package안의 어떤 모듈들이 import 되는지에 대한 list를 명시 하는 것처럼 설명하기도 하고.. 그런데 어떤 package를 보면 아예 __all__ 이란 내장 변수는 없고 그냥 코드만 있기도 하다. 

<__init__.py>

이런 개념을 젤 잘 설명해주는 것이 여기. package도 모듈과 같이 사용될 수 있는데 그 걸 가능하게 해주는 파일이 __init__.py 


<setuptools>
그리고 setuptools 이용하는 방법 : 여기

<eggtestinfo>
아.. 그리고.. 하나 더 이것땜에 이게 뭔가 찾느라 고생 많이 햇는데.. 
전에 de bruijn graph 포스팅 한거에서 파일 받아서 보면 setup.py 에 test_suite이라는 인자(?)가 있는데 이는 아무리 뒤져봐도 setuptools 패키지 설명에서는 찾아볼수가 없었다. 이건 여기에 설명. 




<python egg>
아 이렇게 된 이상 python의 egg를 안 찾아 볼수가 없는데(항상 미뤄왔던 일..)
여기에 잘 나온다. 
요약하자면..
python egg는 리눅스의 rpm 같은 것으로(perl의 cpan에 더 가까워 보인다) easy_install 이라는 프로그램(easy_install은 setuptool package의 한 부분으로 setuptools는 distutils를 기반으로 한다)으로 쉽게 install이 가능케 한다. easy_install 라고 하면(easy_install 도 가능)  알아서 pypi에서 그 package가 있는지 검색해서 install 해준다(구글에서 검색하면 pypi에 어떻게 package를 등록하는지도 나온다). 또 좋은건 dependency 해결해 준단다. 
아.. classifiers 라는게 egg를 만든담에 pypi에 올렸을 때 제대로 된 카테고리에 올라가게 끔 하는거구나.. 
글구 pythonpaste에 대해서도 언급하는데 이건 다음에..

RNA-seq 에 관한 이것 저것

<expression unit> 
RPKM : reads per kilobase of transcript per million of sequenced read
FPKM : expected fragments per kilobase of transcript per million fragments sequenced


음.. 정확한 개념은 ERANGE 논문과 Cufflinks 논문을 참조하면 된다. 


간단한 intro는 Cufflinks 논문 intro에 잘 나와있는데.. RPKM이라는 개념이 gene 단위 레벨에서의 expression unit으로 생각한다면 FPKM은 transcript 단위에서의(하나의 gene에서도 여러 isoform의 transcript가 존재 가능) expression unit으로 생각하면 된다.


내가 여기서 이해가 안갔던 단어가 fragment인데.. 이건 뭐냐면 RPKM은 read 하나하나에 의미를 두는데 반해 FPKM은 paired-read의 한 쌍이 fragment를 의미한다 (별거 아닌건데..). 


자세한 내용은 cufflinks 논문의 supplementary 데이터에 나오는데.. 이해 되면 다시 정리 하도록 하련다.




참고로 KOBIC에서도 NEUMA 라고 NAR에 논문을 낸게 있는데 이것도 보면 좋을듯. 이것 역시 정리되면 다시 올린기로 한다.


<bias> 
EMBL의 자료를 보면 strand를 정보를 살리기 위한 두가지 방법과 그냥 strand를 안살리고 한거랑 tiling array로 비교한 그림이 있는데 bias가 상당하다. random hexamer bias도 있고.. 관련 자료를 좀 훑어봐야 할듯.

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를 보지 않는 이상.