ggbio_methylation.R

library(plyr)
library(ggbio)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(rtracklayer)


#=================================================#
#=================================================#


getSeqlengths <- function(dat)
{
	a <- with(dat, tapply(BP, CHR, max))
	names(a) <- paste("chr", names(a), sep="")
	a <- a[sort(names(a))]
	return(as.numeric(a))
}


#=================================================#
#=================================================#


anno <- read.table("~/repo/methylation/probe_annotation.txt.gz", header=T, sep="\t")
probe <- read.table("~/repo/methylation/rs111482415_methylation.txt.gz", header=T)

dat <- merge(probe, subset(anno, select=c(Name, Type, SNP_Inside, SNP_Outside, Location)), by.x="PROBE", by.y="Name", all.x=TRUE)
dat <- dat[order(dat$CHR, dat$BP), ]
head(dat)
save(dat, file="~/repo/methylation/rs111482415_methylation.RData")


#=================================================#
#=================================================#


load("~/repo/methylation/rs111482415_methylation.RData")

gr <- GRanges(seqnames = paste("chr", dat$CHR, sep=""),
	IRanges(start = dat$BP, width=1),
	strand = "+",
	Type = dat$Type,
	SNP_Inside = dat$SNP_Inside,
	SNP_Outside = dat$SNP_Outside,
	Location = dat$Location,
	beta = dat$EFFECT,
	pval = dat$P,
	probe = dat$PROBE,
	significant = dat$P < 0.05 / nrow(dat)
)
seqlengths(gr) <- getSeqlengths(dat)


#=================================================#
#=================================================#


sig <- GRanges("chr20", IRanges(1, 1130100000))
a <- subsetByOverlaps(gr, sig)
seqlevels(a) <- as.character(unique(seqnames(a)))


meth <- autoplot(a, 
	geom = "point", 
	aes(y = beta, x=start, colour=significant), 
	space.skip = 0.01) + 
	scale_colour_manual(values=c("grey", "black")) +
	geom_hline(yintercept=0, linetype="dotted")
meth




data(genesymbol, package = "biovizBase")
genesymbol["HLA-H"]

genes <- subsetByOverlaps(genesymbol, sig)

plot1 <- autoplot(genes, fill="grey", stat="reduce", facets = strand ~ .) + geom_text(aes(x=start, label=symbol, size=0.5), y=1, angle=45)
plot1

tracks(Effect=meth, Genes=plot1, heights=c(4,2))




txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
isActiveSeq(txdb)[seqlevels(txdb)] <- FALSE
isActiveSeq(txdb) <- c("chr6"=TRUE)


genes <- autoplot(txdb, which = sig)
genes2 <- autoplot(txdb, which = sig, stat="reduce")
genes2
genes3 <- autoplot(txdb, which = genesymbol["HLA-H"], coord="genome", stat="reduce")
genes3


hlah <- subsetByOverlaps(txdb, genesymbol["HLA-H"])

tracks(meth, genes)

hlah <- genesymbol["HLA-H"]
seqlevels(hlah) <- as.character(unique(seqnames(hlah)))
genes4 <- autoplot(hlah)
genes4


p.ideo <- plotIdeogram(genome = "hg19")
wh <- GRanges("chr16", IRanges(30064491, 30081734))
p1 <- autoplot(txdb, which = wh, names.expr = "gene_id")
p2 <- autoplot(txdb, which = wh, stat = "reduce", color = "brown",
               fill = "brown")
tracks(p.ideo, full = p1, reduce = p2, heights = c(1.5, 5, 1)) + ylab("") + theme_tracks_sunset()