Commit 63c8d58f authored by Guillaume Devailly's avatar Guillaume Devailly
Browse files

2 new figures

parent 38eb6585
library(dmrseq)
library(GenomicRanges)
library(UpSetR)
library(eulerr)
# parsing & filtering -----------
dmrline <- readRDS("~/mnt/genotoul_genepi/guillaume/rosepigs/dmrseq_line_gl_.RDS")
difpeak <- readRDS("~/mnt/genotoul_genepi/guillaume/rosepigs/dba_report_peaks_line.RDS")
difcons <- readRDS("~/mnt/genotoul_genepi/guillaume/rosepigs/dba_report_Cons_line.RDS")
dmrline <- dmrline[dmrline$qval <= 0.05]
# crossing -------
length(dmrline) # 514
length(difpeak) # 563
length(difcons) # 1855
sum(countOverlaps(dmrline, difpeak)) # 59
sum(countOverlaps(dmrline, difcons)) # 177
sum(countOverlaps(difcons, difpeak)) # 239
sum(countOverlaps(dmrline, difcons[overlapsAny(difcons, difpeak)])) # 36
fit1 <- euler(c(
"DMRseq" = 514 - 58 - 177 + 36,
"diffBind Consensus" = 1855 - 177 - 239 + 36,
"diffBind peaks" = 563 - 59 - 239 + 36,
"DMRseq&diffBind Consensus" = 177 - 36,
"DMRseq&diffBind peaks" = 59 - 36,
"diffBind Consensus&diffBind peaks" = 239 - 36,
"DMRseq&diffBind Consensus&diffBind peaks" = 36
))
fontsize = 32
svglite::svglite("~/work/dmr_venn.svg")
plot(
fit1,
fills = FALSE,
quantities = list(fontsize = fontsize),
labels = list(fontsize = fontsize),
main = list(label = "Line", fontsize = fontsize),
edges = list(col = "dodgerblue", lwd = 3)
)
dev.off()
# setwd("/home/gdevailly/work/project/rosepigs/rosepigs")
library(data.table)
library(rtracklayer)
library(tximport)
library(stringr)
library(ggplot2)
library(readr)
library(dplyr)
library(visNetwork)
library(purrr)
# metadata --------
metadata <- fread("data/20201104_Genealogie_24anx_ROSEpigs.txt")
# nf_core metadata had the wrong conditions
# we correct our mistake
old_metadata <- fread("data/ROSEpigs_rnaseq_file.csv")
old_metadata$nanim <- str_extract(old_metadata$fastq_1, "FR17MAG[:digit:]*")
old_metadata$nfcore_id <- with(old_metadata, paste0(group, "_R", replicate))
old_metadata <- unique(old_metadata[, c("nanim", "nfcore_id")])
metadata <- merge(x = metadata, y = old_metadata, by.x = "n_nation", by.y = "nanim")
metadata$id <- with(metadata, paste(
line,
condition,
sex,
str_trunc(n_nation, 4, "left", ellipsis = ""),
sep = "_")
) # /!\ nfcore_id has the wrong conditions!
rm(old_metadata)
metadata$family <- rep(LETTERS[1:6], each = 4)
gtf <- import(con = "data/Sus_scrofa.Sscrofa11.1.102.gtf.gz", format = "GFF")
tx2gene <- gtf[gtf$type == "transcript", c("transcript_id", "gene_id")]
tx2gene <- as.data.frame(tx2gene)[ , c("transcript_id", "gene_id")]
files <- file.path("data", "salmon" , metadata$nfcore_id, "quant.sf")
txi <- tximport(files, type = "salmon", tx2gene = tx2gene, countsFromAbundance = "lengthScaledTPM") # browseVignettes("tximport")
# network ------------
reacto <- read_tsv("~/work/project/rosepigs/Participating Molecules [R-HSA-400508].tsv")
genes <- gtf[gtf$type == "gene", ]
table(reacto$Name %in% genes$gene_name)
# FALSE TRUE
# 1 22
genes[genes$gene_id == "genes", ]
nodes <- left_join(reacto, as.data.frame(mcols(genes)[, c("gene_id", "gene_name")]), by = c("Name" = "gene_name"))
nodes_m_fas <- txi$abundance[rownames(txi$abundance) %in% nodes$gene_id, metadata$line == "G11-" & metadata$condition == "fasted"]
nodes_m_fed <- txi$abundance[rownames(txi$abundance) %in% nodes$gene_id, metadata$line == "G11-" & metadata$condition == "fed" ]
nodes_p_fas <- txi$abundance[rownames(txi$abundance) %in% nodes$gene_id, metadata$line == "G11+" & metadata$condition == "fasted"]
nodes_p_fed <- txi$abundance[rownames(txi$abundance) %in% nodes$gene_id, metadata$line == "G11+" & metadata$condition == "fed" ]
nodes_m <- cbind(
nodes,
map_dfr(nodes$gene_id, function(x) {
if(is.na(x)) return(tibble(log2fc = NA, pval = NA))
tibble(
log2fc = log2(median(nodes_m_fed[x, ], na.rm = TRUE) / median(nodes_m_fas[x, ], na.rm = TRUE)),
pval = t.test(nodes_m_fed[x, ], nodes_m_fas[x, ])$p.value
)
})
)
nodes_p <- cbind(
nodes,
map_dfr(nodes$gene_id, function(x) {
if(is.na(x)) return(tibble(log2fc = NA, pval = NA))
tibble(
log2fc = log2(median(nodes_p_fed[x, ], na.rm = TRUE) / median(nodes_p_fas[x, ], na.rm = TRUE)),
pval = t.test(nodes_p_fed[x, ], nodes_p_fas[x, ])$p.value
)
})
)
nodes_m$adjpval <- p.adjust(nodes_m$pval, method = "fdr")
nodes_p$adjpval <- p.adjust(nodes_p$pval, method = "fdr")
# https://reactome.org/PathwayBrowser/#/R-HSA-400508&SEL=R-HSA-400482&DTAB=MT
edges <- tribble(
~from, ~to,
"CTNNB1", "GCG",
"TCF7L2", "GCG",
"CTNNB1", "TCF7L2",
"CDX2", "GCG",
"PAX6", "GCG",
"SEC11A", "GCG",
"SEC11C", "GCG",
"SPCS1", "GCG",
"SPCS2", "GCG",
"SPCS3", "GCG",
"SEC11A", "SEC11C",
"SPCS1", "SPCS2",
"SPCS2", "SPCS3",
"SPCS1", "SPCS3",
"PCSK1", "GCG",
"GNB3", "GCG",
"GNG13", "GCG",
"GNAT3", "GCG",
"GNB3", "GNAT3",
"GNB3", "GNG13",
"GNG13", "GNAT3",
"FFAR4", "GCG",
"GPR119", "GCG",
"FFAR1", "GCG",
"LEP", "GCG",
"GRP", "GCG",
"DPP4", "GCG",
"PAX6", "GIP",
"ISL1", "GIP",
"GATA4", "GIP",
"SEC11A", "GIP",
"SEC11C", "GIP",
"SPCS1", "GIP",
"SPCS2", "GIP",
"SPCS3", "GIP",
"PCSK1", "GIP",
"DPP4", "GIP",
"FFAR1", "GIP",
"GPR119", "GIP"
)
# https://stackoverflow.com/a/18749392
map2color <-function(x,pal,limits=NULL){
if(is.null(limits)) limits=range(x)
pal[findInterval(x,seq(limits[1],limits[2],length.out=length(pal)+1), all.inside=TRUE)]
}
mynodes_m <- select(
nodes_m,
id = Name,
label = Name,
log2fc,
adjpval
) %>%
mutate(
shape = "ellipse",
color.border = "black",
borderWidth = case_when(
adjpval < 0.01 ~ 3,
adjpval < 0.05 ~ 2,
TRUE ~ 0
),
color.background = map2color(
if_else(is.na(log2fc), 0, log2fc),
colorRampPalette(c("red", "white", "blue"))(21),
limits = c(-2, 2)
)
) %>%
mutate(color.background = if_else(is.na(log2fc), "lightgrey", color.background))
mynodes_p <- select(
nodes_p,
id = Name,
label = Name,
log2fc,
adjpval
) %>%
mutate(
shape = "ellipse",
color.border = "black",
borderWidth = case_when(
adjpval < 0.01 ~ 3,
adjpval < 0.05 ~ 2,
TRUE ~ 0
),
color.background = map2color(
if_else(is.na(log2fc), 0, log2fc),
colorRampPalette(c("red", "white", "blue"))(21),
limits = c(-2, 2)
)
) %>%
mutate(color.background = if_else(is.na(log2fc), "lightgrey", color.background))
visNetwork(nodes = mynodes_m, edges = edges, main = "G11-") %>%
visPhysics(barnesHut = list(gravitationalConstant = -300)) %>%
visLayout(improvedLayout = TRUE, randomSeed = 1)
visNetwork(nodes = mynodes_p, edges = edges, main = "G11+") %>%
visPhysics(barnesHut = list(gravitationalConstant = -300)) %>%
visLayout(improvedLayout = TRUE, randomSeed = 1)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment