Wednesday, December 21, 2011

head first design pattern 3

  • adapter pattern & facade pattern
    • adapter pattern : 한 클래스의 인터페이스를 클라이언트에서 사용하고자 하는 다른 인터페이스로 변환한다. 어댑터를 이용하면 인터페이스 호환성 문제 때문에 같이 쓸 수 없는 클래스들을 연결해서 쓸 수 있다.
      • 객체 어댑터 : client 클래스에서는 target 인터페이스를 구성 맴버 변수로 가지고 있고 target 인터페이스를 구현한 구상 클래스인 adapter 클래스가 adaptee를 맴버 변수로 가지고 있어서(구성:composition) client에서 target 변수로 adaptee 객체를 쓸수 있도록 target의 기능에 맞게 adaptee 기능들 wrapping 한다.
      • 클래스 어댑터 :  객체 어댑터가 adapter가 adaptee를 구성 맴버로 갖으면서 target 인터페이스를 구현한 것에 반해 클래스 어댑터는 adapter가 target과 adaptee를 다중 상속한다.
    • facade pattern : 어떤 서브시스템의 일련의 인터페이스에 대한 통합된 인터페이스를 제공한다. 퍼사드에서 고수준 인터페이스를 정의하기 때문에 서브시스템을 더 쉽게 사용할 수 있다.
      • 최소지식원칙을 지키는데 퍼사드 패턴을 사용할 수 있다.
  • template method pattern : 템플릿 메소드 페턴에서는 메소드에서 알고리즘의 골격을 정의한다. 알고리즘의 여러 단계 중 일부는 서브클래스에서 구현할 수 있다. 템플릿 메소드를 이용하면 알고리즘의 구조는 그대로 유지하면서 서브클래스에서 특정 단계를 재정의할 수 있다.
    • hook : 추상클래스에서 선언되는 매소드이지만 일반적으로 아무 것도 들어 있지 않다. 서브클래스들은 이 후크를 사용해서 특정 순서에 특정 메소드를 시행할 수 있다. 다르게 말하면 추상클래스에서 hook라는 매소드를 만들고 탬플릿 메소드에서 실행시키면 필요 없는 서브 클래스들은 hook에서 null을 return 하면 되고 특정 메소드를 해야 하는 서브클래스에서는 hook 메소드를 오버라이드해서 사용.
  • iterator pattern & composite pattern
    • iterator pattern : 이터레이터 패턴은 컬렉션 구현 방법을 노출시키지 않으면서도 그 집합체 안에 들어 있는 모든 항목에 접근할 수 있게 해주는 방법을 제공한다.
      • 예제 : 두 음식점이 합병됐다. 메뉴판을 합쳐야 하는데 메뉴 클래스는 두 음식점이 통일했다. 하지만 두 음식점이 매뉴들를 담고 있던 컬랙션이 다르기 때문에 웨이트리스라는 객체를 만들었을때 메뉴를 출력하려면 각 컬랙션을 가각 for문을 돌려야 하는 현상이 나타난다.
      • 이터레이터 패턴 : iterator 라는 인터페이스를 만들고 각 음식점의 컬랙션에 맞는 구상 iterator 클래스를 만든다. 각 음식점 클래스에서는 이 자신의 메뉴 컬랙션에 맞는 구상 iterator 클래스를 return 한다. 그러면 웨이트리스 객체는 동일한 메소드를 가지고 두 음식의 메뉴 리스트를 출력 가능하다.
    • composite pattern : 컴포지트 패턴을 이용하면 객체들을 트리 구조로 구성하여 부분과 전체를 나타내는 계층구조로 만들 수 있다. 이 패턴을 이용하면 클라이언트에서 개별 객체와 다른 객체들로 구성된 복합 객체를 똑같은 방법으로 다룰 수 있다.



디자인 원칙 : 
  1. 최소 지식 원칙 : 정말 친한 친구하고만 얘기해라. 아무 객체의 매소드를 호출하는 것이 아니라, 특정 객체(1.객체 자신, 2.메소드에 매개변수로 전달된 객체, 3.그 메소드에서 생성하거나 인스턴스를 만든 객체, 4.그 객체에 속하는 구성요소) 의 메소드만 호출한다.
  2. 헐리우드 원칙 : 먼저 연락하지 마라. 우리가 연락하겠다. 이는 고수준 구성요소에서 저수준 구성요소에 하는 소리. 
  3. 클래스를 바꾸는 이유는 한 가지 뿐이어야 한다.

Saturday, December 17, 2011

head first design pattern 2

  • Factory pattern
    • Factory method pattern : factory method pattern에서는 객체를 생성하기 위한 인터페이스를 정의하는데, 어떤 클래스의 인스턴스를 만들지는 서브클래스에서 결정하게 만든다. 팩토리 메소드 패턴을 이용하면 클래스의 인스턴스를 만드는 일을 서브클래스에게 맡기는 것이다.
    • abstract factory pattern : 추상 팩토리 패턴에서는 인터페이스를 이용하여 서로 연관된, 또는 의존하는 객체를 구상 클래스를 지정하지 않고도 생성할 수 있다.
    • factory method pattern은 상송을 통해서 객체를 만들고 abstract factory pattern은 객체 구성(interface의 구현)을 통해서 만든다. 
    • 예제 : PizzaStore을 예로 든다. factory method pattern 은 PizzaStore를 추상클래스로 만들어서 이 클래스의 추상 메소드를 derived class 들이 구현하게 한다. abstract factory pattern에서는 PizzaIngrediantFactory 라는 인터페이스를 만들어서 이를 정의 하는 구상 클래스를 만든다.  이때 PizzaStore의 derived class들은 각기 다른 PizzaIngrediantFactory를 맴버변수로 갖고 있어서 이 PizzaIngrediantFactory를 Pizza 객체를 만들때 인자로 넘겨준다. factory method pattern의 PizzaStore에서는 지역별 Pizza 클래스가 있어서 string(예를 들어 "cheese", "clam") 등을 인자로 넘기는데 반해 abstract factory pattern의 경우 지역별 Pizza 클래스가 있는것이 아니라 종류별 Pizza 클래스(예를 들어 "CheesePizza") 가 있어서 PizzaIngrediantFactory 객체를 인자로 넘긴다.
  • singleton pattern : 싱클턴 패턴은 해당 클래스의 인스턴스가 하나만 만들어지고, 어디서든지 그 인스턴스에 접근할 수 있도록 하기 위한 패턴이다.
    • 필요성 : 스레드 풀, 캐시, 대화상자, 사용자 설정, 레지스트리 설정등을 처리하는 객체는 여러개의 객체가 필요하지 않고 오히려 하나의 객체만 있어야 한다.
    • 단순 답변 : 생성자를 private로 하고 생성자를 호출하는 getInstance 라는 static 함수를 만들고 이 함수가 호출 되었을때 private로 되어 있는 생성자를 호출해서 자신의 객체를 만들고 이를 static 맴버 변수에다 넣어서 return 한다.
    • 문제점 : 멀티스레드에서 접근할 때 하나 이상의 객체가 생성될 가능성이 여전히 존재
    • singleton pattern : 1. getInstance 메소드 자체를 synchromized 키워드를 이용해서 동기화 한다. 혹은 2.static 맴버 변수에다가 프로그램이 시작 하자마자 클래스 객체를 생성해서 집어 넣는다. 혹은 3. DCL(double-checking locking)을 사용한다. 
  • command pattern : 커맨드 패턴을 이용하면 요구 사항을 개체로 캡슐화 할 수 있으며, 매개변수를 써서 여러 가지 다른 요구 사항을 집어넣을 수도 있다. 또한 요청 내역을 큐에 저장하거나 로그로 기록할 수도 있으며, 작업취소 기능도 지원 가능하다.
    • 예제 : 특정 리모컨이 있는데 여기에는 10개의 슬롯이 있고 각 슬롯에는 on/off 스위치가 두개씩 있다. 그리고 undo 라는 스위치가 하나 있다. 이 리모컨을 만능 리모콘이 되도록 코딩한다. 그러니까 어떤 슬롯은 형광등의 스위치 버튼이 되고 어떤 슬롯은 자동문의 버튼이 되는 것. 문제는 어떤 슬롯이 어떤 일을 할지 모른다는것. 결국 필요한건 메소드의 캡슐화
    • command pattern : client는 receiver를 command 객체에 넣어서 이 command 객체를 invoker에게 넘기면 invoker는 command 객체를 실행 시킨다. command 는 interface 이고 execute라는 virtual 함수가 있고 이를 구현한 구상 객체들이 존재한다. 실질적으로 이 구상 객체들이 execute 라는 함수를 구현한다. 이 execute 함수 안에 각 command가 구성 객체로 받는 receiver들의 함수들을 호출하는 것이다. 결국 invoker는 command를 받아서 command 객체의 execute 만 실행 할 뿐 어떤 메소들이 실행되는지는 알 수가 없다. 위 예제를 비유하자면 client는 리모컨을 사용하는 사람이 되겠고 invoker가 리모컨, receiver가 리모컨이 컨트롤 하는 각종 기기들의 function을 정의한 클래스가 되겠고, command는 단지 리모컨의 메소드 캡슐화를 위한 실제적인 receiver 핸들 객체가 되겠다.

Thursday, December 15, 2011

head first design pattern 1

head first design pattern에 정의된 내용.
  • strategy pattern : strategy pattern 은 알고리즘군을 정의(인터페이스)하고 각각을 캡슐화하여 교환해서 사용할 수 있도록 만든다. 이를 이용하면 알고리즘을 사용하는 클라이언트와는 독립적으로 알고리즘을 변경할 수 있다.
    • 예제 : 연못의 오리 게임. Duck 이라는 super class를 만들고 이를 상속 받은 여러 종류의 오리 class가 있다. 원래 Duck class에는 없는 기능을 넣을 려면 어떻게 해야 하는가?
    • 단순 답변 : Duck class에 특정 기능의 함수를 넣어준다. 하지만 이는 Duck class를 상속 받은 클래스 모두에게 이 기능이 생기게 된다. 만약 특정 derived class는 이 기능이 변형되어야 한다면?
    • strategy pattern : 위 문제는 함수 오버라이딩으로 가능하기는 하지만.. 이보다는 변형 가능한 것들을 하나의 인터페이스로 만든 다음 이것을 Duck class 의 맴버 변수로 넣는다. 구현은 이 인터페이스를 상속받는 여러 클래스들에서 구현. 
  • observer pattern : observer pattern 에서는 한 객체(subject)의 상태가 바뀌면 그 객체에 의존하는 다른 객체들한테 연락이 가고 자동으로 내용이 갱신되는 방식으로 일대다 의존성을 정의한다.
    • 예제 : 기상스테이션 예제로 WeatherStation이라는 클래스가 있고 이 클래스에서는 기상변화를 관찰 한다. 기상 정보는 실시간으로 WeatherData라는 클래스의 객체로 넘어가고 이때 WeatherData 객체는 이 업데이트된 정보를 각각의 Display 객체들에게 통보해야 한다. 
    • observer pattern : Subject 라는 인터페이스를 만들고 그 안에는 Observer 라는 interface를 갖고 있다. WeatherData 클래스는 Subject 인터페이스를 구현하고 있다. Observer 는 여러 Display 클래스들에 의해 구현되어 지고 이 Display 클래스들은 WeatherData 클래스 객체를 맴버 변수고 갖고 있어서 생성자가 호출될 때 WeatherData 객체를 넘겨 받아서 넘겨받은 WeatherData의 add 함수를 호출해서 WeatherData의 observer 맴버 변수 배열에 Display 자기 자신을 등록한다.
  • Decorator pattern : decorator pattern에서는 객체에 추가적인 요건을 동적으로 첨가한다. 데코레이터는 서브클래스를 만드는 것을 통해서 긴능을 유연하게 확장할 수 있는 방법을 제공한다.
    • 예제 : 스타버즈라는 커피숍의 주문 시스템을 만든다. 각 매뉴들을 표현하고 가격을 구현해야 한다. 
    • 단순 답변 : Beverage라는 베이스 클래스를 만들고 새로 생기는 메뉴마다 이 슈퍼 클래스를 상속받는 서브 클래스를 만든다. 아니면 베이스 클래스에다가 서브 클래스에서 쓸 맴버 변수들을 포함시킨다. 이러한 경우 서브 클래스가 엄청 많아진다거나 혹은 서브 클래스가 필요하지도 않은 많은 맴버 변수들을 상속받게 된다.
    • decorator pattern : 문제를 파악해보면 DarkRoast, HouseBlend, Decaf는 기본이 되는 음료가 되고 여기에다가 Mocha나 Whip (휘핑크림) 등은 decorator가 되겠다. 모든 음료의 기본이 되는 베이스 가상 클래스인 Beverage를 정의 하고 각종음료(DarkRoast, HouseBlend, Decaf)는 이를 상속하고 구현하는 구상 클래스가 되겠다. 그리고 decorator들(Mocah, Whip) 역시 Beverage 가상 클래스를 상속하게 한다. 이것의 목적은 Beverage로 부터 행동을 상속 받기 위함이 아니라 decorator의 형식과 각종 음료의 형식을 맞취기 위한 것이다. 확장되는 행동들은 decorator 안에 구현하고 decorator 안에는 또한 Beverage를 구성 요소로 갖는다. 결론적으로 기본음료를 decorator가 감싸고 이 감싼 객체를 또다른 decorator가 감싸는 형식으로 구현하는 것이다. 


디자인 원칙
  1. 애플리케이션에서 달라지는 부분을 찾아내고, 달라지지 않는 부분으로부터 분리 시킨다.
  2. 구현이 아닌 인터페이스에 맞춰서 프로그래밍한다.
  3. 상속보다는 구성을 활용한다. 
  4. 클래스는 확장에 대해서는 열려 있어야 하지만 코드 변경에 대해서는 닫혀 있어야 한다.

    Tuesday, October 11, 2011

    GWAS

    GWAS 란?
    다음을 요약해보면..
    일단 GWAS는 genome wide association study의 약자. 뭔말인고 하니 genome wide 하게 association study를 한다는 의미. 이건 또 무슨 말인고 하니 genome wide 하게 DNA를 관측하고 이 DNA의 차이 혹은 변화를 질병과 같은 변수와 연관지여 연구를 한다는 의미.
    사람의 genome은 서로 99.9%가 동일하다. 그래서 차이는 거의 없을 거고 이 차이가 나름 중요한 역할을 할거다(질병이라던지 관련해서). 그럼 이 개인간의 차이 즉 common SNP은 몇개나 있을까? 10 million쯤 있을거라 생각되고.. 그럼 GWAS를 연구할려면 이 10M SNP를 다 detecting 해야 하냐? 아니다. 과학자들이 HapMap project라는 것을 진행해서 SNP의 조합인 haplotype들을 찾아놨기 때문에 돈은 많이 싸졌다. 물론 요즘 NGS 값이 워낙 싸서 더 싸졌다.




    그럼 어떤식으로 GWAS를 진행할까?
    많은 다수 그러니까 수천명 규모에서 수만명까지 특정 질병군(case)와 일반군(control)의 사람이 필요하고(이때 confounding factor인 성, 종족등을 통일(?)시킨다)


    1.각각의 사람을 genotyping 한다. 그러니까 SNP를 찾는다.
    2.quality control of genotype data.
    3.genotype과 phenotype간의 통계적으로 연관이 있는지 확인. chi-squared 같은 걸 이용
    4.연관된 주변의 추가적인 SNP의 genotyp으로 fine mapping of association signal. 그니까 3번 step에서 유의하다고 생각되는 snp 주변의 SNP들을 추가적으로 genotyping 해서 haplotype을 경험적으로다가 만들어 낸다.
    5.그 담에 다른 집단에다가 적용. 이땐 지정된 SNP로다가만 testing 한다.
    6.biological validation




     GWAS tools
    많다. CRAN 의 task view의 genetics에서도 보면 GWAS를 할 수 있는 많은 package를 소개한다. 또한 broad institute 에서 만든 PLINK도 많이 사용한다.

    Friday, September 23, 2011

    R multiple plots in same window

    맨날 까먹어.. 난.. 인터넷이 폐해인가? 아님 내가 원래 붕어인건가?




    check point 1
    "par 함수의 사용"


    일단은 한 윈도우 안에 여러 그림을 그릴려면 par 라는 함수를 사용해야 하고 이건 뭐하는 거냐 하면 query for graphical parameter. 각종 device 마다의 parameter가 있는데 이걸 par 함수의 attribute라고 해야 하나 옵션이라고 해야 하나.. 이런 것들을 이용해서 query로 던지는 것. 특히 mfrow(or mfcol)과 fig 옵션을 사용하면 각 plot의 위치와 범위를 선정 할 수 있다. 다음의 사이트와 블로그에 잘 설명이 되어 있는데 간략하게 설명하자면 


    par(mfrow=c(2,2)) #multiple figures by row로 multiple plot이 몇 x by x 인지 또 각 plot들어가는 순서를 나타낸다.
    par(fig=c(0,0.3,0.5,1)) #앞의 두 vector는 x 축의 coordinate, 뒤 두 vector는 y 축의 coordinate를 의미하는 것으로 0과 1 사이 값을 갖게 된다.아래 그림 참조
    plot(somthing) # par의 fig로 영역을 잡은 곳에 그림을 그림


    %아래 그림은 "불탄五징어" 님 블로그에서 발췌함을 밝힌다.%
    그런데 여기서 fig의 coordinate의 기준이 어디냐? 아래 그림을 보면 NDC space가 fig가 따르는 체계(?)이다. 다음 참조




    check point 2
    "heat map 그리기"


    사실 이건 multiple plot이랑 상관없지만 또 잊어버릴까봐 기록한다.




    check point 3
    "histogram horizontal"


    hist 자체에는 horiz 이라는 parameter가 없는데 이를 해결하는 것이 hist의 return 값을  barplot을 이용하여 그리는 것. 근데 문제가 y축의 label을 쓰는 것이 생각보다 골치 아픈데.. 다음을 참조하면 해결 가능해진다.


    h.temp = hist(something$value, plot=FALSE) #plot 옵션을 F로 줘서 그림은 그리지 않은채 hist의 return 값만 받는다.
    barplot(h.temp$count, space=0,horiz=TRUE) #space는 hist가 원래 space가 0이므로 가능한 비슷하게 그리기 위해 준 옵션이고 horiz를 통해서 rotated hist를 그린다.
    labels = h.temp$breaks #y축에 label을 할 숫자들
    text(par('usr')[1]-3,1:length(labels)-1,srt=90,adj=0.5,labels=labels,xpd=T) #첫 두 숫자는 x,y의 coordinate로 label의 갯수만큼 vector를 생성해야 한다(여기의 경우 y vector를 label수 만큼 지정했기 때문에 x 축의 값은 하나라도 반복되서 적용). srt는 label의 글씨를 90도 기울인 것이고 adj는 label의 가운데와 tick과 정렬시키는 것. 


    Conclusion
    결국 내가 하려던 것은 The DNA methylome of Human Peripheral Blood Mononuclear Cells에 나오는 figure를 그리고 싶었던 것인데.. 최종 코드는 아래와 같다.

    Thursday, September 22, 2011

    UCSC GENE table cdsStartStat, csEndStat SQP type 설명

    아 cdsStartStat과 cdsEndStat의 필드의 데이터가 enum('none', 'unk', 'incmpl', 'cmpl')로 되어 있는지 첨 알았다. 항상 exonStarts와 exonEnds와 같은 형식으로 되어있는줄 알았는데 parsing 하다보니까 알게됐네.. 


    이 enum('none', 'unk', 'incmpl', 'cmpl') 는 뭐야 도대체..


    다음 링크를 따라가 본다.

    • none - no CDS specified from the sequence's data source.
    • unk - unknown - not known if CDS start/end is complete.
    • incmpl - the CDS start/end is incomplete
    • cmpl - the CDS start/end is complete.

    Monday, August 29, 2011

    GFF tools

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


    link가 tool

    Thursday, August 18, 2011

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


     Obs/Exp CpG = Number of CpG * N / (Number of C * Number of G)
    
    
    아 별거 아니였다. 여기서 찾음

    Wednesday, August 17, 2011

    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)



    Friday, August 12, 2011

    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

    Thursday, August 11, 2011

    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 의 메모리를 할당한다. 


    Wednesday, August 10, 2011

    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 결과 정리

    Monday, August 8, 2011

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

    Wednesday, July 27, 2011

    STR (short tandem repeat)

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


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


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


    STRBaseDB논문은 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을 보인단다. 뭐 여튼 난 알고 싶은건 길이니까. 아래 표를 참고하자.



    Monday, July 25, 2011

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

    Thursday, July 21, 2011

    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 인가를 증명한다. 이건 생략

    Monday, July 18, 2011

    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 함을 알 수 있다.

    Saturday, July 9, 2011

    MEDIPS manual 정리 + 대략의 MBD-Seq pipeline

    <MEDIPS setting>
    이게 R package이기 때문에 기본적으로 R을 깔아야 하고 그 밖에 package를 깔아야 하는데 manual 따라하면 된다




    <MEDIPS input file format>
    헤더 없이 tab-separated format.


    chr22   16052499        16052535        +
    chr22   16053992        16054028        -
    chr22   16054824        16054860        +
    chr22   16056675        16056711        -



    <command procedure>

    0.
    fastq quality format 체크. 일반적으로는 fastq quality file format를 알고 시작하지만.. 내가 갖고 있는 보건원 자료는 전혀 그런 단서가 없기에.. 일단 체크 하고 한다. 이는 mapping 프로그램의 mapping quality를 구할때 인자로 이용된다. quality file format에 대해서는 이 논문에서 자세하게 나온다.

    1.
    bowtie 같은 프로그램으로 mapping(여기서는 bowtie 사용).

    2.
    데이터 필터, quality와 duplication read 제거


    3.
    MEDIPS 프로그램 input 파일 만들기. 이게 일종의 BED 파일 형식인데.. 방법은 여러가지가 있을 수 있다. 여기를 참조하면 여러가지 방법이 있을 수 있다. 난 pysam 이용해서 script 만들었다(-u : unambiguously mapped read 만 선택).


    4.여기서 부터가 이제 MEDIPS package 사용에 대한 것이다. 다음은 UCSC 의 hg19에다가 read를 mapping 했다고 가정했을때의 분석 command이다.


    -특정 genome의 sequence를 R로 loading 하는 것은 BSgenome 패키지 이용.
    >library(BSgenome) #BSgenome package loading
    >abailable.genomes() #이용가능한 genome list 확인
    >library(BSgenome.Hsapiens.UCSC.hg19) #만약 이 데이터가 없다면 source("http://bioconductor.org/biocLite.R"); biocLite("데이터")로 데이터 다운로드


    -MEDIPS workflow
    >library(MEDIPS)
    >library(BSgenome.Hspiens.UCSC.hg19)
    > file = system.file("extdata", "MeDIP_hESCs_chr22.txt", package = "MEDIPS") #system.file 은 R의 기본 function 인데 특정 package의 subdirectory의 파일의 경로를 알려줌(여기서 "extdata"가 subdirectory)
    > CONTROL.SET = MEDIPS.readAlignedSequences(BSgenome = "BSgenome.Hsapiens.UCSC.hg19", file = file, numrows = 672866) #여기까지만 하면 CONTROL.SET에는 input file의 정보만 들어있는것
    > CONTROL.SET = MEDIPS.genomeVector(data = CONTROL.SET, bin_size = 50, extend = 400) #bin_size로 binning 한 genome의 각 bin에 extend로 확장한 read들 중 overlapped read 갯수를 할당해서 raw MeDIP-Seq signal로 인식
    > MEDIPS.exportWIG(file = "output.rpm.control.WIG", data = CONTROL.SET, raw = T, descr = "hESCs.rpm") #raw=T 옵션에 의해 raw MeDIP-Seq signal의 rpm(reads per million; xbin*106/n) value를 wig file로 export
    > sr.control = MEDIPS.saturationAnalysis(data = CONTROL.SET, bin_size = 50, extend = 400, no_iterations = 10, no_random_iterations = 1) #read 양이 충분한지 확인하기 위한 분석. 총 read를 random 하게 두 집단으로 나누고 각 집단에서 read를 random하게 뽑아서 두 집단간의 correlation을 측정한다( no_iterations 숫자 만큼 반복, 이때 뽑는 read양을 늘려감), 이때 saturation analysis는 총 read를 두 그룹으로 나눠서 했기 때문에 총 read에 대해서 했다기 보다는 half read에 대해서 한것이므로 total read의 saturation analysis를 위해 두 그룹의 read를 doubling 해서 saturation analysis 를 함, 이것이 estimated saturation analysis. 결과를 이해하기 위해 manual 꼭 참조 바람, 다른 인자들도 고려 해야 함(no_random_interations, empty_bins, rank)
    > MEDIPS.plotSaturation(sr.control) #ploting
    > CONTROL.SET = MEDIPS.getPositions(data = CONTROL.SET, pattern = "CG") #MeDip-Seq 의 signal이 CpG density의 영향을 받기 때문에 CpG 에 대한 정보를 넣어주는 함수. plus stand의 CpG start site를 다 찾는다. CpG 이외의 pattern도 찾을 수 있다. 이때는 plus, minus strand에 대해서 다 찾음(CpG 이외에는 서로 complementary 하지 않으므로)
    > cr.control = MEDIPS.coverageAnalysis(data = CONTROL.SET, extend = 400, no_iterations = 10) #getPositions 함수에 의해 추가된 pattern(여기서는 CpG)이 region(read를 extend 한것)에 의해 얼마나 커버가 되는지 확인하기 위한 함수. no_interations 수만큼 반복(이 때 역시 randomly selected region수를 늘려간다), 이 역시 결과를 해석하기 위해선 manual 참조
    > MEDIPS.plotCoverage(cr.control) #ploting result
    > er.control = MEDIPS.CpGenrich(data = CONTROL.SET) #reference에 비해 MeDIP-Seq에의해 CpG가 얼마나 enrichment되어 있는가를 검사. 전체 genome과 region이 커버하는 영역에 대해 C,G,CpG 갯수를 세서 relative CpG의 값과 observed/expected value를 구하고 또 genome과  region의 커버 영역의 relative CpG의 값과 observed/expected value의 비율을 구함, 경험적으로 봤을때 $enrichment.score.relH의 값이 2.547679가 나옴
    > CONTROL.SET = MEDIPS.couplingVector(data = CONTROL.SET, fragmentLength = 700, func = "count") #CpG density에 따라 MeDIP-Seq의 mC에 대한 affinity가 달라지므로 이를 보정하기 위한 함수. 각 pre-defined genomic bin에 대해 [bin_position-fragmentLength, bin_position+fragmentLength] 범위 안의CpG (getPositions 함수에서 등록한 pattern)의 density를 구한다(fragmentLength는 insertsize). 계산 방법은 단순하게 range 안의 CpG 갯수의 총합을 구하는 방법("count")과 bin과의 거리에 따라 weight를 주는 방법이 있다("linear","exp","log" 혹은 "custom", manual 확인, 여기서 fragmentLength랑 func 테스트 한것 확인
    > MEDIPS.exportWIG(file = "PatternDensity.WIG", data = CONTROL.SET, pattern.density = TRUE, descr = "Pattern.density") #coupling vector를 WIG format으로 export
    > CONTROL.SET = MEDIPS.calibrationCurve(data = CONTROL.SET) #MeDIP-Seq signal intensity와 CpG density의 dependency를 측정. coupling factor들을 regular level로 나누고 각 level안에 속하는 genomic bin의 평균 signal intensity를 x 축으로 평균 coupling factor를 y축으로 해서 calibration curve를 계산한다.
    > MEDIPS.plotCalibrationPlot(data = CONTROL.SET, linearFit = T, plot_chr = "chr22") #plotting calculation curve, 시간 엄청걸림. 그러니까 png() 함수 이용. 아.. 여기서 결정적인 말이 나온다. plot을 하면 low CpG density에서만 linear한 correlation이 있고 higher CpG에서는 마치 그러한 관계가 없는것처럼 나오는데 이는 실제로 higher CpG 에서는 correlation이 없는게 아니라 실제로 methylation이 많이 안되서 마치 관계가 없는것처럼 나온다는 것(이는 bayesian deconvolution 논문에 나온단다) 그래서 higher CpG density에서도 linear dependency가 있는 것으로 가정. linear regression 구하는건 manual 참조.
    > CONTROL.SET = MEDIPS.normalize(data = CONTROL.SET) #calculationCurve 함수에서 linear regression의 a,b(y 절편과 기울기) 로 구해진 estimated x로 normalized signal를 구한다, 그러니까 원래 signal에서 estimated x의 값을 나누어 줌. 아래 식 참조. 여기까지 하면 MEDIPS SET의 모든 slot이 채워지게 된다.

    > MEDIPS.exportWIG(file = "output.rms.control.WIG", data = CONTROL.SET, raw = F, descr = "hESCs.rms") #wig 파일로 normalized relative methylation score를 export
    > file = system.file("extdata", "hg19.chr22.txt", package = "MEDIPS") #ROI(region of interest) 파일. 아래와 같은 형식으로 만들면 된다. 마지막 컬럼은 identifier (ID).

    chr22   25170187        25171687        ENST00000332271
    chr22   18874011        18875511        ENST00000455288
    chr22   18874337        18875837        ENST00000424147
    chr22   19928833        19930333        ENST00000400521
    > promoter = MEDIPS.methylProfiling(data1 = CONTROL.SET, ROI_file = file, math = mean, select = 2) #methylProfiling 함수가 ROI 분석 뿐 아니라 전체 게놈에 걸쳐서 frame_size를 window 사이즈로 step 인자만큼 슬라이딩하면서 methylation profile을 생성한다. select가 1이면 rpm을 기준으로 2이면 rms를 기준으로  math (여기서는 mean) 방식으로 ROI의 MeDIP-Seq signal의 평균값을 계산한다. 다른 인자들 manual로 확인, 특히 tranf (이는 rpm내지는 rms 값을 log2를 취한후 [0,1000] 사이의 값이 되도록 변형. 곧 ROI 중 최대 rms 값을 같는 것을 1000으로 최소값을 0으로 해서 normalize 한다). 그리고 이 함수의 결과 ams(absolute methylation signal) 값도 report가 되는데 그는 다음과 같다. 
    > hist(promoter$coupling, breaks = 100, main = "Promoter CpG densities", xlab = "CpG coupling factors") #ROI의 CpG density에 대해 histogram을 그린것 
    > hist(promoter$rpm_A[promoter$rpm_A != 0], breaks = 100, main = "RPM signals", xlab = "reads/bin") #rpm을 histogram으로 표현한것 
    > hist(promoter$ams_A[promoter$ams_A != 0], breaks = 100, main = "AMS signals", xlab = "absolute methylation score (ams)") #ams을 histogram으로 표현한것 
    > frames.frame500.step250 = MEDIPS.methylProfiling(data1 = CONTROL.SET, frame_size = 500, step = 250, math = mean, select = 2) #ROI 말고 전체 genome에 대해 frame_size를 window로 해서 step 만큼 bp를 슬라이딩해서 methylation profiling을 구할 수있는데 이것이 시간이 많이 걸리니까 결과를 아래 명령어 write.table로 저장
    > write.table(frames.frame500.step250, file = "frames.chr22.meth.txt", sep = "\t", quote = F, col.names = T, row.names = F)
    > frames.frame500.step250 = read.table(file = "frames.chr22.meth.txt", header = T) #이건 이제 저장한 파일을 다시 R의 변수로 loading 하는것


    여기서부터는 하나의 MEDIPS set이 아니라 두개 이상의 MEDIPS set 일때 분석 command line으로 자세한 설명은 생략하고 command만 나열하겠다. 여기서는 단순히 새로운 SET을 추가하는데 중점을 맞추기 때문에 새로운 SET에 대해 quality control이나 그밖에 analysis를 하지 않는다. 하지만 두번째 sample과 input에 대해서도 첫번째와 마찬가지로 모든 analysis를 하길 권장한다.


    > file = system.file("extdata", "MeDIP_DE_chr22.txt", package = "MEDIPS")
    > TREAT.SET = MEDIPS.readAlignedSequences(BSgenome = "BSgenome.Hsapiens.UCSC.hg19", file = file, numrows = 863054)
    > TREAT.SET = MEDIPS.genomeVector(data = TREAT.SET, bin_size = 50, extend = 400)
    > TREAT.SET = MEDIPS.getPositions(data = TREAT.SET, pattern = "CG")
    > TREAT.SET = MEDIPS.couplingVector(data = TREAT.SET, fragmentLength = 700, func = "count")
    > TREAT.SET = MEDIPS.calibrationCurve(data = TREAT.SET)
    > TREAT.SET = MEDIPS.normalize(data = TREAT.SET) #여기까지 두번째 SAMPLE 처리
    > file = system.file("extdata", "Input_StemCells_chr22.txt", package = "MEDIPS")
    > INPUT.SET = MEDIPS.readAlignedSequences(BSgenome = "BSgenome.Hsapiens.UCSC.hg19", file = file, numrows = 352756)
    > INPUT.SET = MEDIPS.genomeVector(data = INPUT.SET, bin_size = 50, extend = 400) #여기까지만  input set에 대한 process를 해야한다. input set은 antibody 처리를 하지 않은 그냥 sonication 까지만 한 DNA fragment이기 때문에 CpG bias를 교정 할 필요가 없다. 그래서 rpm 값을 그냥 사용하면 된다
    > MEDIPS.exportWIG(file = "output.rpm.input.WIG", data = INPUT.SET, raw = T, descr = "INPUT.rpm")
    > sr.input = MEDIPS.saturationAnalysis(data = INPUT.SET, bin_size = 50, extend = 400, no_iterations = 10, no_random_iterations = 1)
    > MEDIPS.plotSaturation(sr.input) #saturation plot을 그리면 MeDIP 데이터 보다 input 데이터가 reproducibility가 천천히 증가함을 볼 수 있다(이는 genomic complexity가 높아졌으므로; 곧 methylation 부분만 시퀀싱 된게 아니라 전체 genome을 시퀀싱 한 결과가 되므로).
    > INPUT.SET = MEDIPS.getPositions(data = INPUT.SET, pattern = "CG") 
    > cr.input = MEDIPS.coverageAnalysis(data = INPUT.SET, extend = 400, no_iterations = 10)
    > MEDIPS.plotCoverage(cr.input)
    > er.input = MEDIPS.CpGenrich(data = INPUT.SET)
    > INPUT.SET = MEDIPS.couplingVector(data = INPUT.SET, fragmentLength = 700, func = "count")
    > INPUT.SET = MEDIPS.calibrationCurve(data = INPUT.SET)
    > png("CalibrationPlotINPUT.png")
    > MEDIPS.plotCalibrationPlot(data = INPUT.SET, linearFit = F, plot_chr = "chr22")
    > dev.off() #여기까지는 그냥 input과 MeDIP 데이터 간의 차이를 보기 위한 plot을 그리기 위한 것이다


    > diff.meth = MEDIPS.methylProfiling(data1 = CONTROL.SET, data2 = TREAT.SET, input = INPUT.SET, frame_size = 500, select = 2) #3개의 data set에 대해 전체 genome에 대해 methylation profile을 구한 것, 이거 시간 엄청 걸린다. 자세한 내용은 manual 참조(무엇 무엇이 계산되는지, 어떤 확률을 쓰는지 등). input 이 있어야 나중에 filtering 할때 background 때의 threshold를 잡게 된다.
    > write.table(diff.meth, file = "diff.meth.txt", sep = "\t", quote = F, col.names = T, row.names = F)
    > diff.meth = read.table(file = "diff.meth.txt", header = T)


    여기서 부터 DMR을 찾고 annotation하는 것인데 manual 참조
    > diff.meth.control = MEDIPS.selectSignificants(frames = diff.meth, control = T, input = T, up = 2, p.value = 1e-04, quant = 0.99) #이 함수에 인자가 많음, 반드시 확인 해야 한다
    > diff.meth.control.merged = MEDIPS.mergeFrames(frames = diff.meth.control) #이 예에서는 frame이 overlap되어 있기 때문에 이를 합치는 작업을 한다. 이는 단순히 coordinate를 나타내는 것 뿐이지 다른 수치들은 전혀 들어가지 않는다.
    > write.table(diff.meth.control, file = "DiffMethyl.Up.hESCs.bed", sep = "\t", quote = F, row.names = F, col.names = F)
    > file = system.file("extdata", "hg19.chr22.txt", package = "MEDIPS")
    > diff.meth.control.annotated = MEDIPS.annotate(diff.meth.control, anno = file)
    > length(unique(diff.meth.control.annotated[, 4]))
    > diff.meth.treat = MEDIPS.selectSignificants(frames = diff.meth, control = F, input = T, up = 2, down = 0.5, p.value = 1e-04, quant = 0.99)
    > write.table(diff.meth.treat, file = "DiffMethyl.Up.DE.bed", sep = "\t", quote = F, row.names = F, col.names = F)
    > file = system.file("extdata", "hg19.chr22.txt", package = "MEDIPS")
    > diff.meth.treat.annotated = MEDIPS.annotate(diff.meth.treat, anno = file)
    > length(unique(diff.meth.treat.annotated[, 4]))




    <caution>
    1.MEDIPS는 bed 파일의 start position이 1부터 시작. 일반적인 bed 파일과는 다르다.
    2.input 파일은 무조건 하나. 여러개 파일이면 합쳐야 한다.
    3.MEDIPS 각 함수의 인자들 전부 체크해야 한다.


    <check list>
    1. flat_400_700.tab 파일이 어떻게 나온건지
    2.couplingVector 함수에서 fragnment=700, func="count"일때 HEP의 bisulphite 데이터와 the best negative correlation이었다는데.. 이게 무슨 의미인지. ==> 아.. 이거 생각해보면 couplingVector가 CpG density에 의해 bias가 생기는 걸 보정하는 것이므로 bisulphite 데이터와 negative correlation의 값이 제일 크게 나오게끔 couplingVector를 생성하는 것이 가장 보정이 잘될 수 있는거라 그런거 같다.

    Friday, July 8, 2011

    MBD-Seq, BS-Seq 고려사항

    0.MBD-Seq에서 insert size를 어떻게 넣을 것인지. MEDIPS를 이용한다고 가정햇을때 single end 라면 insert size를 어떻게 판단할 것인지(spike in sample을 넣는등의 방법 고려), paired-end 라면 left read에다가 right read의 end position을 넣어서 마치 single end 데이터처럼 만들어서 extend 인자를 사용하지 않는식의 방법 고려.

    1.MBD-Seq 할때 MEDIPS로 calibration curve를 그려서 CpG density dependency 체크해서 normalization 할지 결정.

    2.MEDIPS 이용시 ROI 들을 정해서 biomart나 아니면 ucsc 에서 받아서 ROI 별 methylation level를 측정한다.

    2.genome 상의 repeat region 을 어떻게 처리 할지 고려 : YH 논문에서는 local copy number 개념으로 접근해서 repeat region을 제거했고, 또 다른 논문에서는 아예 RepeatMasker로 genome 상의 repeat 부분을 mask 했다.

    3.read duplication 제거

    4.sampling curves for confirm read sufficiency

    Thursday, July 7, 2011

    Bismark manual 정리

    음.. 이 포스팅 역시 보건원 과제에 의한 것으로.. BS-Seq 데이터를 reference에 mapping 할때 쓰는 프로그램인 Bismark의 manual을 정리한 것이다. 일반적으로 BS-Seq 한 논문들 보면 딱히 published 된 프로그램 말고 그냥 reference C->T로 바꾸고 read도 바꾸고 해서 그냥 각종 mapping 프로그램 사용하는데.. 어짜피 내가 짜는 것보단 published 된게 더 나을거 같다 싶어서 가장 최신으로 보이는 Bismark를 사용하기로 한다. 근데 뭐 잘 안맞으면 새로 만들어야 하니까..


    사이트는 여기. 음 manual 보니까 Quick reference가 있다.

    Monday, July 4, 2011

    SRA 데이터 변환 toolkit

    안타깝게도.. 
    오늘 GMI에서는 네이쳐에 십수명의 사람들의 RNA-Seq 분석해서 보니까 보전적인 RNA-modification 있다는 식으로 논문이 났다는데(아직 논문의 abstract 조차 보지 않았다. 금선생한테 들은거.. 내 조만간 얼마전 science에 나온 RNA-editing 조사한 논문이랑 이 논문이랑 비교해서 포스팅 하리라).. 
    본사 마크로젠에서는 RNA-Seq 데이터 조차 찾을 수 없어서.. cufflinks를 이용해 보려 했다만.. 데이터가 없는 관계로 논문의 SRA 데이터 다운 받아서 쓰려는데.. 


    뭐 여튼 한없이 안타까운 맘은 뒤로 하고 SRA 데이터 변환 프로그램이 NCBI에서 제공한다. 혹시나 잊어버리지 않을까 하는 맘에 링크를 남긴다.


    http://www.ncbi.nlm.nih.gov/books/NBK47540/#SRA_Download_Guid_B.5_Converting_SRA_for