Commit 5813a3b1 authored by Safia Saci's avatar Safia Saci
Browse files

Upload New File

parent 7ccc5b23
library("dmrseq")
library("bsseq")
library("rtracklayer")
library("BiocParallel")
library("purrr")
library("stringr")
library("tibble")
register(MulticoreParam(2))
windows_genome <- readRDS("windows_genome.RDS")
get_cov_matrix <- function(files){
cov <- lapply(files, readRDS)
do.call(cbind, cov)
}
cov <- get_cov_matrix(list.files("/work2/genphyse/genepi/guillaume/rosepigs/rds_file/", pattern = "^G11.*" , full.names = TRUE))
dim(cov)
files <- list.files("/work2/genphyse/genepi/guillaume/rosepigs/rds_file", pattern = "^G11.*" , full.names = TRUE)
get_m_matrix <- function(files){
mvalue <- lapply(files, function(file){
import(file, format = "BigWig")$score
})
do.call(cbind, mvalue)
}
mvalue <- get_m_matrix(list.files("/work2/genphyse/genepi/guillaume/rosepigs/bigwig/", pattern = "^G11.*" , full.names = TRUE))
mvalue <- mvalue * cov
pData <- tibble(
path = files,
file = map_chr(strsplit(path, "/", fixed = TRUE), 8),
name = map_chr(strsplit(file, ".", fixed = TRUE), 1),
line = str_sub(file, end = 4),
condition = map_chr(strsplit(file, "_", fixed = TRUE), 2)
)
keep <- (rowSums(cov[ ,pData$line == "G11-"]) > 0) & (rowSums(cov[ ,pData$line == "G11+"]) > 0) &
(rowSums(cov[ ,pData$condition == "fasted"]) > 0) & (rowSums(cov[ ,pData$condition == "fed"]) > 0) &
(rowSums(cov) <= 10000)
window <- windows_genome[keep]
cov_keep <- cov[keep, ]
mvalue_keep <- mvalue[keep, ]
my_bsseq <- BSseq(M = mvalue_keep,
Cov = cov_keep,
gr = window,
sampleNames = pData$name,
rmZeroCov = FALSE)
rm(cov, mvalue, windows_genome, cov_keep, mvalue_keep, window)
pData(my_bsseq) <- pData
#bsseq_test <- chrSelectBSseq(my_bsseq, seqnames = "12")
result_line <- dmrseq(my_bsseq,
testCovariate = "line",
adjustCovariate = "condition",
minNumRegion = 3,
smooth = FALSE)
result_cond <- dmrseq(bsseq_test,
testCovariate = "condition",
adjustCovariate = "line",
minNumRegion = 3,
smooth = FALSE)
saveRDS(result_line, file = "dmrseq_result_line_.RDS")
saveRDS(result_cond, file = "dmrseq_result_cond.RDS")
png("dmrseq_line.png")
plotDMRs(my_bsseq,
extend = 8000,
regions=result_line[1,],
testCovariate="line")
dev.off()
saveRDS(test_result_cond, file = "dmrseq_result_cond.RDS")
png("dmr_test_cond.png")
plotDMRs(bsseq_test_cond,
extend = 8000,
regions=result_cond[1,],
testCovariate="condition")
dev.off()
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