#An example of the analysis run for a data set generated for a D.melanogaster sample by MNase titration #four MNase concentrations were used -- 1.5U, 6.25U, 25U, and 100U #The paired-end sequencing data representing the digestion fragments were aligned to the dm3 genome and stored #in the files test.sample.mn1.5U.bam, test.sample.mn6.25U.bam, test.sample.mn25U.bam, and test.sample.mn100U.bam #in the directory /data/MACC/fly. ## Set the computational environment > library(limma) > library(GenomicRanges) > library(Rsamtools) > library(descr) > library(MACC) > require(parallel) > genome <- "dm3" # only "dm3", "mm9", and "hg19" are currently supported > chrn <- c("chr2L","chr2R","chr3L","chr3R","chr4","chrX","chr2LHet","chr2RHet","chr3LHet","chr3RHet","chrXHet","chrYHet") # use NULL to make calculations for all chromosomes appearing in the input data files > tit.points <- c(1.5,6.25,25,100) > mc.cores <- 12 # the number of cores to use for computations > bin <- 200 # the bin size in base pairs > path2files <- "/data/MACC/fly/" # the path to the data files > setlist <- list( + sample1 = c("test.sample.mn1.5U.bam", "test.sample.mn6.25U.bam", "test.sample.mn25U.bam", "test.sample.mn100U.bam") + ) # the structure describing all the samples to be analyzed > filter.anomalies <- TRUE > scale.genome.to.100Mb <- TRUE > normalize.perBinSize <- TRUE > gen.dat <- read.genomic.data(genome=genome, chrn=chrn, bin=bin, mc.cores=mc.cores) done for 200bp bins > str(gen.dat) List of 4 $ chr.stat:'data.frame': 12 obs. of 3 variables: ..$ chr : chr [1:12] "chr2L" "chr2LHet" "chr2R" "chr2RHet" ... ..$ length : int [1:12] 23011544 368872 21146708 3288761 24543557 2555491 27905053 2517507 1351857 22422827 ... ..$ comment: chr [1:12] "/gbdb/dm3/dm3.2bit" "/gbdb/dm3/dm3.2bit" "/gbdb/dm3/dm3.2bit" "/gbdb/dm3/dm3.2bit" ... $ chrn : Named chr [1:12] "chr2L" "chr2R" "chr3L" "chr3R" ... ..- attr(*, "names")= chr [1:12] "chr2L" "chr2R" "chr3L" "chr3R" ... $ gc.cont :List of 12 ..$ chr2L : num [1:115058] 0.183 0.411 0.401 0.411 0.386 ... ..$ chr2R : num [1:105734] 0.173 0.307 0.277 0.386 0.356 ... .. ... ..$ chrYHet : num [1:1736] 0.208 0.302 0.317 0.421 0.342 ... $ CpG : NULL ## Read-in the data and perform computations > tags <- read.tags(setlist=setlist, path2files=path2files, chrn=gen.dat$chrn, mc.cores=mc.cores, filter.anomalies= filter.anomalies) Reading /store0/jakubm/mnaseseq/package/test.sample.mn1.5U.bam done Reading /store0/jakubm/mnaseseq/package/test.sample.mn6.25U.bam done Reading /store0/jakubm/mnaseseq/package/test.sample.mn25U.bam done Reading /store0/jakubm/mnaseseq/package/test.sample.mn100U.bam done > str(tags) List of 1 $ sample1:List of 4 ..$ test.sample.mn1.5U.bam :List of 12 .. ..$ chr2L : num [1:6077598] 13168425 -13168571 5206729 -5206922 16257619 ... .. ..$ chr2R : num [1:5151088] 5498666 -5498827 13591020 -13591253 15905666 ... .. .. ... .. ..$ chrYHet : num [1:450] 175212 -175701 284209 -284588 151708 ... ..$ test.sample.mn6.25U.bam:List of 12 .. ..$ chr2L : num [1:6891562] 13384937 -13385291 10405460 -10405774 18215294 ... .. ..$ chr2R : num [1:5794382] 17426108 -17426264 10942123 -10942295 402907 ... .. .. ... .. ..$ chrYHet : num [1:488] 114689 -114827 337723 -338062 284164 ... ..$ test.sample.mn25U.bam :List of 12 .. ..$ chr2L : num [1:8007960] 15821701 -15821830 8290271 -8290554 11915653 ... .. ..$ chr2R : num [1:6822558] 9827352 -9827478 5451892 -5452023 11143411 ... .. .. ... .. ..$ chrYHet : num [1:224] 225605 -225833 296203 -296548 4586 ... ..$ test.sample.mn100U.bam :List of 12 .. ..$ chr2L : num [1:8832742] 2854872 -2854998 2098967 -2099081 3239756 ... .. ..$ chr2R : num [1:7801018] 4476030 -4476150 12317772 -12317932 5861755 ... .. .. ... .. ..$ chrYHet : num [1:126] 42326 -42583 79877 -80102 204639 ... > profiles <- generate.profiles(tags,chrn=gen.dat$chrn, mc.cores=mc.cores) done done done done > str(profiles) List of 2 $ frag.pos:List of 1 ..$ sample1:List of 4 .. ..$ test.sample.mn1.5U.bam :List of 12 .. .. ..$ chr2L : num [1:3038799] 13168498 5206802 16257743 17380582 10209288 ... .. .. ..$ chr2R : num [1:2575544] 5498739 13591093 15905738 7041249 2647070 ... .. .. .. ... .. .. ..$ chrYHet : num [1:225] 175628 284282 152049 307152 284244 ... .. ..$ test.sample.mn6.25U.bam:List of 12 .. .. ..$ chr2L : num [1:3445781] 13385218 10405533 18215368 12523591 1521274 ... .. .. ..$ chr2R : num [1:2897191] 17426186 10942222 403174 18246450 10537998 ... .. .. .. ... .. .. ..$ chrYHet : num [1:244] 114758 337989 284237 99979 118552 ... .. ..$ test.sample.mn25U.bam :List of 12 .. .. ..$ chr2L : num [1:4003980] 15821766 8290344 11915727 12666136 16021476 ... .. .. ..$ chr2R : num [1:3411279] 9827415 5451958 11143541 13426362 13632874 ... .. .. .. ... .. .. ..$ chrYHet : num [1:112] 225678 296475 4636 215457 133393 ... .. ..$ test.sample.mn100U.bam :List of 12 .. .. ..$ chr2L : num [1:4416371] 2854935 2099024 3239829 18620408 13385519 ... .. .. ..$ chr2R : num [1:3900509] 4476090 12317852 5861834 10013119 18191068 ... .. .. .. ... .. .. ..$ chrYHet : num [1:63] 42510 79950 204712 215376 114754 ... $ lengths :List of 1 ..$ sample1:List of 4 .. ..$ test.sample.mn1.5U.bam :List of 12 .. .. ..$ chr2L : num [1:3038799] 146 193 197 184 148 233 170 386 232 235 ... .. .. ..$ chr2R : num [1:2575544] 161 233 144 179 144 158 433 169 137 182 ... .. .. .. ... .. .. ..$ chrYHet : num [1:225] 489 379 414 471 418 481 205 433 158 138 ... .. ..$ test.sample.mn6.25U.bam:List of 12 .. .. ..$ chr2L : num [1:3445781] 354 314 147 393 358 159 175 157 279 172 ... .. .. ..$ chr2R : num [1:2897191] 156 172 340 376 168 162 310 163 192 176 ... .. .. .. ... .. .. ..$ chrYHet : num [1:244] 138 339 424 259 495 230 227 363 494 478 ... .. ..$ test.sample.mn25U.bam :List of 12 .. .. ..$ chr2L : num [1:4003980] 129 283 148 153 128 372 172 160 189 178 ... .. .. ..$ chr2R : num [1:3411279] 126 131 203 170 147 124 118 172 256 172 ... .. .. .. ... .. .. ..$ chrYHet : num [1:112] 228 345 99 154 110 386 465 466 369 123 ... .. ..$ test.sample.mn100U.bam :List of 12 .. .. ..$ chr2L : num [1:4416371] 126 114 166 144 100 345 110 149 137 120 ... .. .. ..$ chr2R : num [1:3900509] 120 160 157 144 125 289 98 95 113 128 ... .. .. .. ... .. .. ..$ chrYHet : num [1:63] 257 225 201 119 157 132 100 135 110 446 ... > tags.in.bins <- count.tags.in.bins(profs=profiles$frag.pos, chr.stat=gen.dat$chr.stat, chrn=gen.dat$chrn, mc.cores=mc.cores, bin=bin, scale.genome.to.100Mb=scale.genome.to.100Mb, normalize.perBinSize=normalize.perBinSize) chr2L chr2R chr3L chr3R chr4 chrX chr2LHet chr2RHet chr3LHet chr3RHet chrXHet chrYHet > str(tags.in.bins) List of 1 $ sample1:List of 3 ..$ counts:List of 4 .. ..$ test.sample.mn1.5U.bam :List of 12 .. .. ..$ chr2L : num [1:115058] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..$ chr2R : num [1:105734] 0 0 0 0 0 0 0 0 0 0 ... .. .. .. ... .. .. ..$ chrYHet : num [1:1736] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ test.sample.mn6.25U.bam:List of 12 .. .. ..$ chr2L : num [1:115058] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..$ chr2R : num [1:105734] 0 0 0 0 0 0 0 0 0 0 ... .. .. .. ... .. .. ..$ chrYHet : num [1:1736] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ test.sample.mn25U.bam :List of 12 .. .. ..$ chr2L : num [1:115058] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..$ chr2R : num [1:105734] 0 0 0 0 0 0 0 0 0 0 ... .. .. .. ... .. .. ..$ chrYHet : num [1:1736] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ test.sample.mn100U.bam :List of 12 .. .. ..$ chr2L : num [1:115058] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..$ chr2R : num [1:105734] 0 0 0 0 0 0 0 0 0 0 ... .. .. .. ... .. .. ..$ chrYHet : num [1:1736] 0 0 0 0 0 0 0 0 0 0 ... ..$ pos :List of 12 .. ..$ chr2L : num [1:115058] 100 300 500 700 900 1100 1300 1500 1700 1900 ... .. ..$ chr2R : num [1:105734] 100 300 500 700 900 1100 1300 1500 1700 1900 ... .. .. ... .. ..$ chrYHet : num [1:1736] 100 300 500 700 900 1100 1300 1500 1700 1900 ... ..$ values:List of 4 .. ..$ test.sample.mn1.5U.bam :List of 12 .. .. ..$ chr2L : num [1:115058] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..$ chr2R : num [1:105734] 0 0 0 0 0 0 0 0 0 0 ... .. .. .. ... .. .. ..$ chrYHet : num [1:1736] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ test.sample.mn6.25U.bam:List of 12 .. .. ..$ chr2L : num [1:115058] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..$ chr2R : num [1:105734] 0 0 0 0 0 0 0 0 0 0 ... .. .. .. ... .. .. ..$ chrYHet : num [1:1736] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ test.sample.mn25U.bam :List of 12 .. .. ..$ chr2L : num [1:115058] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..$ chr2R : num [1:105734] 0 0 0 0 0 0 0 0 0 0 ... .. .. .. ... .. .. ..$ chrYHet : num [1:1736] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ test.sample.mn100U.bam :List of 12 .. .. ..$ chr2L : num [1:115058] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..$ chr2R : num [1:105734] 0 0 0 0 0 0 0 0 0 0 ... .. .. ..... .. .. ..$ chrYHet : num [1:1736] 0 0 0 0 0 0 0 0 0 0 ... > MACC.res <- macc(tags.in.bins=tags.in.bins, tit.points=tit.points, gc.cont=gen.dat$gc.cont, chr.stat=gen.dat$chr.stat, chrn=gen.dat$chrn,bin=bin,mc.cores=mc.cores, CpG=gen.dat$CpG) slope for sample1 was computed MACC for sample1 was computed > str(MACC.res) List of 1 $ :'data.frame': 648323 obs. of 10 variables: ..$ chr : Factor w/ 12 levels "chr2L","chr2R",..: 1 1 1 1 1 1 1 1 1 1 ... ..$ start : int [1:648323] 1 201 401 601 801 1001 1201 1401 1601 1801 ... ..$ end : int [1:648323] 200 400 600 800 1000 1200 1400 1600 1800 2000 ... ..$ mn1.5 : num [1:648323] 0 0 0 0 0 0 0 0 0 0 ... ..$ mn6.25: num [1:648323] 0 0 0 0 0 0 0 0 0 0 ... ..$ mn25 : num [1:648323] 0 0 0 0 0 0 0 0 0 0 ... ..$ mn100 : num [1:648323] 0 0 0 0 0 0 0 0 0 0 ... ..$ Pooled: num [1:648323] 0 0 0 0 0 0 0 0 0 0 ... ..$ Slope : num [1:648323] NA NA NA NA NA NA NA NA NA NA ... ..$ MACC : num [1:648323] NA NA NA NA NA NA NA NA NA NA ...