|
4
|
1 #!/usr/bin/Rscript --vanilla
|
|
|
2 library(getopt)
|
|
|
3 spec <- matrix(c("help", "h", 0, "logical", "view this help", "segfile1", "s", 1, "character", "seg file 1", "segfile2", "t", 1, "character", "seg file 2", "output", "o", 1, "character", "output file", "fdr", "f", 2, "character", paste("fdr method [", paste(p.adjust.methods, collapse="|"), "]", sep=""),
|
|
|
4 "reference", "r", 2, "character", "reference to use [b37|hg19|GRCh37|mm9|NCBIM37|mm10|GRCm38|dm3|BDGP5]", "annot", "a", 2, "character", "annotation to add [both|gene|cpg]", "processes", "p", 2, "integer", "number of cluster instances to open [1]"), ncol=5, byrow=T)
|
|
|
5 opt <- getopt(spec)
|
|
|
6
|
|
|
7 # set default options
|
|
|
8 if(is.null(opt$fdr)) opt$fdr <- "BH"
|
|
|
9 if(is.null(opt$reference)) opt$reference <- "hg19"
|
|
|
10 if(is.null(opt$annot)) opt$annot <- "both"
|
|
|
11 if(is.null(opt$processes)) opt$processes <- 1
|
|
|
12
|
|
|
13 # check if any invalid options
|
|
|
14 if(! opt$annot %in% c("both", "gene", "cpg") || !opt$fdr %in% p.adjust.methods || !opt$reference %in% c("b37", "hg19", "GRCh37", "mm9", "NCBIM37", "mm10", "GRCm38", "dm3", "BDGP5")) {
|
|
|
15 opt$help <- 1
|
|
|
16 }
|
|
|
17
|
|
|
18 # print help file if any incorrect
|
|
|
19 if(!is.null(opt$help) || is.null(opt$segfile1) || is.null(opt$output) || is.null(opt$output) || is.na(opt$processes)) {
|
|
|
20 self = commandArgs()[1];
|
|
|
21 cat(getopt(spec, usage=T, command="differential_methylation.R"))
|
|
|
22 q(status=1)
|
|
|
23 }
|
|
|
24
|
|
|
25 library(snow)
|
|
|
26 cl <- makeCluster(opt$p, type = "MPI")
|
|
|
27
|
|
|
28 segfile1 <- read.table(opt$segfile1, sep="\t", as.is=T, head=T, quote = "")
|
|
|
29 segfile2 <- read.table(opt$segfile2, sep="\t", as.is=T, head=T, quote = "")
|
|
|
30
|
|
|
31 rownames(segfile1) <- paste(segfile1[,2], ":", segfile1[,3], "-", segfile1[,4], sep="")
|
|
|
32 rownames(segfile2) <- paste(segfile2[,2], ":", segfile2[,3], "-", segfile2[,4], sep="")
|
|
|
33
|
|
|
34 common_reg <- intersect(rownames(segfile1), rownames(segfile2))
|
|
|
35 prop_test <- function(x, n, fdr) {
|
|
|
36 library(abind)
|
|
|
37 ESTIMATE <- x/n
|
|
|
38 DELTA <- ESTIMATE[,1L] - ESTIMATE[,2L]
|
|
|
39 YATES <- pmin(abs(DELTA)/rowSums(1/n), 0.5)
|
|
|
40 p <- rowSums(x)/rowSums(n)
|
|
|
41 df <- 1
|
|
|
42 x <- abind(x, n - x, along=3)
|
|
|
43 E <- abind(n*p, n*(1-p), along=3)
|
|
|
44 STATISTIC <- rowSums((abs(x - E) - YATES)^2/E)
|
|
|
45 STATISTIC[is.na(STATISTIC)] <- 0
|
|
|
46 PVAL <- pchisq(STATISTIC, df, lower.tail = FALSE)
|
|
|
47 return(data.frame(X_Squared = round(STATISTIC, 3), P.value= round(PVAL, 6), P.adjusted = round(p.adjust(PVAL, method=fdr), 6)))
|
|
|
48 }
|
|
|
49 results <- prop_test(cbind(segfile1[common_reg, "Methylated"], segfile2[common_reg, "Methylated"]), cbind(segfile1[common_reg, "Total"], segfile2[common_reg, "Total"]), opt$fdr)
|
|
|
50
|
|
|
51 sample1 <- segfile1[1, 1]
|
|
|
52 sample2 <- segfile2[1, 1]
|
|
|
53 tab_out <- data.frame(ID = paste(sample1, "vs", sample2, sep="."), segfile1[common_reg, c(2:4)], segfile1[common_reg, c("Methylated", "Total", "FractionMethylated")],
|
|
|
54 segfile2[common_reg, c("Methylated", "Total", "FractionMethylated")], results, DiffProp = round(segfile1[common_reg, "FractionMethylated"] - segfile2[common_reg, "FractionMethylated"], 6), stringsAsFactors=F)
|
|
|
55 colnames(tab_out)[c(5:10)] <- c(paste(sample1, "Methylated", sep="."), paste(sample1, "Total", sep="."), paste(sample1, "Proportion", sep="."), paste(sample2, "Methylated", sep="."),
|
|
|
56 paste(sample2, "Total", sep="."), paste(sample2, "Proportion", sep="."))
|
|
|
57
|
|
|
58 rm(segfile1)
|
|
|
59 rm(segfile2)
|
|
|
60 gc()
|
|
|
61
|
|
|
62 add_annot <- function(tab, annotation, genome) {
|
|
|
63 # find closest feature
|
|
|
64 if(genome == "hg19" || genome == "GRCh37" || genome == "b37") {
|
|
|
65 genome <- "hg19"
|
|
|
66 dataset <- "hsapiens_gene_ensembl"
|
|
|
67 biomart <- "ensembl"
|
|
|
68 host <- "www.biomart.org"
|
|
|
69 attributes <- c("ensembl_gene_id", "hgnc_symbol", "refseq_mrna", "start_position", "end_position")
|
|
|
70 } else if(genome == "mm9" || genome == "NCBIM37") {
|
|
|
71 genome <- "mm9"
|
|
|
72 dataset <- "mmusculus_gene_ensembl"
|
|
|
73 biomart <- "ENSEMBL_MART_ENSEMBL"
|
|
|
74 host <- "may2012.archive.ensembl.org"
|
|
|
75 attributes <- c("ensembl_gene_id", "mgi_symbol", "refseq_mrna", "start_position", "end_position")
|
|
|
76 } else if(genome == "mm10" || genome == "GRCm38") {
|
|
|
77 genome <- "mm10"
|
|
|
78 dataset <- "mmusculus_gene_ensembl"
|
|
|
79 biomart <- "ensembl"
|
|
|
80 host <- "www.biomart.org"
|
|
|
81 attributes <- c("ensembl_gene_id", "mgi_symbol", "refseq_mrna", "start_position", "end_position")
|
|
|
82 } else if(genome == "dm3" || genome == "BDGP5") {
|
|
|
83 genome <- "dm3"
|
|
|
84 dataset <- "dmelanogaster_gene_ensembl"
|
|
|
85 biomart <- "ensembl"
|
|
|
86 host <- "www.biomart.org"
|
|
|
87 attributes <- c("ensembl_gene_id", "flybasecgid_gene", "refseq_mrna", "start_position", "end_position")
|
|
|
88 }
|
|
|
89
|
|
|
90 e_to_U = c("GL000191.1" = "1_gl000191_random", "GL000192.1" = "1_gl000192_random", "GL000193.1" = "4_gl000193_random", "GL000194.1" = "4_gl000194_random",
|
|
|
91 "GL000195.1" = "7_gl000195_random", "GL000196.1" = "8_gl000196_random", "GL000197.1" = "8_gl000197_random", "GL000198.1" = "9_gl000198_random",
|
|
|
92 "GL000199.1" = "9_gl000199_random", "GL000200.1" = "9_gl000200_random", "GL000201.1" = "9_gl000201_random", "GL000202.1" = "11_gl000202_random",
|
|
|
93 "GL000203.1" = "17_gl000203_random", "GL000204.1" = "17_gl000204_random", "GL000205.1" = "17_gl000205_random", "GL000206.1" = "17_gl000206_random",
|
|
|
94 "GL000207.1" = "18_gl000207_random", "GL000208.1" = "19_gl000208_random", "GL000209.1" = "19_gl000209_random", "GL000210.1" = "21_gl000210_random",
|
|
|
95 "GL000211.1" = "Un_gl000211", "GL000212.1" = "Un_gl000212", "GL000213.1" = "Un_gl000213", "GL000214.1" = "Un_gl000214",
|
|
|
96 "GL000215.1" = "Un_gl000215", "GL000216.1" = "Un_gl000216", "GL000217.1" = "Un_gl000217", "GL000218.1" = "Un_gl000218",
|
|
|
97 "GL000219.1" = "Un_gl000219", "GL000220.1" = "Un_gl000220", "GL000221.1" = "Un_gl000221", "GL000222.1" = "Un_gl000222",
|
|
|
98 "GL000223.1" = "Un_gl000223", "GL000224.1" = "Un_gl000224", "GL000225.1" = "Un_gl000225", "GL000226.1" = "Un_gl000226",
|
|
|
99 "GL000227.1" = "Un_gl000227", "GL000228.1" = "Un_gl000228", "GL000229.1" = "Un_gl000229", "GL000230.1" = "Un_gl000230",
|
|
|
100 "GL000231.1" = "Un_gl000231", "GL000232.1" = "Un_gl000232", "GL000233.1" = "Un_gl000233", "GL000234.1" = "Un_gl000234",
|
|
|
101 "GL000235.1" = "Un_gl000235", "GL000236.1" = "Un_gl000236", "GL000237.1" = "Un_gl000237", "GL000238.1" = "Un_gl000238",
|
|
|
102 "GL000239.1" = "Un_gl000239", "GL000240.1" = "Un_gl000240", "GL000241.1" = "Un_gl000241", "GL000242.1" = "Un_gl000242",
|
|
|
103 "GL000243.1" = "Un_gl000243", "GL000244.1" = "Un_gl000244", "GL000245.1" = "Un_gl000245", "GL000246.1" = "Un_gl000246",
|
|
|
104 "GL000247.1" = "Un_gl000247", "GL000248.1" = "Un_gl000248", "GL000249.1" = "Un_gl000249")
|
|
|
105
|
|
|
106 # function to conver ucsc chroms to ensembl
|
|
|
107 ensembl <- function(chr, genome) {
|
|
|
108 # ensembl does not us chr and M is MT
|
|
|
109 chr <- sub("^chr", "", chr)
|
|
|
110 if(genome == "dm3" || genome == "BDGP5") {
|
|
|
111 chr <- sub("^M$", "dmel_mitochondrion_genome", chr)
|
|
|
112 } else {
|
|
|
113 chr <- sub("^M$", "MT", chr)
|
|
|
114 }
|
|
|
115 chr[which(chr %in% e_to_U)] <- names(e_to_U)[match(chr[which(chr %in% e_to_U)], e_to_U)]
|
|
|
116 return(chr)
|
|
|
117 }
|
|
|
118
|
|
|
119 # function to conver ensembl chroms to ucsc
|
|
|
120 ucsc <- function(chr) {
|
|
|
121 chr <- sub("^chr", "", chr)
|
|
|
122 tmp <- which(chr %in% names(c(e_to_U, dmel_mitochondrion_genome = "M", MT = "M")))
|
|
|
123 chr[tmp] <- c(e_to_U, dmel_mitochondrion_genome = "M", MT = "M")[chr[tmp]]
|
|
|
124 paste("chr", chr, sep="")
|
|
|
125 }
|
|
|
126
|
|
|
127 # function to get the gene annotation by chromosome
|
|
|
128 get_genes <- function(chr) {
|
|
|
129 library(biomaRt)
|
|
|
130 mart <- useMart(biomart=biomart, dataset=dataset, host=host)
|
|
|
131 tab <- getBM(attributes = attributes, filters = "chromosome_name", values = chr, mart = mart)
|
|
|
132 if(any(is.na(tab))) for(i in 1:ncol(tab)) tab[is.na(tab[,i]), i] <- ""
|
|
|
133 mult_ids <- names(which(table(tab$ensembl_gene_id) > 1))
|
|
|
134 rem <- NULL
|
|
|
135 for(i in mult_ids) {
|
|
|
136 index <- which(tab$ensembl_gene_id == i)
|
|
|
137 refseq <- which(tab$ensembl_gene_id == i & tab$refseq_mrna != "")
|
|
|
138 if(length(refseq) > 0) {
|
|
|
139 rem <- c(rem, setdiff(index, refseq[1]))
|
|
|
140 } else {
|
|
|
141 rem <- c(rem, index[-1])
|
|
|
142 }
|
|
|
143 }
|
|
|
144 if(length(rem) > 0) tab <- tab[-rem,]
|
|
|
145 colnames(tab)[4:5] <- c("feature_start", "feature_end")
|
|
|
146 return(tab)
|
|
|
147 }
|
|
|
148
|
|
|
149 get_cpg <- function(chr) {
|
|
|
150 library(rtracklayer)
|
|
|
151 options(stringsAsFactors=F)
|
|
|
152 session <- browserSession()
|
|
|
153 genome(session) <- genome
|
|
|
154 query <- ucscTableQuery(session, "CpG Islands", GRangesForUCSCGenome(genome, chr))
|
|
|
155 cpg <- getTable(query)
|
|
|
156 cpg <- cpg[, c("name", "perCpg", "perGc", "chromStart", "chromEnd")]
|
|
|
157 colnames(cpg) <- c("cpg", "cpg_perCpg", "cpg_perGc", "cpg_start", "cpg_end")
|
|
|
158 return(cpg)
|
|
|
159 }
|
|
|
160 chr <- tab[1, 2]
|
|
|
161
|
|
|
162 # get the gene info for this chrom
|
|
|
163 if(annotation == "gene") {
|
|
|
164 tab <- cbind(tab, ensembl_gene_id = "", id = "", refseq_mrna = "", feature_start = "", feature_end = "", Distance_To_Feature = "", stringsAsFactors=F)
|
|
|
165 annotData <- get_genes(ensembl(chr, genome))
|
|
|
166 colnames(tab)[16] <- attributes[2]
|
|
|
167 } else {
|
|
|
168 tab <- cbind(tab, cpg = "", cpg_perCpg = "", cpg_perGc = "", cpg_start = "", cpg_end = "", Distance_To_cpg = "", stringsAsFactors=F)
|
|
|
169 annotData <- get_cpg(ucsc(chr))
|
|
|
170 }
|
|
|
171 if(nrow(annotData) == 0) return(tab)
|
|
|
172
|
|
|
173 starts <- t(matrix(cbind(-1, as.numeric(annotData[, 4])), ncol=2, byrow=F))
|
|
|
174 ends <- t(matrix(cbind(-1, as.numeric(annotData[, 5])), ncol=2, byrow=F))
|
|
|
175
|
|
|
176 # calculate the distances from the features to the regions
|
|
|
177 dist_start_start <- matrix(cbind(tab$loc.start, 1), ncol=2, byrow=F) %*% starts
|
|
|
178 dist_start_end <- matrix(cbind(tab$loc.start, 1), ncol=2, byrow=F) %*% ends
|
|
|
179 dist_end_start <- matrix(cbind(tab$loc.end, 1), ncol=2, byrow=F) %*% starts
|
|
|
180 dist_end_end <- matrix(cbind(tab$loc.end, 1), ncol=2, byrow=F) %*% ends
|
|
|
181
|
|
|
182 # determine which regions overlap at least 1 feature
|
|
|
183 sum_signs <- abs(sign(dist_start_start) + sign(dist_start_end) + sign(dist_end_start) + sign(dist_end_end))
|
|
|
184 regions <- which(sum_signs != 4, arr.ind=TRUE)
|
|
|
185 if(length(regions) > 0) {
|
|
|
186 overlap <- sort(unique(regions[,1]))
|
|
|
187 non_overlap <- c(1:nrow(tab))[-overlap]
|
|
|
188 } else {
|
|
|
189 overlap <- NULL
|
|
|
190 non_overlap <- c(1:nrow(tab))
|
|
|
191 }
|
|
|
192
|
|
|
193 # reduce to regions with no overlaping feqature
|
|
|
194 if(length(overlap) > 0) {
|
|
|
195 dist_start_start <- matrix(dist_start_start[non_overlap,], ncol=ncol(dist_start_start))
|
|
|
196 dist_start_end <- matrix(dist_start_end[non_overlap,], ncol=ncol(dist_start_end))
|
|
|
197 dist_end_start <- matrix(dist_end_start[non_overlap,], ncol=ncol(dist_end_start))
|
|
|
198 dist_end_end <- matrix(dist_end_end[non_overlap,], ncol=ncol(dist_end_end))
|
|
|
199 }
|
|
|
200
|
|
|
201 rm(sum_signs)
|
|
|
202 gc()
|
|
|
203
|
|
|
204 # extract the annot for the regions with overlaping features
|
|
|
205 if(length(overlap) > 0) {
|
|
|
206 annot <- sapply(overlap, function(x, y) {
|
|
|
207 x <- regions[which(regions[,1] == x), 2]
|
|
|
208 sub("^//$", "", c(paste(annotData[x,1], collapse="//"), paste(annotData[x,2], collapse="//"), paste(annotData[x,3], collapse="//"), paste(annotData[x,4], collapse="//"), paste(annotData[x,5], collapse="//")))
|
|
|
209 }, regions)
|
|
|
210 annot <- as.data.frame(t(annot), stringsAsFactors=F)
|
|
|
211 annot <- cbind(annot, "Overlap", stringsAsFactors=F)
|
|
|
212 colnames(annot) <- c(colnames(annotData), tail(colnames(tab),1))
|
|
|
213 }
|
|
|
214 rm(regions)
|
|
|
215 gc()
|
|
|
216
|
|
|
217 # for non the non-overlaps the distance of the closest features to each region
|
|
|
218 if(length(non_overlap) > 0) {
|
|
|
219 clst_pts <- matrix(0, ncol=4, nrow=length(non_overlap))
|
|
|
220 clst_pts[,1] <- max.col(-abs(dist_start_start), "last")
|
|
|
221 clst_pts[,2] <- max.col(-abs(dist_start_end), "last")
|
|
|
222 clst_pts[,3] <- max.col(-abs(dist_end_start), "last")
|
|
|
223 clst_pts[,4] <- max.col(-abs(dist_end_end), "last")
|
|
|
224 dist <- matrix(0, ncol=4, nrow=length(non_overlap))
|
|
|
225 if(length(clst_pts[,1]) == 1) {
|
|
|
226 dist[,1] <- dist_start_start[1, clst_pts[,1]]
|
|
|
227 dist[,2] <- dist_start_end[1, clst_pts[,2]]
|
|
|
228 dist[,3] <- dist_end_start[1, clst_pts[,3]]
|
|
|
229 dist[,4] <- dist_end_end[1, clst_pts[,4]]
|
|
|
230
|
|
|
231 # extract the annot for the regions with non-overlaping features
|
|
|
232 clst_all <- max.col(-abs(dist))
|
|
|
233 dist_to_feat <- dist[1, clst_all]
|
|
|
234 clst <- clst_pts[1, clst_all]
|
|
|
235 } else {
|
|
|
236 dist[,1] <- dist_start_start[cbind(seq(clst_pts[,1]), clst_pts[,1])]
|
|
|
237 dist[,2] <- dist_start_end[cbind(seq(clst_pts[,2]), clst_pts[,2])]
|
|
|
238 dist[,3] <- dist_end_start[cbind(seq(clst_pts[,3]), clst_pts[,3])]
|
|
|
239 dist[,4] <- dist_end_end[cbind(seq(clst_pts[,4]), clst_pts[,4])]
|
|
|
240
|
|
|
241 # extract the annot for the regions with non-overlaping features
|
|
|
242 clst_all <- max.col(-abs(dist))
|
|
|
243 dist_to_feat <- dist[cbind(seq(clst_all), clst_all)]
|
|
|
244 clst <- clst_pts[cbind(seq(clst_all), clst_all)]
|
|
|
245 }
|
|
|
246 annot_non_overlap <- cbind(annotData[clst, ], dist_to_feat)
|
|
|
247 colnames(annot_non_overlap) <- c(colnames(annotData), tail(colnames(tab),1))
|
|
|
248 }
|
|
|
249
|
|
|
250 rm(dist_start_start)
|
|
|
251 rm(dist_start_end)
|
|
|
252 rm(dist_end_start)
|
|
|
253 rm(dist_end_end)
|
|
|
254 gc()
|
|
|
255
|
|
|
256 if(length(overlap) > 0) tab[overlap, colnames(annot)] <- annot
|
|
|
257 if(length(non_overlap) > 0) tab[non_overlap, colnames(annot_non_overlap)] <- annot_non_overlap
|
|
|
258 return(tab)
|
|
|
259 }
|
|
|
260
|
|
|
261 tab_list <- list()
|
|
|
262 for(i in unique(tab_out[,2])) tab_list[[i]] <- tab_out[which(tab_out$chrom == i),]
|
|
|
263 rm(tab_out)
|
|
|
264 gc()
|
|
|
265 if(opt$annot == "both" || opt$annot == "gene") {
|
|
|
266 annotation <- "gene"
|
|
|
267 tab_list <- clusterApplyLB(cl, tab_list, add_annot, annotation, opt$reference)
|
|
|
268 }
|
|
|
269
|
|
|
270 if(opt$annot == "both" || opt$annot == "cpg") {
|
|
|
271 annotation <- "cpg"
|
|
|
272 tab_list <- clusterApplyLB(cl, tab_list, add_annot, annotation, opt$reference)
|
|
|
273 }
|
|
|
274
|
|
|
275 tab_out <- NULL
|
|
|
276 for(i in 1:length(tab_list)) tab_out <- rbind(tab_out, tab_list[[i]])
|
|
|
277
|
|
|
278 cat("'", file=opt$output)
|
|
|
279 write.table(tab_out[,c(1:13, 15:26, 14)], opt$output, row.names=F, col.names=T, quote=F, sep="\t", append=T)
|
|
|
280 stopCluster(cl)
|
|
|
281 q(status=0)
|
