bowtie option

bowtie aligner option을 번역해 놓는다. 한번 훑을 겸 해서.


The bowtie aligner

• -v/-n/-e/-l (-I/-X/--fr/--rf/--ff 옵션은 paired end를 위한 것)  옵션:  alignment의 기준이 되는 옵션
• -k/-a/-m/-M/--best/--strata는 위 기준을 통과한 alignment에서 report을 얼마만큼 할것인지를 정하는 옵션

예로 이해하는 option

• -n 2 -l 28 -e 70 : read의 왼쪽에서 28bp까지 2까지 mismatch가 허용되고 mismatch 된 bp의 phred quality score의 합이 70을 넘지 말아야 한다(n은 seed에서의 mismatch 갯수, l은 seed의 길이). 만약 위와 같은 케이스가 많을 경우 mismatch가 적은 것이 첫번째 기준이 되고 두번째는 mismatch의 phred score 값이 적은 것이 기준이 된다.
• -v 는 quality value는 무시한채로 전체 read에서의 mismatch 허용 갯수를 지정하는 옵션. -n과는 상호 배타적.
• ./bowtie -a -v 2 e_coli --suppress 1,5,6,7 -c ATCGCGA : -v옵션으로 read 전체에서 2개까지의 mismatch를 허용했고  --suppress 옵션으로 output의 1,5,6,7 컬럼을 나타내지 못하도록 했고 -a 옵션으로 -v 기준을 통과하는 모든 output을 report 하게 한다.
• -k 3 : 결과에서 아무거나 3개를 출력하게 한다 default 가 -k 1.
• -a --best : 전체 결과를 출력하는데 위에서 말한 기준대로 alignment 좋은 순서대로 출력
• --strata 는 -n 옵션일 때 seed에서의 mismatch가 기준이 되고 -v일때는 전체 read에서 mismatch가 기준이 되서 계층을 나눈다. 
• -a --best --strata 하면 전체에서 alignment 좋은 순서대로 나열해서 첫번째 계층에 있는 alignment 만 출력. strata를 쓸려면 best가 명시 되어야 한다.
• -a -m 3 : 출력이 3개를 넘어가면 그 결과는 출력되지 않는다. -m 은 unique를 찾을때 쓰기 좋다.
• -a -m 1 --best --strata : 전체 결과를 좋은 순서대로 나열한 뒤에 첫번째 계층의 결과만 출력하는데 그 수가 1을 넘지 않아야 한다. weaker form of uniqueness
• -a -m 1 : strong form of uniqueness
• -I 60 -X 100 --fr 
• paired-end alignment 에서는 --strata랑 --best는 허용되지 않는다.

bowtie -I 157 -X 257 --fr -p 3 --chunkmbs 500 -m 10 --max maxRead --un unmappedRead  --solexa1.3-quals ../total_cdHit/ref/total_cdHit_1.0 -1 result_trim_1.filtered -2 result_trim_2.filtered > bowtie.result

# paired end data를 돌릴때 명령어. -m을 10으로 했다는 건 align했을때 candidate 가 10개 넘어가면 그 read는 maxRead 파일에다가 기록하고 mapping이 안된 read는 unmappedRead 파일에 기록한다. paired end mapping 시 insert size는 최소 157보다 크고 257 보다는 작은 것만 추린다. thread는 3개를 만들어서 parallel 하게 돌리고 각 thread마다 500mb 의 메모리를 할당한다. 

RNA-Seq denovo assembly pipeline 1

problem
RNA-Seq 데이터를 de novo assembly 할려고 하는데 가장 최근에 나온 것이 trinity. 사실 이게 돌아가면 끝나는 문제이긴 한데.. 이게 의외로 시간도 많이 걸리고 안돌아 가는 경우가 생긴다(그 원인에 대해서는 아직 조사하지 않았음). 


solution
그래서 차선책으로 생각한 것이 다음 논문의 "additive Multiple-k" 방법. false positive가 존재할 확률에 대한 정확한 estimation이 없기 때문에 잘못된 transcript를 assembly 할 경우가 생기겠지만 일단은 sensitity를 높이는데 목적을 둔다.


procedure
내 붕어랑 동급인지라 기록을 안해놓으면 다 까먹어서 명령어만 쭉 나열해 본다. 중간 중간 in-house script(python으로 시작하는 건 전부)가 존재하고 프로그램들의 option은 아직 최적화 되어 있지 않다. draft라 생각하면 되겠다. 시간되면 이어놔야지.

1.  perl popoolation_1.017/basic-pipeline/trim-fastq.pl --input1 Pum2_1.fastq --input2 Pum2_2.fastq --output result_trim #velvet을 사용할 것이기 때문에 quality trim의 영향을 많이 받는다. 해서 velvet을 하기 전에 popoolation 패키지의 프로그램을 사용해서 trimming 한다.
2. python /home/bi/code/mergeWOmemory4velvet.py result_trim_1 result_trim_2 Pum2_combine.fa 71 #velvet의 input파일을 생성하기 위해 paired end 데이터를 합치는데 multiple-k 방법을 사용하기 때문에 trimmed read의 길이가 최대 k의 길이보다 긴 것들만 골라낸다.
3. velveth Pum2_velveth 41,71,4 -shortPaired Pum2_combine.fasta 
4. python runVelvet.py Pum2 206 #velvetg를 돌리는 batch file. insert size 만 고려했다.
5. python combineContigs.py Pum2 #velvet을 위처럼 돌리면 여러 폴더고 각각의 k-mer 사이즈마다의 contigs가 생기는데 이를 하나의 multifasta 합치는 것.
6. cd-hit-est -i contigs.fa -o result_cdHit -c 0.95 -n 8
7. bowtie #read를 cd-hit으로 만든 cluster의 ref seq에다가 mapping 한다. 이는 RPKM을 구하기 위한 것
8. python stats.py #assembly 결과 정리

2011 Asian Institute in Statistical Genetics and Genomics

오늘부터 한국 유전체 학회에서 진행하고 있는 워크샵을 다녀왔다. 내가 들은 건 임상 유전학. 결론부터 말하면 만족스러운 강의였다. 


발표자는 성균관의대 기창석 박사님, 서울의대 성문우 박사님.


기창석 박사님의 강의는 멘델의 유전법칙에서 부터 가계도를 그리고 보는 방법, genetic disease의 타입과 어떻게 질환이 나타났을때 genetics disease인지 판단하는지 genetics disease의 detection 하는 실험적인 방법과 genetics disease의 진단과 관련 유전자의 탐색의 실례까지. 


서울의대 성문우 박사님의 경우 임상유전학에서의 사용되는 기본 용어와 각종 실험 방법, 관련 데이터 베이스, risk estimation, risk model 과 genetic counseling 에 대한 이야기까지.


그럼 과연 무엇이 만족스러웠느냐? 
만족을 했다는 것은 기대가 충족되었다는 것이니까.. 내가 어떤 생각으로 NGS 분석 기법에 대한 강의가 아닌 임상 유전학을 선택했는가에 만족의 기준이 있다. 
난 bioinformatician(물론 아직까지 자랑스럽게 이야기 할 정도는 아니지만) 이라는 직업으로 생물학적 데이터를 다루기는 하나 실제 그 데이터가 사용되는 최전선의 일보다는 데이터 프로세싱, 관련 툴 셋업, 툴의 이용, 툴이 돌아가는 알고리즘에 초점이 맞춰서 분석 과정을 선택하고 진행하고 있다.  생물학적으로 이야기 해서는 SNP calling, DEG, DMR 찾기 정도라고 할까. 그럼 이 프로세싱 결과, 내가 찾은 SNP, DEG는 실질적으로 어떻게 쓰일까? 그리고 내가 실제 환자들을 접하고 genetic disease라고 판단했을 때 어떤식으로 실험을 계획할 수 있을까? 란 의문이 많이 들었다. 그 의문의 이유는 결국 최종의 의미를 갖는 일은(물론 모든 일이 의미가 있긴하다) 실질적인 데이터의 적용 단계라는 생각에서였다. 다시 말하면 세부사항이 아닌 좀 더 큰 그림으로 하나의 완성된 결과물을 보고 싶었다.
이번 강의가 바로 그런 생각을 충분히 만족시킬만한 내용들로 꾸며져 있었다. 특히나 기창석 박사님께서 마지막에 들은 실례들은 그러한 결과물들을 느끼기에 충분한 내용들이 였으며, 성문우 박사님의 마지막 risk estimation 부분은 응용사례분만 아니라 좀더 원론적인(내가 말하는 원론이란 실제 판단에 대한 근거에 대한 설명이랄까) 내용을 커버함으로써 만족도를 높였다.


근데 들으면서 계속 꺼림칙 생각이 떠나지 않았는데..  


무슨 생각이 들었냐?
아무리 sequencing이 싸져서 bioinformatics가 붐을 이룬다 해도 결국 대학 병원 연구소의 테크니션 정도의 취급 밖에 되지 않을 수도 있겠구나 라는 생각이 들었다. 결국 환자를 대할 수 있는 것은 의사이고.. bioinformatics는 어쩌면 그저 pcr 정도의 일이 될 수 있겠구나라는 생각. 


그래서 어쩌란 거냐?
두가지다. 

1. 지금과 같은 비슷한 일을 한다고 가정한다면.. 의사랑 동등한 위치가 되는 거다. 이런 표현하는 걸 상당히 싫어하지만.. 학위나 스펙으로다가 나도 책임을 질 수 있는 위치가 되는 거다.
2. 다른 하나는 google 같은 회사를 세우거나 그런 곳에서 일하는 것. 뭔 말인고 하니 google이 쌓여가는 웹 문서를 기반으로 성공한 기업처럼, 쌓여가는 유전체(일단은 내가 있는 업계때문에 생물학적 데이터를 유전체에만 국한한다)를 활용하여 비지니스를 만들 수 있는 회사에서 일하는 것을 말하는 것이다. 유전체가 단순 의료쪽에만 적용 될 것은 결코 아닐 것이라 생각한다. 여러가지 응용 사업이 생길 수 있고 특히나 오늘 강의를 들으면서 문뜩 생각한 것이 결혼 정보업체에서 회원들의 유전체 정보를 바탕으로 짝이 될 수 있는 candidate를 뽑는 것으로 한 예가 될 수 있겠다(기준은 단순하게만 생각한다면 결혼시 자식이 유전적 질환이 적게 걸릴 확률) .

결론은 뭐 싱겁지만.. 강의 맘에 들엇다.


아 그리고 한가지..
왜 굳이 genetic test를 통해서 disease를 규명해야 하는 것인가?
이거에 대한 박사님의 대답이 참으로 인상 깊었는데.. 내가 강의를 듣기 전에 갖고 있던 생각은 막연했다.  설령 치료가 안되는 유전 질환이라도 원인을 알고 있는 것이랑 모르고 있는 것이랑은 potential이 다르다. 이것이 고작 내가 할 수 있는 대답이였다. 박사님은 뭐라했냐? 
2가지를 이야기 했다. 하나는 의외로 한 종류 일 것이라고 생각한 유전 질환이 알고 보면 여러 subtype이 있을 수 있다. 그래서 유전적으로 정확한 원인 유전자를 규명하는 것은 매우 중요한 일이다. 다른 하나는 환자뿐만 아니라 환자가 속한 가족의 삶의 질이 달려있다. 환자의 유전질환이 autosomal dominant 냐 autosomal recessive냐를 확인함으로써 다른 가족들의 disease 발병의 확률 정도를 알 수 있고 그들의 불확신에서 오는 불안을 해소할 수 있다. 
명쾌하다.

STR (short tandem repeat)

갑자기 뜬금없이 STR에 대해 조사하게 되었다. 원래 모르던 거니 알아보면서 정리 해보자.


STR이 뭐냐? 위키에 따르면 두개 이상의 nucleotide가 연속적으로 반복되어 있는것. 보통 non-coding intron에서 나타난단다. STRP (STR Polymorphism)은 사람마다 STR loci에 repeat의 수가 다름으로서 생기는 다형성. 법의학적으로 굉장히 많이 사용된단다. 


여기까지만 하고 DB부터 검색해 보자.
구글부터 검색을 해보자면.. 여러개 나온다. ATCC, YHRD, ENFSI, STRBase. STRBase가 웹페이지는 허접해 보이나 왠지 좀 academic 해 보이는 관계로 이걸 위주로 보겠다.


STRBase의 DB논문은 2001년에 NAR에 나왔다. 이 DB는 NIST(National Institute of Standards and Technology)에서 관리하는 DB로 뭐 STR DNA marker, 관련 논문 정리, 관련 연구 기관, STR multiplex kits 등 STR 과 관련된 정보는 다 모아놓은 DB라고 할 수 있겠다. 여기까지가 abstract 내용이고 


STRBase에서 제공하는 intro를 보자면 forensic DNA typing의 역사서 부터 시작해서 어떤 케이스에 STR이 사용되는지 실험kit, 기계, base call signal, 등등 아주 간략하게 그림만 나와있는데 마지막 슬라이드가 흥미롭다. CODIS(COmbined Dna Index System)라고 FBI의 DB가 있단다. 아래 그림에 나와있는 13 개의 dna loci 에 대해 법의학적으로 DNA profiling을 한다고 한다. 

아 근데 이게 좀.. 음 STR이 microsatellite의 하나인 것 같은데.. 위키에서도 보면 두 개의 정의가 merge하도록 suggestion 되어 있고 다음 pdf를 봐도 microsatellite = STS + STR로 되어 있다. 음.. 나와 같은 고민을 한사람이 있었다. 여기 참조뭐.. 여튼


여기까지가 대충 한번 훑어 본거고. 그럼 내가 알아봐야 하는건 무엇인가. 이 STR의 길이의 variation과 maximum length. 음 근데 이거 STRBase에 잘 나와있다. 특히나 cat, dog, cattle의 경우 굉장히 편리하게 table로 정리해놨다. 음.. 사이즈가 대략 maximum 몇백 하는군..


하지만 내가 궁극적으로 찾아야 하는건 뭐냐. 식물에서의 STR. googling 하자마자 나온게 재밌게도 대마초의 STR에 관한 논문. STR이 법의학적으로 사용이 많이 되서 그런지 내가 찾아보면서 본 논문들은 거의 Forensic Science International 에서 나온 듯. abstract를 보자면 93개의 대마를 5개의 STR marker 로 profiling 했단다. 이 5 loci 에서 총 79 개의 allele이 확인됐다.  그담에 AMOVA로 genetic variation을 찾고 fibre crop accession과 drug accession의 variation이 어떻고 뭐하고 뭐한다는데.. 뭐래는 건지 모르겠다(결론은 drug종과 fibre종간의 차이를 확실하게 구분짓지는 못하나 drug 종간에는 사용할만하다 뭐 이런거 같은). 여튼 확인하고자 하는 건 STR의 maximum length. 아래 표가 그것을 의미하는게 아닐까 한다. 얼추 2~300 정도 하는거 같네.

그 다음에 드는 생각이 대마도 있는데 arabidopsis는 없을까? 역시 없을리 없다. 다음 논문을 보자(그닥 최신건 아닌데.. 아마도 다른 논문이 더 있지 않을까 싶지만..). 


어이쿠야.. 인삼에 관한 것도 있다. 중국애들이 낸 논문인데 2개의 microsatellite loci(CT 12, CA 33) 에 대해 동양의 ginseng과 미국거랑 different allele pattern을 보인단다. 뭐 여튼 난 알고 싶은건 길이니까. 아래 표를 참고하자.

gene fusion tools (DeFuse)

RNA-Seq 분석의 일환으로 요즘 특히나 gene fusion을 많이 찾기 때문에 이를 또 넘어갈 수 없기에 일단을 프로그램을 돌려본다는 생각으로 셋팅한다. 

구글에 검색하면 바로 나오는 것이 여기. gene fusion을 찾는 각종 툴을 모아놓은 것.

여기서 난 아무래도 가장 최신에 논문이 나온 DeFuse를 이용하려 한다.

일단은 논문도 안보고 실행에 초점을 맞춰서 manual을 정리해 본다.


<Prerequisite Setting>
DeFuse를 사용하기 위해서 선행적으로 셋팅되어 있어야 하는 것이 bowtie, blat, faToTwoBit, R의 kernlab library.  그리고 각종 데이터를 받고 나서 DeFuse의 scripts 폴더 안에 config.txt 파일의 내용을 변경해줘야 한다(프로그램과 데이터의 path를 지정). 이 config.txt 파일내의 모든 경로는 fully 절대 경로여야 한다(뭐 꼭 그런건 아닌데 여튼 이게 편하다). 그리고 bin을 경로를 적을시에는 실행 파일까지 적어야 한다(/usr/local/bin/R 이런식으로). 그리고 절대 라인 끝에 space가 들어가면 안된다(아.. 이게 완전 피곤한 건데.. 왜 에러가 나는지 몰랐는데 얘네가 짜놓은 script가 config.txt를 파싱해서 정보를 취하는데 line 마다 strip을 하지 않고 그냥 바로 load 하기에 에러가 난다).


그런 담에 create_reference_dataset.pl 을 실행시키는데(아마도 bowtie 를 이용하는데 genome index 파일을 만들고 뭐 이런 저런 일을 하는 script인 갑다)  에러가 났을시에 config.txt 에 명시되어 있는 dataset_directory 폴더에 생성된 파일들을 확인해본다. 여차하면 reference.fa.index도 다시 지워야 한다 (무슨 무슨 파일이 생성되는지 잘 확인 요망. 에러시 trackback 하면서 고쳐나가야 해서).

maximum Matching in bipartite graph

cufflinks를 보다보면 이런 말이 나온다. "It does so by reducing the comparative assembly problem to a problem in maximum matching in bipartite graphs".

아.. comparative genome assembly가 무엇이냐. 생소한 표현이다 assembly면 de novo assembly지 뭔 comparative냐. 아.. 근데 문득 드는 생각이.. 혹시 이거 mapping 해서 assembly 하는거니까 그걸 comparative genome assembly 라고 하는거 아니야? 음.. 맞는거 같다. 특히나 이 사이트를 보니까 더욱 확실한거 같다. comparative genome assembly 이거 그냥 mapping assembly를 뜻하는 거다. 근데 찾다 보니 이런 논문도 나오더라. 역시나 이 논문도 CBCB에서 나온거다. 이 논문에 완전 정확하게 comparative genome assembly의 정의를 내린다. AMOS 프로젝트에서 나온 program을 소개하는 논문인데 이건 패스하자. 이게 중요한게 아니다.

그리고 또 bipartite graph는 뭐냐? 이건 쉽게 끝날 문제가 아닌듯 보인다.. 일단은 여기를 보자. 인트로다. 그거 다 봤으면 여기를 보자. 이거 다 보고 정리하겠다. (단순 graph theory에 대한 건 여기를 보자)


The marriage Problem and Matchings
결혼이란 문제로 maximum matching problem을 설명하자면 ..
아래 그림과 같이 각각 n 명(여기서는 6명)의 싱글 남녀가 있다고 할 때 점은 남여 개체를 의미하고 선은 상호간의 관심이 있는 상대를 의미한다. 예를 들어 Ann 이라는 여성은 Adam과 Bob에게 관심이 있다. 이 때 match가 됐다는 것(matching)은 각 개체가 오직 한사람의 상대와 이어짐을 의미 한다. 아래 빨간 선이 matching을 의미 한다(물론 다른 경우의 수도 있다). 특히 이때 모든 개체가 빠짐 없이 matching이 되면 이를 perfect matching이라 한다(아래 그림). 

그러나 일반적인 문제에서는 위와 같이 perfect matching 이 어려운데 이 때 최대의 match를 끌어내는 것이 maximum matching problem이다.


Hall's Theorem
Hall 이란 사람이 perfect matching 이 되기 위한 조건을 정리 했는데 다음과 같다
Hall's Theorem : 두 클래스 X,Y를 갖는 bipartite graph G=(V,E) 가 주어졌을 때 이 G가 perfect matching이 되기 위해서는 클래스 X의 subset인 S에 대해 |Adj(S)| >= |S|를 만족해야 한다. Adj(S) 란건 S에서 부터 연결된 클래스 Y의 개체를 의미. 다르게 이야기 하면 X의 subset S랑 연결 선이 있는 클래스 Y의 사람 수가 최소한 S의 사람수 이상은 되야 된다는 이야기.


그런데 이건 정리지 뭐 이거 자체가 알고리즘이거나 하진 않다.  그래서 여기서 이야기 하길 Dijkstra's algorithm에서 효과적인 알고리즘은 찾는단다. 단계적으로 matching을 개선해서 결국 optimal solution을 찾는다는 건데.. 어쩌라는 거냐. 다음 예를 보자.


아래와 같은 상황이 있다고 하자. 안타깝게도 Dorothy, Evelyn, Carl, Erik은 짝이 없다. 이거 어떻게 해결할거냐. 일단 Dorothy 부터 보자. Dorothy에서 연결될 수 있는 짝이 Bob이고 그렇게 되면 Bob에 연결된 Ann을 Adam에 match 시켜야 하고 다시 또 Adam에 연결되었던 Beth는 Carl과 연결시킨다. 이렇게 함으로써 기존의 match M에서 새로운 match M' (|M'| = |M| +1) 이 생기게 된다.


여기서 중요한건 alternative path (= a-path)가 있었다는 것.  음 이게 뭐냐면.. unmatch vertex 에서 시작해서 unmatch vertex로의 path. 즉 길. 근데 왜 alternative이냐 이 길에서 홀수 번째 edge는 현재 match M에서 사용되지 않고 있는 edge이고 짝수 번째 edge는 M에 사용되는 edge인데 이 짝수 번째 edge 대신 홀수 번째 edge를 match로 여기면 이게 새로운 match M'가 되기 때문이다.

이것에 대한 pseudo code는 다음과 같다.

여기서 핵심은 3번째 라인 a-path를 찾는것. 여기서는 modified Breadth-First-Search를 언급하나 이건 패스.


마지막에 a-path가 없을때까지 improve를 시켰다면 이것이 과연 optimal maximum match 인가를 증명한다. 이건 생략

using k-mer distribution in de novo assembly

assembly를 하기 전에 k-mer의 depth별 frequency로 genome에 대한 정보를 얻을 수 있는데 그에 대해 간략하게 이야기 하려한다.

Genome assembly 를 하기 전에 순수 read를 이용하여 genome의 대략적인 size를 예측할 수 있다. 전체 read를 특정 k-mer size (bp)로 hashing을 해서 각 k-mer의 depth 별 frequency (=depth)를 계산했을 때, 가장 frequency가 높게 나온 depth (M), read length (L), k-mer length (K), sequencing depth (N) 간에는 다음과 같은 관계가 성립하게 된다.

M = N * (L – K + 1) / L

여기서 K와 L은 주어진 값이고 M은 계산으로 아래와 같은 plot으로 계산이 가능하기 때문에 N, 즉 sequencing depth를 구할 수 있게 된다. 이때 구해진 sequencing depth로 전체 read의 length로 나누어 주게 되면 assembly를 하고 있는 genome의 대략적인 genome size를 예측할 수 있는 것이다. 

아래 그림이 그 예. BGI의 giant panda 논문에서 발췌

또한 위의 k-mer의 depth 별 frequency를 이용하여 genome의 heterozygosity를 판단할 수 있다. 아래 그림에서 보듯이 k-mer의 depth 별 frequency에 대해 plotting을 했을 때 위의 그림과 같이 하나의 peak가 아닌 2개의 peak가 보일 시에 genome이 highly heterozygous 하다고 볼 수 있다. 이는 두 haploid가 동일하다면(=homozygosity) 각각의 haploid에서의 k-mer plot이 동일해야 하기 때문에 하나의 peak가 나타나는 형태여야 하지만 두 haploid가 heterozygous 하다면 서로 다른 위치에서 peak가 생성되기 때문에 bi-modal distribution이 생성되게 된다.
아래 그림(potato genome sequencing 논문에서 발췌)이 그 예로 RH의 경우 higher heterozygous 함을 알 수 있다.