Commit 090c3745 authored by Holger Brandl's avatar Holger Brandl
Browse files

gene loss count debugging for liver proteomics

parent 8b8232f6
......@@ -22,14 +22,15 @@ opts = docopt(doc, commandArgs(TRUE))
## load some packages first because of name-space order
library(knitr)
#source("https://bioconductor.org/biocLite.R")
#biocLite("DESeq2")
library(DESeq2)
#load_pack(readr)
#devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.23/R/core_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.34/R/core_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.34/R/ggplot_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.34/R/bio/diffex_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.35/R/core_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.35/R/ggplot_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.35/R/bio/diffex_commons.R")
## process input arguments
......@@ -168,10 +169,11 @@ if(!all(contrastSamples %in% designSamples)){
dds = DESeqDataSetFromMatrix(countMatrix, exDesign, formula(as.formula(paste("~", designFormula))))
#dds = estimateSizeFactors(dds)
#dds = estimateDispersions(dds)
dds = DESeq(dds)
dds = DESeq(dds, parallel=T)
session::save.session("fc_debug.dat");
# session::restore.session("fc_debug.dat");
install_package("session")
session::save.session(".fc_debug.dat");
# session::restore.session(".fc_debug.dat");
########################################################################################################################
......@@ -201,6 +203,18 @@ sizeFactors(dds) %>%
#' The estimation of dispersions 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.
plotDispEsts(dds, main="Dispersion plot")
## DEBUG-START
if(F){
## batch factor variable importance
ddsred <- DESeq(dds, test="LRT", reduced=~sequence_start_date, parallel=T)
summary(batchImpact)
results(ddsred) %>% as_df %>% count(padj< 0.01)
}
## DEBUG-END
## todo add comparison to non-corrected model here using LTR as detailed out in bioinfo/rnaseq_diffex.md:16
......@@ -217,6 +231,7 @@ plotDispEsts(dds, main="Dispersion plot")
# Regularized log transformation for clustering/heatmaps, etc
rld = rlogTransformation(dds)
#varstab = vst(dds) ## vst is supposed to be much 1k x faster
#plotPCA(rld, intgroup = "sample")
plotPCA(rld, intgroup = contrastAttribute)
......@@ -252,6 +267,28 @@ distMatrix %>% d3heatmap(xaxis_height=1)
## see http://www.bioconductor.org/help/workflows/rnaseqGene/
## DEBUG-START
if(F){
load_pack("genefilter")
topVarGenes <- head(order(rowVars(assay(rld)),decreasing=TRUE),20)
mat <- assay(rld)[ topVarGenes, ]
mat <- mat - rowMeans(mat)
df <- as.data.frame(colData(rld)[,c("cell","dex")])
#pheatmap(mat, annotation_col=df)
pheatmap(mat, annotation_col=exDesign)
}
## DEBUG-END
## DEBUG-START
if(F){
## use pca explorer to explore valiadate
source("https://bioconductor.org/biocLite.R")
biocLite("pcaExplorer")
}
## DEBUG-END
########################################################################################################################
#' # Perform Differential Expression Analysis
......@@ -534,7 +571,7 @@ degs = deAnnot %>% filter(is_hit)
#' For how many genes did Deseq did not assign a p-value (due to zero counts or outliers; see section 1.5.3 os [DEseq manual](https://bioconductor.org/packages/release/bioc/vignettes/DESeq/inst/doc/DESeq.pdf))
degs %>% count(sample_1, sample_2, is.na(pvalue)) %>% mutate(contrast=paste(sample_1, "vs", sample_2)) %>%
ggplot(aes(contrast, n)) + geom_bar(stat="identity") + ggtitle("Genes with NA-palue")
ggplot(aes(contrast, n)) + geom_bar(stat="identity") + ggtitle("Genes with NA-palue") + coord_flip()
# degs %>%
......
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