ART_methylation/run_medips.R at master · markziemann/ART_methylation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
# this part is run on a big server with lots of disk/memory

library("BSgenome")
library("MEDIPS")
library(BSgenome.Hsapiens.UCSC.hg19)
BSgenome="BSgenome.Hsapiens.UCSC.hg19"
uniq=1e-3
extend=300
shift=0
ws=100
# modify this
#chr.select="chr22"

# list all files
bams <- list.files(".",pattern="bam$")

# read em all in
xx <- sapply(bams,function(x) {
  MEDIPS.createSet(file = x,
  BSgenome = BSgenome, extend = extend, shift = shift, uniq = uniq,
  window_size = ws)
})

names(xx) <- gsub(".bam","",names(xx))

# need to read in the sample groups based on "medical help to conceive"
s <- read.table("samplegroups.txt")

# separate two groups
mset1 <- xx[which(names(xx) %in% s[which(s$V2=="no"),1] )]
mset2 <- xx[which(names(xx) %in% s[which(s$V2=="yes"),1] )]

## CBMC first - select "C" samples first
mset1 <- mset1[grep("C",names(mset1))]
mset2 <- mset2[grep("C",names(mset2))]

# prepare the coupling set - not sure what that is but what the hey
CS <- MEDIPS.couplingVector(pattern = "CG", refObj = xx[[1]])

# run the contrast
mr.edgeR <- MEDIPS.meth(
  MSet1 = mset1,
  MSet2 = mset2,
  CSet = CS,
  p.adj = "fdr",
  diff.method = "edgeR",
  MeDIP = TRUE,
  CNV = FALSE,
  minRowSum = 10)

mr.edgeR.c <- mr.edgeR[order(mr.edgeR$edgeR.p.value),]

write.table(mr.edgeR.c, file="mr.edgeR.c.tsv",sep="\t",quote=FALSE)

# WBC
# separate two groups
mset1 <- xx[which(names(xx) %in% s[which(s$V2=="no"),1] )]
mset2 <- xx[which(names(xx) %in% s[which(s$V2=="yes"),1] )]

## CBMC first - select "C" samples first
mset1 <- mset1[grep("W",names(mset1))]
mset2 <- mset2[grep("W",names(mset2))]

# prepare the coupling set - not sure what that is but what the hey
CS <- MEDIPS.couplingVector(pattern = "CG", refObj = xx[[1]])

# run the contrast
mr.edgeR <- MEDIPS.meth(
  MSet1 = mset1,
  MSet2 = mset2,
  CSet = CS,
  p.adj = "fdr",
  diff.method = "edgeR",
  MeDIP = TRUE,
  CNV = FALSE,
  minRowSum = 10)

mr.edgeR.w <- mr.edgeR[order(mr.edgeR$edgeR.p.value),]

write.table(mr.edgeR.w, file="mr.edgeR.w.tsv",sep="\t",quote=FALSE)

save.image("castillo.Rdata")