Methylation Main Data QC Analysis
Aaron Kovacs & Mark Ziemann
2021-10-04
This code is available at https://github.com/aaronsk7/guppy-methylation
This script performs differential methylation analysis using MethylKit.
library("R.utils")## Loading required package: R.oo## Loading required package: R.methodsS3## R.methodsS3 v1.8.1 (2020-08-26 16:20:06 UTC) successfully loaded. See ?R.methodsS3 for help.## R.oo v1.24.0 (2020-08-26 16:11:58 UTC) successfully loaded. See ?R.oo for help.##
## Attaching package: 'R.oo'## The following object is masked from 'package:R.methodsS3':
##
## throw## The following objects are masked from 'package:methods':
##
## getClasses, getMethods## The following objects are masked from 'package:base':
##
## attach, detach, load, save## R.utils v2.11.0 (2021-09-26 08:30:02 UTC) successfully loaded. See ?R.utils for help.##
## Attaching package: 'R.utils'## The following object is masked from 'package:utils':
##
## timestamp## The following objects are masked from 'package:base':
##
## cat, commandArgs, getOption, inherits, isOpen, nullfile, parse,
## warningslibrary("parallel")
library("reshape2")
# library("kableExtra") having some trouble installing dependancies on bio2
library("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary("RColorBrewer")
library("GenomicRanges")## Loading required package: stats4## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which, which.max, which.min## Loading required package: S4Vectors##
## Attaching package: 'S4Vectors'## The following objects are masked from 'package:dplyr':
##
## first, rename## The following object is masked from 'package:base':
##
## expand.grid## Loading required package: IRanges##
## Attaching package: 'IRanges'## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice## The following object is masked from 'package:R.oo':
##
## trim## Loading required package: GenomeInfoDblibrary("limma")##
## Attaching package: 'limma'## The following object is masked from 'package:BiocGenerics':
##
## plotMAlibrary("methylKit")##
## Attaching package: 'methylKit'## The following object is masked from 'package:dplyr':
##
## selectlibrary("gplots")##
## Attaching package: 'gplots'## The following object is masked from 'package:IRanges':
##
## space## The following object is masked from 'package:S4Vectors':
##
## space## The following object is masked from 'package:stats':
##
## lowesslibrary("seqinr")##
## Attaching package: 'seqinr'## The following object is masked from 'package:limma':
##
## zscore## The following object is masked from 'package:dplyr':
##
## count## The following object is masked from 'package:R.oo':
##
## getNameImport biscuit data with a function
Define and test the function
vcfz <- "meth_data/test-R1.bam.vcf.gz.cg.bed.gz"
biscuit2methylkit <- function(vcfz) {
vcf <- gsub(".gz$","",vcfz)
gunzip(vcfz)
tdata <- read.table(vcf)
gzip(vcf)
chrBase <- paste(tdata$V1,tdata$V2,sep=".")
chr <- tdata$V1
base <- tdata$V2
strand="F"
coverage <- tdata$V5
freqC <- tdata$V4 *100
freqT <- 100-freqC
mk <- data.frame( chrBase , chr , base, strand, coverage, freqC, freqT )
mk
mkname <- paste(vcf,".mk",sep="")
mkname <- gsub(".bam.vcf.gz.cg.bed","",mkname)
write.table(mk,file=mkname,quote=FALSE,sep="\t",row.names = FALSE)
}
biscuit2methylkit(vcfz)
getwd()## [1] "/home/ziemannm/guppy_methylseq/main_data1"myobj=methRead("meth_data/test-R1.mk",
sample.id="test1",
assembly="hg18",
treatment=1,
context="CpG",
mincov = 10 )## Received single location.## Reading file.myobj## chr start end strand coverage numCs numTs
## 1 KK215283.1 113 113 + 126 3 123
## 2 KK215283.1 2071 2071 + 10 8 2
## 3 KK215283.1 2111 2111 + 13 11 2
## 4 KK215283.1 2139 2139 + 13 13 0
## 5 KK215283.1 4814 4814 + 10 8 2
## 6 KK215283.1 5196 5196 + 10 10 0
## 7 KK215283.1 5241 5241 + 10 7 3
## 8 KK215283.1 5882 5882 + 10 10 0
## 9 KK215283.1 5885 5885 + 10 10 0
## 10 KK215283.1 7324 7324 + 10 9 1
## 11 KK215283.1 8493 8493 + 10 0 10
## 12 KK215283.1 8605 8605 + 10 9 1
## 13 KK215283.1 8617 8617 + 11 10 1
## 14 KK215283.1 8814 8814 + 10 9 1
## 15 KK215283.1 11557 11557 + 15 13 2
## 16 KK215283.1 11638 11638 + 12 12 0
## 17 KK215283.1 12468 12468 + 12 12 0
## 18 KK215283.1 12497 12497 + 10 10 0
## 19 KK215283.1 12506 12506 + 10 9 1
## 20 KK215283.1 12531 12531 + 11 10 1
## 21 KK215283.1 12560 12560 + 11 9 2
## 22 KK215283.1 12575 12575 + 12 10 2
## 23 KK215283.1 12580 12580 + 12 7 5
## 24 KK215283.1 12586 12586 + 12 7 5
## 25 KK215283.1 15071 15071 + 12 8 4
## 26 KK215283.1 16099 16099 + 10 10 0
## 27 KK215283.1 16120 16120 + 11 10 1
## 28 KK215283.1 16149 16149 + 11 11 0
## 29 KK215283.1 16175 16175 + 13 10 3
## 30 KK215283.1 16188 16188 + 10 7 3
## 31 KK215283.1 16647 16647 + 10 8 2
## 32 KK215283.1 16675 16675 + 10 7 3
## 33 KK215283.1 16980 16980 + 10 10 0
## 34 KK215283.1 17746 17746 + 11 9 2
## 35 KK215283.1 17771 17771 + 10 8 2
## 36 KK215283.1 21884 21884 + 10 6 4
## 37 KK215283.1 22800 22800 + 12 12 0
## 38 KK215283.1 26664 26664 + 10 7 3
## 39 KK215283.1 27193 27193 + 11 10 1
## 40 KK215283.1 27205 27205 + 10 9 1
## 41 KK215283.1 27224 27224 + 10 9 1
## 42 KK215283.1 27234 27234 + 10 8 2
## 43 KK215283.1 27978 27978 + 10 9 1
## 44 KK215283.1 29546 29546 + 13 13 0
## 45 KK215283.1 29628 29628 + 11 6 5
## 46 KK215283.1 30232 30232 + 13 12 1
## 47 KK215283.1 30277 30277 + 10 9 1
## 48 KK215283.1 30776 30776 + 10 4 6
## 49 KK215283.1 30909 30909 + 11 9 2
## 50 KK215283.1 31777 31777 + 10 10 0
## 51 KK215283.1 31779 31779 + 10 10 0
## 52 KK215283.1 31950 31950 + 11 11 0
## 53 KK215283.1 31961 31961 + 11 11 0
## 54 KK215283.1 31970 31970 + 10 7 3
## 55 KK215283.1 32001 32001 + 10 10 0
## 56 KK215283.1 32004 32004 + 10 10 0
## 57 KK215283.1 32030 32030 + 10 9 1
## 58 KK215283.1 32038 32038 + 11 11 0
## 59 KK215283.1 32042 32042 + 11 11 0
## 60 KK215283.1 32092 32092 + 10 8 2
## 61 KK215283.1 32101 32101 + 11 10 1
## 62 KK215283.1 32703 32703 + 10 7 3
## 63 KK215283.1 32958 32958 + 10 1 9
## 64 KK215283.1 32963 32963 + 10 3 7
## 65 KK215283.1 32970 32970 + 11 2 9
## 66 KK215283.1 32976 32976 + 10 1 9
## 67 KK215283.1 32979 32979 + 10 3 7
## 68 KK215283.1 32986 32986 + 10 3 7
## 69 KK215283.1 32989 32989 + 10 3 7
## 70 KK215283.1 32994 32994 + 10 2 8
## 71 KK215283.1 33021 33021 + 11 1 10
## 72 KK215283.1 33050 33050 + 11 2 9
## 73 KK215283.1 33081 33081 + 11 1 10
## 74 KK215283.1 33091 33091 + 12 0 12
## 75 KK215283.1 33150 33150 + 10 0 10
## 76 KK215283.1 34458 34458 + 11 10 1
## 77 KK215283.1 34484 34484 + 11 9 2
## 78 KK215283.1 34932 34932 + 10 10 0
## 79 KK215283.1 34937 34937 + 10 10 0
## 80 KK215283.1 37099 37099 + 10 10 0
## 81 KK215283.1 37831 37831 + 13 10 3
## 82 KK215283.1 38700 38700 + 10 9 1
## 83 KK215283.1 38720 38720 + 11 11 0
## 84 KK215283.1 38786 38786 + 11 11 0
## 85 KK215283.1 38856 38856 + 10 10 0
## 86 KK215283.1 39566 39566 + 10 10 0
## 87 KK215283.1 39576 39576 + 11 8 3
## 88 KK215283.1 40597 40597 + 13 4 9
## 89 KK215283.1 41681 41681 + 10 10 0
## 90 KK215283.1 42120 42120 + 10 9 1
## 91 KK215283.1 42134 42134 + 10 7 3
## 92 KK215283.1 42189 42189 + 10 9 1
## 93 KK215283.1 42294 42294 + 10 9 1
## 94 KK215283.1 42309 42309 + 14 13 1
## 95 KK215283.1 42337 42337 + 13 12 1
## 96 KK215283.1 42368 42368 + 13 11 2
## 97 KK215283.1 42371 42371 + 13 12 1
## 98 KK215283.1 42374 42374 + 11 10 1
## 99 KK215283.1 42399 42399 + 11 10 1
## 100 KK215283.1 42404 42404 + 12 10 2
## 101 KK215283.1 42407 42407 + 12 11 1
## 102 KK215283.1 42473 42473 + 12 11 1
## 103 KK215283.1 42504 42504 + 11 11 0
## 104 KK215283.1 42533 42533 + 10 7 3
## 105 KK215283.1 43073 43073 + 10 10 0
## 106 KK215283.1 43091 43091 + 11 11 0
## 107 KK215283.1 43094 43094 + 10 10 0
## 108 KK215283.1 43104 43104 + 10 10 0
## 109 KK215283.1 44555 44555 + 11 9 2
## 110 KK215283.1 44606 44606 + 11 11 0
## 111 KK215283.1 44618 44618 + 10 8 2
## 112 KK215283.1 44694 44694 + 14 0 14
## 113 KK215283.1 44710 44710 + 11 7 4
## 114 KK215283.1 44789 44789 + 10 9 1
## 115 KK215283.1 44990 44990 + 11 10 1
## 116 KK215283.1 45027 45027 + 11 10 1
## 117 KK215283.1 45074 45074 + 11 11 0
## 118 KK215283.1 45377 45377 + 10 7 3
## 119 KK215283.1 45391 45391 + 10 8 2
## 120 KK215283.1 46352 46352 + 11 5 6
## 121 KK215283.1 46403 46403 + 10 5 5
## 122 KK215283.1 46459 46459 + 10 0 10
## 123 KK215283.1 47013 47013 + 10 5 5
## 124 KK215283.1 47586 47586 + 10 7 3
## 125 KK215283.1 47672 47672 + 10 6 4
## 126 KK215283.1 47675 47675 + 10 6 4
## 127 KK215283.1 47677 47677 + 10 6 4
## 128 KK215283.1 47685 47685 + 10 6 4
## 129 KK215283.1 47690 47690 + 10 3 7
## 130 KK215283.1 47876 47876 + 10 10 0
## 131 KK215283.1 48080 48080 + 10 6 4
## 132 KK215283.1 49229 49229 + 10 6 4
## 133 KK215283.1 55651 55651 + 10 0 10
## 134 KK215283.1 56955 56955 + 10 8 2
## 135 KK215283.1 58768 58768 + 10 1 9
## 136 KK215283.1 58770 58770 + 10 0 10
## 137 KK215283.1 59050 59050 + 11 0 11
## 138 KK215283.1 59058 59058 + 10 0 10
## 139 KK215283.1 59363 59363 + 11 0 11
## 140 KK215283.1 62888 62888 + 13 13 0
## 141 KK215283.1 62927 62927 + 11 11 0
## 142 KK215283.1 63931 63931 + 12 12 0
## 143 KK215283.1 65776 65776 + 10 10 0
## 144 KK215283.1 65860 65860 + 13 11 2
## 145 KK215283.1 65896 65896 + 11 9 2
## 146 KK215283.1 65905 65905 + 11 9 2
## 147 KK215283.1 65912 65912 + 11 9 2
## 148 KK215283.1 66304 66304 + 10 10 0
## 149 KK215283.1 66963 66963 + 10 9 1
## 150 KK215283.1 67009 67009 + 12 0 12
## 151 KK215283.1 67014 67014 + 12 12 0
## 152 KK215283.1 67048 67048 + 14 14 0
## 153 KK215283.1 67118 67118 + 11 11 0
## 154 KK215283.1 67139 67139 + 10 8 2
## 155 KK215283.1 67566 67566 + 10 9 1
## 156 KK215283.1 67599 67599 + 11 10 1
## 157 KK215283.1 67605 67605 + 11 10 1
## 158 KK215283.1 67610 67610 + 11 9 2
## 159 KK215283.1 67614 67614 + 11 10 1
## 160 KK215283.1 67619 67619 + 10 8 2
## 161 KK215283.1 67623 67623 + 10 8 2
## 162 KK215283.1 67627 67627 + 10 8 2
## 163 KK215283.1 68352 68352 + 11 10 1
## 164 KK215283.1 68354 68354 + 11 10 1
## 165 KK215283.1 69173 69173 + 12 12 0
## 166 KK215283.1 69176 69176 + 12 10 2
## 167 KK215283.1 69202 69202 + 11 10 1
## 168 KK215283.1 69228 69228 + 10 10 0
## 169 KK215283.1 70298 70298 + 10 8 2
## 170 KK215283.1 70918 70918 + 12 10 2
## 171 KK215283.1 71266 71266 + 11 10 1
## 172 KK215283.1 71441 71441 + 10 9 1
## 173 KK215283.1 71472 71472 + 11 9 2
## 174 KK215283.1 71961 71961 + 11 11 0
## 175 KK215283.1 71989 71989 + 12 12 0
## 176 KK215283.1 71996 71996 + 12 12 0
## 177 KK215283.1 72001 72001 + 11 11 0
## 178 KK215283.1 72041 72041 + 10 10 0
## 179 KK215283.1 72057 72057 + 12 12 0
## 180 KK215283.1 72062 72062 + 10 10 0
## 181 KK215283.1 72863 72863 + 10 9 1
## 182 KK215283.1 74099 74099 + 11 11 0
## 183 KK215283.1 78570 78570 + 10 7 3
## 184 KK215283.1 78600 78600 + 10 6 4
## 185 KK215283.1 78650 78650 + 10 8 2
## 186 KK215283.1 78659 78659 + 12 9 3
## 187 KK215283.1 78675 78675 + 12 9 3
## 188 KK215283.1 78687 78687 + 10 7 3
## 189 KK215283.1 78823 78823 + 10 9 1
## 190 KK215283.1 78981 78981 + 10 0 10
## 191 KK215283.1 78991 78991 + 11 9 2
## 192 KK215283.1 79162 79162 + 10 10 0
## 193 KK215283.1 81657 81657 + 10 10 0
## 194 KK215283.1 81665 81665 + 10 10 0
## 195 KK215283.1 81717 81717 + 10 9 1
## 196 KK215283.1 82702 82702 + 11 9 2
## 197 KK215283.1 82727 82727 + 13 11 2
## 198 KK215283.1 82743 82743 + 12 9 3
## 199 KK215283.1 82745 82745 + 10 8 2
## 200 KK215283.1 82816 82816 + 10 8 2
## 201 KK215283.1 82929 82929 + 10 6 4
## 202 KK215283.1 82956 82956 + 12 7 5
## 203 KK215283.1 87012 87012 + 10 6 4
## 204 KK215283.1 87126 87126 + 10 8 2
## 205 KK215283.1 87238 87238 + 10 7 3
## 206 KK215283.1 87240 87240 + 10 8 2
## 207 KK215283.1 87250 87250 + 12 10 2
## 208 KK215283.1 87265 87265 + 13 9 4
## 209 KK215283.1 87325 87325 + 12 11 1
## 210 KK215283.1 87330 87330 + 12 12 0
## 211 KK215283.1 87340 87340 + 11 9 2
## 212 KK215283.1 87345 87345 + 11 11 0
## 213 KK215283.1 87356 87356 + 10 10 0
## 214 KK215283.1 87962 87962 + 11 11 0
## 215 KK215283.1 88241 88241 + 11 9 2
## 216 KK215283.1 89147 89147 + 10 7 3
## 217 KK215283.1 89152 89152 + 11 7 4
## 218 KK215283.1 89172 89172 + 11 8 3
## 219 KK215283.1 89736 89736 + 10 2 8
## 220 KK215283.1 89792 89792 + 10 2 8
## 221 KK215283.1 89800 89800 + 11 4 7
## 222 KK215283.1 89805 89805 + 11 5 6
## 223 KK215283.1 89814 89814 + 12 4 8
## 224 KK215283.1 89816 89816 + 12 6 6
## 225 KK215283.1 89965 89965 + 10 6 4
## 226 KK215283.1 90163 90163 + 11 5 6
## 227 KK215283.1 90172 90172 + 10 7 3
## 228 KK215283.1 94049 94049 + 11 7 4
## 229 KK215283.1 97557 97557 + 10 8 2
## 230 KK215283.1 104649 104649 + 10 10 0
## 231 KK215283.1 104656 104656 + 11 11 0
## 232 KK215283.1 104668 104668 + 11 10 1
## 233 KK215283.1 105036 105036 + 11 10 1
## 234 KK215283.1 105061 105061 + 11 10 1
## 235 KK215283.1 105610 105610 + 10 8 2
## 236 KK215283.1 105633 105633 + 11 9 2
## 237 KK215283.1 105657 105657 + 10 10 0
## 238 KK215283.1 108913 108913 + 11 10 1
## 239 KK215283.1 108937 108937 + 10 10 0
## 240 KK215283.1 108963 108963 + 11 11 0
## 241 KK215283.1 109818 109818 + 10 6 4
## 242 KK215283.1 110802 110802 + 10 7 3
## 243 KK215283.1 112149 112149 + 10 9 1
## 244 KK215283.1 112850 112850 + 11 7 4
## 245 KK215283.1 113246 113246 + 11 10 1
## 246 KK215283.1 113252 113252 + 11 10 1
## 247 KK215283.1 113613 113613 + 10 9 1
## 248 KK215283.1 113614 113614 + 10 7 3
## 249 KK215283.1 113619 113619 + 10 6 4
## 250 KK215283.1 114345 114345 + 12 10 2
## 251 KK215283.1 114794 114794 + 10 9 1
## 252 KK215283.1 116970 116970 + 10 9 1
## 253 KK215283.1 118884 118884 + 12 11 1
## 254 KK215283.1 118891 118891 + 10 10 0
## 255 KK215283.1 119085 119085 + 10 9 1
## 256 KK215283.1 119098 119098 + 10 9 1
## 257 KK215283.1 127585 127585 + 10 10 0
## 258 KK215283.1 131561 131561 + 10 9 1
## 259 KK215283.1 134319 134319 + 11 10 1
## 260 KK215283.1 134366 134366 + 12 9 3
## 261 KK215283.1 134373 134373 + 10 7 3
## 262 KK215283.1 135191 135191 + 13 5 8
## 263 KK215283.1 135194 135194 + 11 6 5
## 264 KK215283.1 135199 135199 + 13 8 5
## 265 KK215283.1 135217 135217 + 12 6 6
## 266 KK215283.1 135221 135221 + 11 6 5
## 267 KK215283.1 135224 135224 + 10 7 3
## 268 KK215283.1 135225 135225 + 12 8 4
## 269 KK215283.1 135256 135256 + 10 3 7
## 270 KK215283.1 135259 135259 + 10 7 3
## 271 KK215283.1 135287 135287 + 11 6 5
## 272 KK215283.1 135288 135288 + 10 2 8
## 273 KK215283.1 135294 135294 + 11 7 4
## 274 KK215283.1 135295 135295 + 11 7 4
## 275 KK215283.1 135299 135299 + 10 5 5
## 276 KK215283.1 135300 135300 + 10 5 5
## 277 KK215283.1 135305 135305 + 10 5 5
## 278 KK215283.1 135313 135313 + 10 4 6
## 279 KK215283.1 135328 135328 + 10 2 8
## 280 KK215283.1 135831 135831 + 10 4 6
## 281 KK215283.1 135914 135914 + 11 9 2
## 282 KK215283.1 137237 137237 + 10 10 0
## 283 KK215283.1 138047 138047 + 10 9 1
## 284 KK215283.1 138054 138054 + 10 9 1
## 285 KK215283.1 138059 138059 + 10 8 2
## 286 KK215283.1 138631 138631 + 10 6 4
## 287 KK215283.1 138667 138667 + 10 6 4
## 288 KK215283.1 138672 138672 + 10 3 7
## 289 KK215283.1 139876 139876 + 10 10 0
## 290 KK215283.1 139882 139882 + 10 10 0
## 291 KK215283.1 139893 139893 + 10 9 1
## 292 KK215283.1 140049 140049 + 10 10 0
## 293 KK215283.1 140051 140051 + 10 10 0
## 294 KK215283.1 140061 140061 + 11 11 0
## 295 KK215283.1 140064 140064 + 11 11 0
## 296 KK215283.1 140066 140066 + 11 7 4
## 297 KK215283.1 140076 140076 + 11 11 0
## 298 KK215283.1 140079 140079 + 11 11 0
## 299 KK215283.1 140081 140081 + 11 11 0
## 300 KK215283.1 140091 140091 + 11 7 4
## 301 KK215283.1 140096 140096 + 11 7 4
## 302 KK215283.1 140228 140228 + 10 9 1
## 303 KK215283.1 140858 140858 + 11 9 2
## 304 KK215283.1 140888 140888 + 10 9 1
## 305 KK215283.1 140902 140902 + 11 10 1
## 306 KK215283.1 140924 140924 + 14 13 1
## 307 KK215283.1 140980 140980 + 11 11 0
## 308 KK215283.1 141012 141012 + 11 11 0
## 309 KK215283.1 141086 141086 + 10 10 0
## 310 KK215283.1 141160 141160 + 11 10 1
## 311 KK215283.1 141170 141170 + 12 12 0
## 312 KK215283.1 141177 141177 + 11 10 1
## 313 KK215283.1 141184 141184 + 10 10 0
## 314 KK215283.1 141192 141192 + 10 10 0
## 315 KK215283.1 141198 141198 + 10 10 0
## 316 KK215283.1 141206 141206 + 11 11 0
## 317 KK215283.1 141211 141211 + 11 11 0
## 318 KK215283.1 142124 142124 + 10 5 5
## 319 KK215283.1 142156 142156 + 10 3 7
## 320 KK215283.1 142171 142171 + 11 5 6
## 321 KK215283.1 144017 144017 + 10 9 1
## 322 KK215283.1 144088 144088 + 10 10 0
## 323 KK215283.1 145044 145044 + 12 12 0
## 324 KK215283.1 145089 145089 + 12 10 2
## 325 KK215283.1 145104 145104 + 12 10 2
## 326 KK215283.1 147395 147395 + 10 3 7
## 327 KK215283.1 147436 147436 + 11 5 6
## 328 KK215283.1 148438 148438 + 10 9 1
## 329 KK215283.1 148452 148452 + 10 10 0
## 330 KK215283.1 148455 148455 + 10 10 0
## 331 KK215283.1 151569 151569 + 10 3 7
## 332 KK215283.1 152962 152962 + 10 10 0
## 333 KK215283.1 156476 156476 + 11 8 3
## 334 KK215283.1 156481 156481 + 11 8 3
## 335 KK215283.1 156503 156503 + 12 7 5
## 336 KK215283.1 156518 156518 + 10 7 3
## 337 KK215283.1 156882 156882 + 10 8 2
## 338 KK215283.1 156886 156886 + 10 6 4
## 339 KK215283.1 156938 156938 + 12 7 5
## 340 KK215283.1 156951 156951 + 11 9 2
## 341 KK215283.1 156967 156967 + 11 5 6
## 342 KK215283.1 156978 156978 + 12 6 6
## 343 KK215283.1 156980 156980 + 11 7 4
## 344 KK215283.1 157180 157180 + 10 8 2
## 345 KK215283.1 157183 157183 + 10 6 4
## 346 KK215283.1 157281 157281 + 10 5 5
## 347 KK215283.1 157285 157285 + 10 5 5
## 348 KK215283.1 160536 160536 + 10 7 3
## 349 KK215283.1 160576 160576 + 13 8 5
## 350 KK215283.1 160596 160596 + 11 7 4
## 351 KK215283.1 160626 160626 + 15 13 2
## 352 KK215283.1 160636 160636 + 15 6 9
## 353 KK215283.1 160646 160646 + 13 10 3
## 354 KK215283.1 162729 162729 + 10 6 4
## 355 KK215283.1 163100 163100 + 10 10 0
## 356 KK215283.1 164088 164088 + 10 9 1
## 357 KK215283.1 164100 164100 + 10 8 2
## 358 KK215283.1 164160 164160 + 10 9 1
## 359 KK215283.1 165919 165919 + 10 7 3
## 360 KK215283.1 165928 165928 + 10 5 5
## 361 KK215283.1 170217 170217 + 10 10 0
## 362 KK215283.1 170235 170235 + 11 11 0
## 363 KK215283.1 170244 170244 + 11 11 0
## 364 KK215283.1 170247 170247 + 11 11 0
## 365 KK215283.1 170255 170255 + 10 10 0
## 366 KK215283.1 170349 170349 + 12 12 0
## 367 KK215283.1 170357 170357 + 12 11 1
## 368 KK215283.1 170360 170360 + 12 12 0
## 369 KK215283.1 170366 170366 + 11 11 0
## 370 KK215283.1 170956 170956 + 10 1 9
## 371 KK215283.1 173240 173240 + 13 4 9
## 372 KK215283.1 173258 173258 + 13 7 6
## 373 KK215283.1 173270 173270 + 13 5 8
## 374 KK215283.1 173296 173296 + 13 6 7
## 375 KK215283.1 173311 173311 + 11 4 7
## 376 KK215283.1 173315 173315 + 12 5 7
## 377 KK215283.1 173327 173327 + 11 6 5
## 378 KK215283.1 173338 173338 + 11 7 4
## 379 KK215283.1 173341 173341 + 11 7 4
## 380 KK215283.1 173346 173346 + 10 6 4
## 381 KK215283.1 173349 173349 + 10 7 3
## 382 KK215283.1 173355 173355 + 11 7 4
## 383 KK215283.1 173675 173675 + 10 8 2
## 384 KK215283.1 174662 174662 + 10 9 1
## 385 KK215283.1 174749 174749 + 10 10 0
## 386 KK215283.1 177479 177479 + 10 6 4
## 387 KK215283.1 178687 178687 + 10 7 3
## 388 KK215283.1 178696 178696 + 11 7 4
## 389 KK215283.1 182396 182396 + 11 7 4
## 390 KK215283.1 182423 182423 + 10 6 4
## 391 KK215283.1 182886 182886 + 10 5 5
## 392 KK215283.1 182890 182890 + 10 7 3
## 393 KK215283.1 183275 183275 + 10 3 7
## 394 KK215283.1 187169 187169 + 10 10 0
## 395 KK215283.1 192320 192320 + 11 4 7
## 396 KK215283.1 192351 192351 + 11 6 5
## 397 KK215283.1 192396 192396 + 13 4 9
## 398 KK215283.1 192432 192432 + 12 5 7
## 399 KK215283.1 192441 192441 + 11 4 7
## 400 KK215283.1 192449 192449 + 13 7 6
## 401 KK215283.1 193447 193447 + 11 8 3
## 402 KK215283.1 193738 193738 + 10 6 4
## 403 KK215283.1 197639 197639 + 10 9 1
## 404 KK215283.1 197640 197640 + 10 9 1
## 405 KK215283.1 197644 197644 + 10 10 0
## 406 KK215283.1 197645 197645 + 11 10 1
## 407 KK215283.1 197674 197674 + 13 12 1
## 408 KK215283.1 197678 197678 + 13 13 0
## 409 KK215283.1 197683 197683 + 12 12 0
## 410 KK215283.1 197685 197685 + 12 12 0
## 411 KK215283.1 197697 197697 + 11 10 1
## 412 KK215283.1 197704 197704 + 10 10 0
## 413 KK215283.1 197720 197720 + 10 10 0
## 414 KK215283.1 199955 199955 + 10 0 10
## 415 KK215283.1 199959 199959 + 10 4 6
## 416 KK215283.1 200127 200127 + 10 1 9
## 417 KK215283.1 200180 200180 + 10 1 9
## 418 KK215283.1 200215 200215 + 10 0 10
## 419 KK215283.1 200587 200587 + 11 6 5
## 420 KK215283.1 200628 200628 + 12 7 5
## 421 KK215283.1 200639 200639 + 11 9 2
## 422 KK215283.1 200663 200663 + 11 8 3
## 423 KK215283.1 201801 201801 + 10 9 1
## 424 KK215283.1 201805 201805 + 10 9 1
## 425 KK215283.1 203495 203495 + 12 11 1
## 426 KK215283.1 203500 203500 + 11 10 1
## 427 KK215283.1 203570 203570 + 10 5 5
## 428 KK215283.1 203645 203645 + 11 11 0
## 429 KK215283.1 203724 203724 + 11 11 0
## 430 KK215283.1 203777 203777 + 10 9 1
## 431 KK215283.1 204841 204841 + 10 6 4
## 432 KK215283.1 206562 206562 + 10 1 9
## 433 KK215283.1 206637 206637 + 10 0 10
## 434 KK215283.1 206678 206678 + 11 1 10
## 435 KK215283.1 206681 206681 + 12 0 12
## 436 KK215283.1 206686 206686 + 12 0 12
## 437 KK215283.1 206703 206703 + 10 0 10
## 438 KK215283.1 206712 206712 + 11 1 10
## 439 KK215283.1 206719 206719 + 10 1 9
## 440 KK215283.1 206727 206727 + 10 0 10
## 441 KK215283.1 208819 208819 + 11 10 1
## 442 KK215283.1 208839 208839 + 10 7 3
## 443 KK215283.1 208901 208901 + 10 0 10
## 444 KK215283.1 208911 208911 + 10 5 5
## 445 KK215283.1 209476 209476 + 12 9 3
## 446 KK215283.1 209577 209577 + 10 7 3
## 447 KK215283.1 209602 209602 + 11 10 1
## 448 KK215283.1 209608 209608 + 12 11 1Now format all the data into methylkit format
This might take a few minutes - get a cuppa.
basenames <- c("Clear2F-R1","Clear2F-R2","Clear2F-R3","Foundation","Green3F-R1","Green3F-R2","Green3F-R3",
"Lilac4F-R1","Lilac4F-R2","Lilac4F-R3")
for (i in basenames){
bedgz <- gsub("$","\\.bam\\.vcf\\.gz\\.cg\\.bed\\.gz",i)
bedgz <- paste("meth_data/",bedgz,sep = "")
mk <- gsub("$","\\.mk",i)
mk <- paste("meth_data/",mk,sep = "")
message(mk)
if(file.exists(mk)){
message(paste(mk,"exists"))
}else{
biscuit2methylkit(vcfz = bedgz)
}}## meth_data/Clear2F-R1.mk## meth_data/Clear2F-R1.mk exists## meth_data/Clear2F-R2.mk## meth_data/Clear2F-R2.mk exists## meth_data/Clear2F-R3.mk## meth_data/Clear2F-R3.mk exists## meth_data/Foundation.mk## meth_data/Foundation.mk exists## meth_data/Green3F-R1.mk## meth_data/Green3F-R1.mk exists## meth_data/Green3F-R2.mk## meth_data/Green3F-R2.mk exists## meth_data/Green3F-R3.mk## meth_data/Green3F-R3.mk exists## meth_data/Lilac4F-R1.mk## meth_data/Lilac4F-R1.mk exists## meth_data/Lilac4F-R2.mk## meth_data/Lilac4F-R2.mk exists## meth_data/Lilac4F-R3.mk## meth_data/Lilac4F-R3.mk existsmyfiles2 <- list.files("meth_data",pattern = "mk$",full.names = TRUE)
myfiles2## [1] "meth_data/Clear2F-R1.mk" "meth_data/Clear2F-R2.mk"
## [3] "meth_data/Clear2F-R3.mk" "meth_data/Foundation.mk"
## [5] "meth_data/Green3F-R1.mk" "meth_data/Green3F-R2.mk"
## [7] "meth_data/Green3F-R3.mk" "meth_data/Lilac4F-R1.mk"
## [9] "meth_data/Lilac4F-R2.mk" "meth_data/Lilac4F-R3.mk"
## [11] "meth_data/test-R1.mk"AllSamples <- as.list(myfiles2[c(grep("Clear",myfiles2) , grep("Green",myfiles2), grep("Lilac",myfiles2), grep("Foundation",myfiles2))])
AllSamples## [[1]]
## [1] "meth_data/Clear2F-R1.mk"
##
## [[2]]
## [1] "meth_data/Clear2F-R2.mk"
##
## [[3]]
## [1] "meth_data/Clear2F-R3.mk"
##
## [[4]]
## [1] "meth_data/Green3F-R1.mk"
##
## [[5]]
## [1] "meth_data/Green3F-R2.mk"
##
## [[6]]
## [1] "meth_data/Green3F-R3.mk"
##
## [[7]]
## [1] "meth_data/Lilac4F-R1.mk"
##
## [[8]]
## [1] "meth_data/Lilac4F-R2.mk"
##
## [[9]]
## [1] "meth_data/Lilac4F-R3.mk"
##
## [[10]]
## [1] "meth_data/Foundation.mk"myobj <- methRead(AllSamples,
sample.id = list("clear1","clear2","clear3","green1","green2","green3","Lilac1","Lilac2","Lilac3","Foundation"),
assembly = "hg19",
treatment = c(0,0,0,0,0,1,1,1,1,1),
context = "CpG",
mincov = 3 )## Received list of locations.## Reading file.
## Reading file.
## Reading file.
## Reading file.
## Reading file.
## Reading file.
## Reading file.
## Reading file.
## Reading file.
## Reading file.Simple stats
And prepare for differential methylation analysis.
lapply(myobj,function(x){
getMethylationStats(x,plot=FALSE,both.strands=FALSE)
getMethylationStats(x,plot=TRUE,both.strands=FALSE)
getCoverageStats(x,plot=TRUE,both.strands=FALSE)
})## methylation statistics per base
## summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 60.00 83.33 73.34 100.00 100.00
## percentiles:
## 0% 10% 20% 30% 40% 50% 60% 70%
## 0.00000 16.66667 50.00000 66.66667 75.00000 83.33333 90.90909 100.00000
## 80% 90% 95% 99% 99.5% 99.9% 100%
## 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000
## methylation statistics per base
## summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 63.64 84.21 73.47 94.44 100.00
## percentiles:
## 0% 10% 20% 30% 40% 50% 60% 70%
## 0.00000 20.00000 54.54545 70.00000 77.77778 84.21053 88.88889 92.30769
## 80% 90% 95% 99% 99.5% 99.9% 100%
## 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000
## methylation statistics per base
## summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 62.50 83.33 72.98 93.33 100.00
## percentiles:
## 0% 10% 20% 30% 40% 50% 60% 70%
## 0.00000 20.00000 54.54545 68.75000 77.77778 83.33333 88.23529 91.66667
## 80% 90% 95% 99% 99.5% 99.9% 100%
## 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000
## methylation statistics per base
## summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 63.64 83.33 72.87 92.86 100.00
## percentiles:
## 0% 10% 20% 30% 40% 50% 60% 70%
## 0.00000 20.00000 54.54545 69.23077 77.77778 83.33333 87.50000 91.66667
## 80% 90% 95% 99% 99.5% 99.9% 100%
## 94.73684 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000
## methylation statistics per base
## summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 62.50 83.33 73.35 94.12 100.00
## percentiles:
## 0% 10% 20% 30% 40% 50% 60% 70%
## 0.00000 20.00000 54.54545 69.23077 77.77778 83.33333 88.88889 92.30769
## 80% 90% 95% 99% 99.5% 99.9% 100%
## 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000
## methylation statistics per base
## summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 62.50 82.14 72.06 92.31 100.00
## percentiles:
## 0% 10% 20% 30% 40% 50% 60% 70%
## 0.00000 20.00000 54.54545 68.75000 76.47059 82.14286 86.66667 90.47619
## 80% 90% 95% 99% 99.5% 99.9% 100%
## 93.75000 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000
## methylation statistics per base
## summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 63.64 84.21 73.39 93.33 100.00
## percentiles:
## 0% 10% 20% 30% 40% 50% 60% 70%
## 0.00000 20.00000 55.55556 70.00000 78.57143 84.21053 88.88889 92.30769
## 80% 90% 95% 99% 99.5% 99.9% 100%
## 96.15385 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000
## methylation statistics per base
## summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 61.54 83.33 72.36 92.86 100.00
## percentiles:
## 0% 10% 20% 30% 40% 50% 60% 70%
## 0.00000 18.75000 53.33333 66.66667 76.92308 83.33333 87.50000 90.90909
## 80% 90% 95% 99% 99.5% 99.9% 100%
## 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000
## methylation statistics per base
## summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 62.50 83.33 73.18 100.00 100.00
## percentiles:
## 0% 10% 20% 30% 40% 50% 60% 70%
## 0.00000 16.66667 50.00000 66.66667 77.77778 83.33333 88.88889 94.44444
## 80% 90% 95% 99% 99.5% 99.9% 100%
## 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000
## methylation statistics per base
## summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 68.18 88.00 75.45 95.45 100.00
## percentiles:
## 0% 10% 20% 30% 40% 50% 60% 70%
## 0.00000 20.00000 57.69231 75.00000 83.33333 88.00000 91.30435 94.44444
## 80% 90% 95% 99% 99.5% 99.9% 100%
## 96.15385 100.00000 100.00000 100.00000 100.00000 100.00000 100.00000
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
##
## [[5]]
## NULL
##
## [[6]]
## NULL
##
## [[7]]
## NULL
##
## [[8]]
## NULL
##
## [[9]]
## NULL
##
## [[10]]
## NULLmyobj <- filterByCoverage(myobj,lo.count=1,lo.perc=NULL,
hi.count=NULL,hi.perc=99.9)
meth <- unite(myobj, destrand=FALSE)## uniting...remove(myobj)
head(meth)## chr start end strand coverage1 numCs1 numTs1 coverage2 numCs2 numTs2
## 1 KK215283.1 652 652 + 3 3 0 7 7 0
## 2 KK215283.1 772 772 + 3 3 0 9 7 2
## 3 KK215283.1 775 775 + 4 2 2 9 7 2
## 4 KK215283.1 846 846 + 4 4 0 7 6 1
## 5 KK215283.1 1662 1662 + 3 2 1 11 10 1
## 6 KK215283.1 1693 1693 + 4 2 2 7 7 0
## coverage3 numCs3 numTs3 coverage4 numCs4 numTs4 coverage5 numCs5 numTs5
## 1 10 8 2 16 15 1 7 5 2
## 2 11 11 0 12 10 2 7 7 0
## 3 10 9 1 12 8 4 6 6 0
## 4 4 3 1 9 9 0 7 7 0
## 5 3 2 1 9 6 3 5 5 0
## 6 4 3 1 12 11 1 4 4 0
## coverage6 numCs6 numTs6 coverage7 numCs7 numTs7 coverage8 numCs8 numTs8
## 1 16 14 2 12 12 0 13 11 2
## 2 18 17 1 15 15 0 11 9 2
## 3 18 17 1 15 14 1 12 10 2
## 4 13 9 4 8 7 1 4 3 1
## 5 8 8 0 7 6 1 6 5 1
## 6 8 7 1 9 5 4 5 4 1
## coverage9 numCs9 numTs9 coverage10 numCs10 numTs10
## 1 8 8 0 16 15 1
## 2 11 10 1 10 9 1
## 3 10 8 2 8 6 2
## 4 3 3 0 11 11 0
## 5 4 4 0 12 11 1
## 6 5 4 1 12 12 0# plots
getCorrelation(meth,plot=TRUE)## clear1 clear2 clear3 green1 green2 green3
## clear1 1.0000000 0.8010335 0.7972176 0.8089492 0.8004106 0.8102666
## clear2 0.8010335 1.0000000 0.8514886 0.8636977 0.8549608 0.8652150
## clear3 0.7972176 0.8514886 1.0000000 0.8594294 0.8506624 0.8617598
## green1 0.8089492 0.8636977 0.8594294 1.0000000 0.8631310 0.8739499
## green2 0.8004106 0.8549608 0.8506624 0.8631310 1.0000000 0.8647369
## green3 0.8102666 0.8652150 0.8617598 0.8739499 0.8647369 1.0000000
## Lilac1 0.8100390 0.8646971 0.8612553 0.8733649 0.8642651 0.8753751
## Lilac2 0.7959160 0.8496160 0.8472358 0.8580093 0.8493268 0.8606130
## Lilac3 0.7811309 0.8341844 0.8299141 0.8427774 0.8334585 0.8438452
## Foundation 0.8236716 0.8803310 0.8714570 0.8894215 0.8790812 0.8872556
## Lilac1 Lilac2 Lilac3 Foundation
## clear1 0.8100390 0.7959160 0.7811309 0.8236716
## clear2 0.8646971 0.8496160 0.8341844 0.8803310
## clear3 0.8612553 0.8472358 0.8299141 0.8714570
## green1 0.8733649 0.8580093 0.8427774 0.8894215
## green2 0.8642651 0.8493268 0.8334585 0.8790812
## green3 0.8753751 0.8606130 0.8438452 0.8872556
## Lilac1 1.0000000 0.8602507 0.8432553 0.8865881
## Lilac2 0.8602507 1.0000000 0.8283482 0.8654853
## Lilac3 0.8432553 0.8283482 1.0000000 0.8595416
## Foundation 0.8865881 0.8654853 0.8595416 1.0000000
clusterSamples(meth, dist="correlation", method="ward", plot=TRUE)## The "ward" method has been renamed to "ward.D"; note new "ward.D2"
##
## Call:
## hclust(d = d, method = HCLUST.METHODS[hclust.method])
##
## Cluster method : ward.D
## Distance : pearson
## Number of objects: 10hc = clusterSamples(meth, dist="correlation", method="ward", plot=FALSE)## The "ward" method has been renamed to "ward.D"; note new "ward.D2"PCASamples(meth, screeplot=TRUE)
PCASamples(meth)
Session information
For reproducibility
sessionInfo()## R version 4.0.1 (2020-06-06)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS/LAPACK: /ceph-g/opt/miniconda3/envs/R40/lib/libopenblasp-r0.3.9.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] seqinr_4.2-8 gplots_3.1.1 methylKit_1.14.2
## [4] limma_3.44.3 GenomicRanges_1.40.0 GenomeInfoDb_1.24.2
## [7] IRanges_2.22.2 S4Vectors_0.26.1 BiocGenerics_0.34.0
## [10] RColorBrewer_1.1-2 dplyr_1.0.0 reshape2_1.4.4
## [13] R.utils_2.11.0 R.oo_1.24.0 R.methodsS3_1.8.1
##
## loaded via a namespace (and not attached):
## [1] mclust_5.4.7 Rcpp_1.0.7
## [3] bdsmatrix_1.3-4 mvtnorm_1.1-2
## [5] lattice_0.20-41 Rsamtools_2.4.0
## [7] Biostrings_2.56.0 gtools_3.9.2
## [9] digest_0.6.25 R6_2.4.1
## [11] plyr_1.8.6 emdbook_1.3.12
## [13] evaluate_0.14 coda_0.19-4
## [15] ggplot2_3.3.5 pillar_1.4.4
## [17] zlibbioc_1.34.0 rlang_0.4.11
## [19] data.table_1.12.8 Matrix_1.2-18
## [21] bbmle_1.0.24 qvalue_2.20.0
## [23] rmarkdown_2.3 splines_4.0.1
## [25] BiocParallel_1.22.0 stringr_1.4.0
## [27] fastseg_1.34.0 RCurl_1.98-1.5
## [29] munsell_0.5.0 DelayedArray_0.14.1
## [31] compiler_4.0.1 numDeriv_2016.8-1.1
## [33] rtracklayer_1.48.0 xfun_0.15
## [35] pkgconfig_2.0.3 htmltools_0.5.0
## [37] SummarizedExperiment_1.18.2 tidyselect_1.1.0
## [39] tibble_3.0.1 GenomeInfoDbData_1.2.3
## [41] matrixStats_0.61.0 XML_3.99-0.8
## [43] crayon_1.3.4 GenomicAlignments_1.24.0
## [45] MASS_7.3-51.6 bitops_1.0-7
## [47] grid_4.0.1 gtable_0.3.0
## [49] lifecycle_0.2.0 magrittr_1.5
## [51] scales_1.1.1 KernSmooth_2.23-17
## [53] stringi_1.4.6 XVector_0.28.0
## [55] ellipsis_0.3.1 generics_0.0.2
## [57] vctrs_0.3.1 tools_4.0.1
## [59] ade4_1.7-18 Biobase_2.48.0
## [61] glue_1.4.1 purrr_0.3.4
## [63] yaml_2.2.1 colorspace_1.4-1
## [65] caTools_1.18.2 knitr_1.28