GFF tools

tophat이나 cufflinks 돌릴때 annotation 파일이 필요한데(뭐 꼭 필요한건 아니다. 다만 있으면 더 좋다는 거지).. manual에는 gtf 만 받는거 같지만 update가 되서 gff 파일형식도 받는다. 다만 가끔 annotation이 gbk 파일로 되어 있는데 이를 gff 바꾸기 위한 tool이 gff tools이 되겠다.


이 link가 tool

CpG island 를 구할때 이용하는 Obs/Exp CpG 공식

 Obs/Exp CpG = Number of CpG * N / (Number of C * Number of G)

아 별거 아니였다. 여기서 찾음

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 발병의 확률 정도를 알 수 있고 그들의 불확신에서 오는 불안을 해소할 수 있다. 
명쾌하다.