-
Notifications
You must be signed in to change notification settings - Fork 26
/
mutation_overlap.R
75 lines (47 loc) · 1.29 KB
/
mutation_overlap.R
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
rm(list = ls())
gc()
library(data.table)
library(stringr)
library(vcfR)
library(dplyr)
vcf_read_one <- read.vcfR("file1.vcf",
verbose = FALSE )
vcf_one = vcfR::getFIX(vcf_read_one) |> as.data.frame() |> setDT()
vcf_one$scenario = "One"
vcf_read_two <- read.vcfR("file2.vcf",
verbose = FALSE )
vcf_two = vcfR::getFIX(vcf_read_two) |> as.data.frame() |> setDT()
vcf_two$scenario = "Two"
rm(vcf_read_one, vcf_read_two)
x = rbind(vcf_one, vcf_two)
rm(vcf_one, vcf_two)
y = x[, c("CHROM", "POS", "REF", "ALT", "scenario"), with = FALSE]
y$mut = paste(y$CHROM, y$POS, y$REF, y$ALT, sep = ":")
z = y[, by = mut, .(
scenarios = paste(scenario, collapse = "+")
)]
z = z[order(z$scenarios), ]
fwrite(
z, "All_mutation-overlaps.csv",
row.names = FALSE, quote = TRUE, sep = ","
)
library(ggvenn)
library(ggplot2)
y = split(y, y$scenario)
y = list(
'One' = y$one$mut,
'Two' = y$two$mut
)
gr = ggvenn(
y,
fill_color = c("#43ae8d", "#ae4364")
) +
coord_equal(clip = "off")
ggsave(
plot = gr, filename = "Plots/All-overlap-plot.pdf", device = cairo_pdf,
width = 8, height = 8, units = "in"
)
ggsave(
plot = gr, filename = "Plots/All-overlap-plot.jpeg",
width = 8, height = 8, units = "in", dpi = 600
)