Commit 7ea79694 authored by Lena Hersemann's avatar Lena Hersemann

added DESeq2 option betaprior with default = TRUE; corrected volcano plot...

added DESeq2 option betaprior with default = TRUE; corrected volcano plot labelling; removed default of geneInfo option
parent f852f099
......@@ -17,7 +17,8 @@ Options:
--design <formula> Design fomula for DeSeq with contrast attribute at the end [default: condition]
--lfc <lfc_cutoff> Just report genes with abs(lfc) > lfc_cutoff as hits [default: 1.0]
--ensembl_db <ensembl_db> Ensebmbl db to be used. If not specified inferred from data. TODO implement NONE here!!
--geneInfo <info_file> Gene information file downloaded from the ensembl website if bioMart version is not present [default: bioMart version will be used]
--geneInfo <info_file> Gene information file downloaded from the ensembl website if bioMart version is not present
--betaprior <boolean> Apply betaprior when running DESeq2 [default: TRUE]
'
#commandArgs <- function(x) c("--design", "'batch+condition'", "--contrasts", "../contrasts.txt star_counts_matrix_reduced_sampleset.txt", "basic_design_reduced_sampleset.txt")
......@@ -46,6 +47,8 @@ devtools::source_url("https://git.mpi-cbg.de/bioinfo/ngs_tools/raw/v4/dge_workfl
count_matrix_file = opts$count_matrix
design_matrix_file = opts$design_matrix
contrasts_file = opts$contrasts
use_betaPrior = as.logical(opts$betaprior)
gene_info_file = opts$geneInfo
resultsBase = if (str_length(opts$out) > 0) paste0(opts$out, ".") else ""
......@@ -203,7 +206,9 @@ if (! all(contrastSamples %in% designSamples)) {
dds = DESeqDataSetFromMatrix(countMatrix, exDesign, formula(as.formula(paste("~", designFormula))))
#dds = estimateSizeFactors(dds)
#dds = estimateDispersions(dds)
dds = DESeq(dds, parallel = T)
# dds = DESeq(dds, parallel = T)
dds = DESeq(dds, parallel = T, betaPrior = use_betaPrior)
......@@ -589,7 +594,7 @@ hitCounts = filter(deResults, is_hit) %>%
rename(hits = n) %>%
## complement hit counts with directed contrasts without any diffex genes
right_join(merge(mutate(contrasts),data_frame(c1_overex=c(T,F)), by=NULL)) %>% mutate(hits=if_else(is.na(hits), as.integer(0), hits)) %>%
merge(data.frame(c1_overex = c(T, F), x_pos = c(- maxX + 2, maxX - 2)))
merge(data.frame(c1_overex = c(T, F), x_pos = c(maxX + 2, -maxX - 2)))
hitCounts %<>% mutate(y_pos = quantile(- log10(deResults$pvalue), 0.99, na.rm = TRUE))
......@@ -641,7 +646,7 @@ ensembl_dataset= if (!is.null(opts$ensembl_db)) paste0(opts$ensembl_db, "_gene_e
# Load transcriptome annotations needed for results annotation
if (is.null(opts$geneInfo)){
if (is.null(gene_info_file)){
geneInfo = quote({
## mart = biomaRt::useDataset("drerio_gene_ensembl", mart = biomaRt::useMart("ensembl"))??
# mart = biomaRt::useDataset(guess_mart(countData$ensembl_gene_id), mart = biomaRt::useMart("ensembl"))
......@@ -656,7 +661,7 @@ if (is.null(opts$geneInfo)){
tbl_df
}) %>% cache_it("geneInfo")
} else {
geneInfo = read_tsv(opts$geneInfo)
geneInfo = read_tsv(gene_info_file)
colnames(geneInfo)[1] = "ensembl_gene_id"
}
......@@ -987,6 +992,6 @@ pheatmap(mat, annotation_col = column2rownames(exDesign, "replicate") %>% select
#' Further informtion on [FPKM and TPM](http://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/)
session::save.session(".featcounts_deseq_mf.dat");
session::save.session(".featcounts_deseq_mf.dat");
# session::restore.session(".featcounts_deseq_mf.dat");
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