diff --git a/dge_workflow/featcounts_deseq.R b/dge_workflow/featcounts_deseq.R index 0c81e14117d5a6dd323a8587ab30a49e06506d30..72baaeccc6b6fd9b3f4c6c7beb303ab1160e567c 100755 --- a/dge_workflow/featcounts_deseq.R +++ b/dge_workflow/featcounts_deseq.R @@ -19,14 +19,15 @@ opts <- docopt(doc, commandArgs(TRUE)) #opts <- docopt(doc, "--contrasts ~/MPI-CBG_work/P5_DESeq/dba_contrasts_human.txt ~/MPI-CBG_work/P5_DESeq/countMatrix_human.txt") ## opts <- docopt(doc, "countMatrix.txt") +## opts <- docopt(doc, "star_counts_matrix.txt") require(knitr) require(DESeq2) -devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.12/R/core_commons.R") -devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.12/R/ggplot_commons.R") -devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.12/R/bio/diffex_commons.R") +devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.14/R/core_commons.R") +devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.14/R/ggplot_commons.R") +devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.14/R/bio/diffex_commons.R") require.auto(DT) @@ -149,7 +150,7 @@ summary(results(dds)) #' The dispersion plot shows how the estimates are shrunk from the gene wise values (black dots) toward the fitted estimates, with the final values used in testing being the blue dots. #' The dispersion can be understood as the square of the coefficient of biological variation. So, if a gene's expression typically differs from replicate to replicate sample by 20%, this gene's dispersion is: .20^2 = .04. ## The function estimateDispersions performs three steps. First, it estimates, for each gene and each (replicated) condition, a dispersion value, then, it fits, for each condition, a curve through the estimates. Finally, it assigns to each gene a dispersion value, using either the estimated or the fitted value. -dds <- estimateDispersions(dds) +#dds <- estimateDispersions(dds) plotDispEsts(dds, main="Dispersion plot") ######################################################################################################################## @@ -192,7 +193,6 @@ heatmap.2(distMatrix, labRow=colnames(labelcntData), labCol=colnames(labelcntDat ## extract all de-sets according to our contrasts model deResults <- adply(contrasts, 1, splat(function(sample_1, sample_2){ - # browser() results(dds, contrast=c("condition", sample_1, sample_2)) %>% rownames2column("ensembl_gene_id") %>% as.data.frame() %>% @@ -263,8 +263,8 @@ baseMeanPerLvl <- sapply( levels(dds$condition), function(lvl) rowMeans( counts( ## add base means to diffßex summary deResults <- baseMeanPerLvl %>% - gather(sample, norm_count, -ensembl_gene_id) %>% - merge(.,., by="ensembl_gene_id", suffixes=c("_1", "_2")) %>% + gather(sample, norm_count, -ensembl_gene_id) %>% + merge(.,., by="ensembl_gene_id", suffixes=c("_1", "_2")) %>% # filter(ac(sample_1)<ac(sample_2)) %>% # add diffex status merge(deResults)