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

Sunday, July 3, 2011

리눅스 커널 내부구조

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