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