SAM format flag cases

paired end read가 mapping 되었을때,

• 163 (0x80,0x20,0x2,0x1) 과 83 (0x40,0x10,0x2,0x1)
• 99 (0x40,0x20,0x2,0x1) 과 147 (0x80,0x10,0x2,0x1)

paired end의 한쪽만 align 되었을 때

• 73 (0x40,0x8,0x1) 과 133 (0x80,0x4,0x1)
• 89 (0x40,0x10,0x8,0x1) 과 181 (0x80,0x20,0x10,0x4,0x1)

ncbi genome file extension

NCBI의 genome DB에서 데이터 받으면 여러가지 확장자의 파일들이 있는데 그 설명은 다음 사이트에 잘 나와있다.

.asn	genome record in asn.1 format 
.faa	protein sequences in fasta format, text file
.ffn	protein coding portions of the genome segments
.fna	genome fasta sequence
.frn	rna coding portions of the genome segments
.gbk	genome in genbank file format 
.gff	genome features
.ptt	protein table
.rnt	rna table
.rpt	summary report
.val	binary file (genome project?)

.gb     Genbank?
.gpff   Genbank protein

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 하면서 고쳐나가야 해서).