이게 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 사용).
bowtie 같은 프로그램으로 mapping(여기서는 bowtie 사용).
2.
데이터 필터, quality와 duplication read 제거
3.
데이터 필터, 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를 생성하는 것이 가장 보정이 잘될 수 있는거라 그런거 같다.
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
> 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$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를 생성하는 것이 가장 보정이 잘될 수 있는거라 그런거 같다.
