Mercurial > repos > fcaramia > methylation_analysis_bismark

comparison methylation_analysis/differential_methylation.R @ 4:282edadee017 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author fcaramia
date Mon, 03 Dec 2012 18:26:25 -0500
parents
children
comparison
equal deleted inserted replaced
3:d0fc260f6dc9 4:282edadee017
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)