Commit c1557a82 authored by Safia Saci's avatar Safia Saci
Browse files

Upload New File

parent 516a97e3
# 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)
# load a gtf file
gtf <- import(con = "data/Sus_scrofa.Sscrofa11.1.102.gtf.gz", format = "GFF")
# load a txi table
load("data/txi.RData")
# csv file with appetite regulating hormones
reacto <- read_tsv("data/Participating Molecules [R-HSA-400508].tsv")
genes <- gtf[gtf$type == "gene", ]
table(reacto$Name %in% genes$gene_name)
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" ]
# applicate student test for G11-
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
)
})
)
# applicate student test for G11+
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
)
})
)
# adjust the pvalues with with FDR method
nodes_m$adjpval <- p.adjust(nodes_m$pval, method = "fdr")
nodes_p$adjpval <- p.adjust(nodes_p$pval, method = "fdr")
Markdown is supported
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