mutDis.plot <- function(plot.data, sample = "sample", chr = NULL, color = NULL,
min.mut = 6, max.dis = 1000) {
if (is.null(color)) {
col = RColorBrewer::brewer.pal(n = 6, name = "Set1")
} else {
col = color
}
names(col) = c("C>A", "C>G", "C>T", "T>A", "T>C", "T>G")
build = plot.data$build[1]
genome.opts = c("hg19", "hg18", "hg38")
if (!build %in% genome.opts) {
stop("Available reference builds: hg18, hg19, hg38")
}
if (build == "hg19") {
chr.lens = c(249250621, 243199373, 198022430, 191154276, 180915260,
171115067, 159138663, 146364022, 141213431, 135534747,
135006516, 133851895, 115169878, 107349540, 102531392,
90354753, 81195210, 78077248, 59128983, 63025520, 48129895,
51304566, 155270560, 59373566)
chr.arm = c(1.25e+08, 93300000, 9.1e+07, 50400000, 48400000,
6.1e+07, 59900000, 45600000, 4.9e+07, 40200000, 53700000,
35800000, 17900000, 17600000, 1.9e+07, 36600000, 2.4e+07,
17200000, 26500000, 27500000, 13200000, 14700000, 60600000,
12500000)
} else if (build == "hg18") {
chr.lens = c(247249719, 242951149, 199501827, 191273063, 180857866,
170899992, 158821424, 146274826, 140273252, 135374737,
134452384, 132349534, 114142980, 106368585, 100338915,
88827254, 78774742, 76117153, 63811651, 62435964, 46944323,
49691432, 154913754, 57772954)
chr.arm = c(124300000, 93300000, 91700000, 50700000, 47700000,
60500000, 59100000, 45200000, 51800000, 40300000, 52900000,
35400000, 1.6e+07, 15600000, 1.7e+07, 38200000, 22200000,
16100000, 28500000, 27100000, 12300000, 11800000, 59500000,
11300000)
} else if (build == "hg38") {
chr.lens = c(248956422, 242193529, 198295559, 190214555, 181538259,
170805979, 159345973, 145138636, 138394717, 133797422,
135086622, 133275309, 114364328, 107043718, 101991189,
90338345, 83257441, 80373285, 58617616, 64444167, 46709983,
50818468, 156040895, 57227415)
chr.arm = c(123400000, 93900000, 90900000, 5e+07, 48800000,
59800000, 60100000, 45200000, 4.3e+07, 39800000, 53400000,
35500000, 17700000, 17200000, 1.9e+07, 36800000, 25100000,
18500000, 26200000, 28100000, 1.2e+07, 1.5e+07, 6.1e+07,
10400000)
} else {
stop("Available reference builds: hg18, hg19, hg38")
}
if (is.null(chr)) {
seq = c(1:24)
seq0 = c(1:22, "X", "Y")
} else {
plot.data$pos.updated = plot.data$pos
seq0 = gsub(pattern = "chr", replacement = "", x = chr, fixed = TRUE)
seq = gsub(pattern = "X", replacement = "23", x = seq0, fixed = TRUE)
seq = gsub(pattern = "Y", replacement = "24", x = seq, fixed = TRUE)
seq = as.numeric(seq)
plot.data = plot.data[which(plot.data$seq %in% seq), ]
for (i in 1:length(chr.lens)) {
if (!i %in% seq) {
chr.lens[i] = 0
}
}
}
xlim = c(0, sum(chr.lens))
chr.lens.sum = cumsum(chr.lens)
chr.lens.sum = c(0, chr.lens.sum)
plot.data$pos.updated = plot.data$pos + chr.lens.sum[plot.data$seq]
chr.abline = chr.lens[seq]
chr.abline.sum = cumsum(chr.abline)
chr.abline.sum = c(0, chr.abline.sum)
chr.axis = seq
for (i in 1:length(seq)) {
chr.axis[i] = (chr.abline.sum[i] + chr.abline.sum[i + 1])/2
}
par(mai = c(1, 1, 1, 1.5))
plot.data$dis.log10 = log10(plot.data$dis)
plot(plot.data$pos.updated, plot.data$dis.log10, col = col[plot.data$alteration],
pch = 16, cex = 0.5, xlim = xlim, axes = F, box(), mgp = c(2,
1, 0), main = paste("Rainfall plot for", sample, sep = " "),
xlab = "Chromosomes", ylab = "Intermutation distance (bp, log10)")
abline(h = 2, col = "red", lty = 2)
abline(v = chr.abline.sum, col = "grey")
axis(2, las = 1, lwd.tick = 0.5)
axis(1, at = chr.axis, labels = c(seq0), tick = FALSE, mgp = c(1,
0.5, 0))
legend("right", names(col), col = col, pch = 16, inset = c(-0.15,
0), bty = "n", xpd = TRUE)
arrows.point = c()
num = dim(plot.data)[1] - 5
katPoint <- matrix(nrow = num, ncol = 8)
i = 1
mutnum = 1
Cmutnum = 0
for (i in 1:num) {
if (plot.data$ref[i] %in% c("C", "G")){
Cmutnum = Cmutnum + 1
}
if (plot.data$dis[i + 1] <= max.dis) {
mutnum = mutnum + 1
} else {
if (mutnum >= min.mut) {
len = plot.data$pos[i] - plot.data$pos[i - mutnum +
1] + 1
chr.n = gsub(pattern = "chr", replacement = "", x = plot.data$chr[i],
fixed = TRUE)
chr.n = gsub(pattern = "X", replacement = "23", x = chr.n,
fixed = TRUE)
chr.n = gsub(pattern = "Y", replacement = "24", x = chr.n,
fixed = TRUE)
chr.n = as.numeric(chr.n)
if (plot.data$pos[i] <= chr.arm[chr.n]) {
arm = paste(chr.n, "p", sep = "")
} else if (plot.data$pos[i - mutnum + 1] >= chr.arm[chr.n]) {
arm = paste(chr.n, "q", sep = "")
} else {
arm = paste(chr.n, "p, ", chr.n, "q", sep = "")
}
wei = round(Cmutnum/mutnum,3)
katPoint[i, 1:8] = c(sample, plot.data$chr[i], plot.data$pos[i -
mutnum + 1], plot.data$pos[i], arm, len, mutnum, wei)
arrows.point = (plot.data$pos.updated[i] + plot.data$pos.updated[i -
mutnum + 1])/2
if (wei < 0.8) {
arrows(arrows.point, 0, arrows.point, 0.5, length = 0.1,
angle = 30, code = 2, col="grey")
} else {
arrows(arrows.point, 0, arrows.point, 0.5, length = 0.1,
angle = 30, code = 2, col="black")
}
}
mutnum = 1
Cmutnum = 0
}
}
message(paste(dim(na.omit(katPoint))[1], "potential kataegis events identified. Use katPoint() function to see detail.",
sep = " "))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.