OAKL to the Browser
Working in Rmd, here is getting some DMLs into IGV format. TLDR - bedfile
Methylkit-02
Now I need to download the files
curl -O http://owl.fish.washington.edu/halfshell/bu-serine-wd/18-04-29/zr2096_[1-10]_dedup.sorted.bam
These were generated in Jupyter like:
%%bash
find /Volumes/Serine/wd/18-04-27/zr2096_*R1* \
| xargs basename -s _s1_R1_val_1.fq.gz | xargs -I{} /Applications/bioinfo/Bismark_v0.19.0/bismark \
--path_to_bowtie /Applications/bioinfo/bowtie2-2.3.4.1-macos-x86_64 \
--genome /Volumes/Serine/wd/18-03-15/genome \
--score_min L,0,-1.2 \
-p 4 \
--non_directional \
-1 /Volumes/Serine/wd/18-04-27/{}_s1_R1_val_1.fq.gz \
-2 /Volumes/Serine/wd/18-04-27/{}_s1_R2_val_2.fq.gz \
2> bismark.err
%%bash
/Applications/bioinfo/Bismark_v0.19.0/deduplicate_bismark \
--bam -p \
*.bam \
2> dedup.err
%%bash
find *deduplicated.bam \
| xargs basename -s _s1_R1_val_1_bismark_bt2_pe.deduplicated.bam | xargs -I{} /Applications/bioinfo/samtools-1.3.1/samtools \
sort {}_s1_R1_val_1_bismark_bt2_pe.deduplicated.bam -o {}_dedup.sorted.bam
Reading in some files:
file.list=list( '/Users/sr320/Desktop/oakl/zr2096_1_dedup.sorted.bam',
'/Users/sr320/Desktop/oakl/zr2096_2_dedup.sorted.bam',
'/Users/sr320/Desktop/oakl/zr2096_3_dedup.sorted.bam',
'/Users/sr320/Desktop/oakl/zr2096_4_dedup.sorted.bam',
'/Users/sr320/Desktop/oakl/zr2096_5_dedup.sorted.bam',
'/Users/sr320/Desktop/oakl/zr2096_6_dedup.sorted.bam',
'/Users/sr320/Desktop/oakl/zr2096_7_dedup.sorted.bam',
'/Users/sr320/Desktop/oakl/zr2096_8_dedup.sorted.bam',
'/Users/sr320/Desktop/oakl/zr2096_9_dedup.sorted.bam',
'/Users/sr320/Desktop/oakl/zr2096_10_dedup.sorted.bam'
)
Bringing in bam files, setting context, AND coverage
myobj = processBismarkAln(location = file.list, sample.id = list("1","2","3","4","5","6","7","8","9","10"), assembly = "v3", read.context="CpG", mincov=10, treatment = c(0,0,0,0,0,1,1,1,1,1))
save(myobj, file = "../analyses/myobj_180427files")
load("../analyses/myobj_180427files")
getCoverageStats(myobj[[4]],plot=TRUE,both.strands=FALSE)
getMethylationStats(myobj[[1]],plot=FALSE,both.strands=FALSE)
## methylation statistics per base
## summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 83.33 93.33 85.11 100.00 100.00
## percentiles:
## 0% 10% 20% 30% 40% 50% 60%
## 0.00000 55.55556 80.00000 85.71429 90.47619 93.33333 96.29630
## 70% 80% 90% 95% 99% 99.5% 99.9%
## 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000
## 100%
## 100.00000
getCoverageStats(myobj[[4]],plot=TRUE,both.strands=FALSE)
meth=unite(myobj)
getCorrelation(meth,plot=TRUE)
## 1 2 3 4 5 6 7
## 1 1.0000000 0.4486002 0.4560082 0.4404980 0.4331196 0.4820232 0.4782540
## 2 0.4486002 1.0000000 0.5205769 0.5071152 0.4994732 0.4561182 0.5355796
## 3 0.4560082 0.5205769 1.0000000 0.5078281 0.5042024 0.4610765 0.5485566
## 4 0.4404980 0.5071152 0.5078281 1.0000000 0.4999247 0.4532686 0.5325240
## 5 0.4331196 0.4994732 0.5042024 0.4999247 1.0000000 0.4484705 0.5300285
## 6 0.4820232 0.4561182 0.4610765 0.4532686 0.4484705 1.0000000 0.5713013
## 7 0.4782540 0.5355796 0.5485566 0.5325240 0.5300285 0.5713013 1.0000000
## 8 0.4427913 0.5140576 0.5175793 0.5032337 0.4955144 0.4546524 0.5355500
## 9 0.4579275 0.5216281 0.5311598 0.5135430 0.5025602 0.4660443 0.5469018
## 10 0.4309082 0.5077382 0.5148965 0.5017601 0.4975297 0.4342481 0.5280633
## 8 9 10
## 1 0.4427913 0.4579275 0.4309082
## 2 0.5140576 0.5216281 0.5077382
## 3 0.5175793 0.5311598 0.5148965
## 4 0.5032337 0.5135430 0.5017601
## 5 0.4955144 0.5025602 0.4975297
## 6 0.4546524 0.4660443 0.4342481
## 7 0.5355500 0.5469018 0.5280633
## 8 1.0000000 0.5192701 0.4981650
## 9 0.5192701 1.0000000 0.5222376
## 10 0.4981650 0.5222376 1.0000000