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

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; x_bin*10⁶/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를 생성하는 것이 가장 보정이 잘될 수 있는거라 그런거 같다.

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

Bismark manual 정리

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


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

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

리눅스 커널 내부구조

이것 역시 부서 예산으로 샀다. 워낙에 평들이 좋아서 특히나 쉽게 써져 있다길래 샀는데.. 난 이정도도 이해 못하는 수준인 갑다. 메모리부분이랑 파일시스템은 그럭저럭 어떻게든 보겠는데 태스크 관리는 전혀 뭔소린지 모르겠다.. (아 그래도 인터럽트이니 트랩이니 시그널이니 뭐 이런 구분이 좀 모호 했는데 책을 보고 나면 개념이 선다.) 저자가 엄청 쉽게 쓸려고 노력했음에도 내 수준이 아직 많이 부족한지라.. 이 정도도 이해 못하고.. 좌절한다.. 아.. 뭔책을 봐도 이해 못할거 같은 겁이 나면서.. 뭐 리눅스는 내 평생 안고 가야 하는건데.. 이거 원 언제 한번 속시원함을 느껴 볼려나.. 여튼 내가 본 리눅스 책 중 가장 쉬운 책. 하지만 여전히 이해 안가는 책.

RNA-Seq mapping-first 방식 procedure

cufflinks의 사용법(단순 command 및 약간의 tip) 을 정리한다.
다음은 단지 cufflinks의 manual을 정리한 것이다.
설명 format은 
숫자.명령어 (예: 1.tophat -r 50 ref read1.fastq read2.fastq) 다음 output 및 그 밖의 사항 설명 




1.bowtie-build  ref_seq.fa  ref_seq 
이는 tophat이 mapping을 할때 bowtie를 사용하는데 bowtie의 reference에 대한 index 파일을 생성한다. bowtie manual에 따르면 4개의 파일이 생성되야 하는데 실제론 6개의 파일이 생성된다(정확한 이유는 모르겠음).


2.tophat  -r  50  -o  tophat_result  ref_seq  read1.fq  read2.fq
이를 실행하면 tophat_result 폴더에 tophat의 결과파일이 생성된다.  accepted_hits.bam 파일이 read를 mapping 한 결과 파일이며 이것들을 정리한 것이 deletions.bed, insertion.bed, junction.bed. junction.bed 파일을 보자면 마지막 두 column이 junction의 size와 시작점을 나타낸다. junction의 시작점은 island의 시작점(그러니까 read가 overlap된 bundle)이 아니라 junction을 포함한 read의 시작점이다.


2.1.samtools  index  accepted_hits.bam
read가 align된 것을 samtools의 tview로 보기 위해서는 bam 파일을 indexing 해야 한다.


2.2.samtools  tview  accepted_hit.bam  ref_seq
tview 통해 read의 alignment를 확인할 수 있다.


3.cufflinks  -o  cufflink_result  tophat_result/accepted_hits.bam
결과는 cufflink_result 폴더 안에 transcripts.gtf 라는 파일이 생성. gtf¹ 파일 포멧은 다음 링크에 설명되어 있다.


3.1.cuffmerge  -s  ref_seq.fa  assemblies.txt
위 명령어를 실행하면 merged_asm/merged.gtf 생성. 원래 cufflinks manual에 보면 human의 여러 tissue의 샘플로 위의 3번까지 진행하는데 각 tissue 별로 3번까지 진행해서 나온결과를 merge 한다.


3.2.cuffcompare  -s  ref_seq.fa  -r  known_annotation.gtf  merged_asm/merged.gtf  
이는 이미 알려진 annotation 데이터(여기서는 known_annotation.gtf) 파일과 비교해서 novel한 gene이나 isoform을 찾는다.


위에까지는 그냥 그냥 하는거고 differentially expressed and regulated gene을 찾는 방법 2가지를 workflow를 설명한다. 하나는 novel 한 gene을 안찾고 그냥 기존의 annotation file로다가 찾겠다는 거고, 다른 하나는 novel한 gene까지 찾아가겠다는걸로 후자 역시 기존의 annotation file 이용가능하다. 이건 그냥 cufflinks manual 보자




---------------------------reference-----------------------------
1.gff 파일 : gff 파일은 general feature format의 약자로 여러가지 버젼이 있는데 그중 cufflnks에서는 gtf2(general trascfer format)랑 gff3 만 인식한단다. 


-gtf2 : column은 tab으로 분리 되어야 함