Commit 531372ff authored by mmassaviol's avatar mmassaviol
Browse files

Update params, and add params for DE tools

parent c4c81153
deseq2_threads: 1
deseq2_threads: 4
fitType: mean
betaPrior: True
testType: LRT
deseq2_tx2gene: TRUE
deseq2_annotations: ""
......
......@@ -39,7 +39,11 @@ dds <- dds[keep,]
dds$condition <- factor(dds$condition, levels = unique(groups))
dds <- DESeq(dds)
if (snakemake@config[["deseq2_testType"]] == "LRT"){
dds <- DESeq(dds, fitType = snakemake@config[["deseq2_fitType"]], betaPrior = as.logical(snakemake@config[["deseq2_betaPrior"]]), test=snakemake@config[["deseq2_testType"]], reduced= ~ 1)
}else{
dds <- DESeq(dds, fitType = snakemake@config[["deseq2_fitType"]], betaPrior = as.logical(snakemake@config[["deseq2_betaPrior"]]), test=snakemake@config[["deseq2_testType"]])
}
res <- results(dds)
save.image(snakemake@output[["r_data"]])
......
......@@ -41,6 +41,26 @@
value: "",
label: "Annotation file (gtf or gff) : ",
},
{
name: deseq2_fitType,
type: radio,
selected: "parametric",
choices: [parametric: parametric, local: local, mean: mean],
label: "Type of fitting of dispersions to the mean intensity",
},
{
name: "deseq2_betaPrior",
type: "checkbox",
value: FALSE,
label: "betaPrior: whether or not to put a zero-mean normal prior on the non-intercept coefficients. By default, the beta prior is used only for the Wald test, but can also be specified for the likelihood ratio test.",
},
{
name: deseq2_testType,
type: radio,
selected: "LRT",
choices: [LRT: LRT, Wald: Wald],
label: "Test type: Use either Wald significance tests, or the likelihood ratio test on the difference in deviance between a full and reduced model formula",
},
],
},
],
......
edger_threads: 1
edger_threads: 4
edger_tx2gene: TRUE
edger_annotations: ""
normfact: TMM
Dispersion: 0
testType: exactTest
edger_outputs:
rdata: edger/data.RData
\ No newline at end of file
......@@ -7,12 +7,9 @@ library(rhdf5)
library(GenomicFeatures)
samples = snakemake@config[["samples"]]
#samples = c("DRR016125","DRR016126","DRR016127","DRR016128","DRR016129","DRR016130","DRR016131","DRR016132")
#files = file.path("/home/mbb/Documents/results2/kallisto/quant",samples,"abundance.h5")
files = snakemake@input[["counts"]]
names(files) = samples
groups = snakemake@config[["groups"]]
#groups = c("a","a","a","a","b","b","b","b")
if (snakemake@config[["edger_tx2gene"]] == TRUE){
TxDb <- makeTxDbFromGFF(file = snakemake@config[["edger_annotations"]])
......@@ -32,29 +29,40 @@ countCheck <- countsPerMillion > 1
keep <- which(rowSums(countCheck) >= 2)
dgelist <- dgelist[keep,]
dgelist <- calcNormFactors(dgelist, method="TMM")
dgelist <- calcNormFactors(dgelist, method=snakemake@config[["edger_normfact"]])
#plotMDS(dgelist)
#'grep' is a string matching function.
designMat <- model.matrix(~groups)
colnames(designMat)[2] = "groups"
de.com <- c()
dgelist$samples$group = designMat[,2]
dgelist <- estimateGLMCommonDisp(dgelist, design=designMat)
dgelist <- estimateGLMTrendedDisp(dgelist, design=designMat)
dgelist <- estimateGLMTagwiseDisp(dgelist, design=designMat)
#plotBCV(dgelist)
if (snakemake@config[["edger_testType"]] == "exactTest"){
if (snakemake@config[["edger_dispersion"]] == 0){
dgelist <- edgeR::estimateDisp(dgelist, designMat)
de.com <- edgeR::exactTest(dgelist)
}else{
de.com <- edgeR::exactTest(dgelist, dispersion=snakemake@config[["edger_dispersion"]])
}
} else if (snakemake@config[["edger_testType"]] == "glmLRT"){
if (snakemake@config[["edger_dispersion"]] == 0){
dgelist <- edgeR::estimateDisp(dgelist, designMat)
fit <- edgeR::glmFit(dgelist, designMat)
} else {
fit <- edgeR::glmFit(dgelist, designMat, dispersion=snakemake@config[["edger_dispersion"]])
}
de.com <- edgeR::glmLRT(fit,coef=2)
}
fit <- glmFit(dgelist, designMat)
lrt <- glmLRT(fit)
options(digits=4)
save.image(snakemake@output[["r_data"]])
write.csv(lrt$table, snakemake@output[["de_table"]])
padj<- p.adjust(de.com$table$PValue, method="bonferroni")
res <-data.frame(cbind(de.com$table$logFC/log(2),de.com$table$PValue, padj))
colnames(res) <- c("log2FoldChange", "pvalue", "padj")
rownames(res) <- names(keep)
#edgeR_result <- topTags(lrt)
rm(txi,tx2gene)
#deGenes <- decideTestsDGE(lrt, p=0.001)
#deGenes <- rownames(lrt)[as.logical(deGenes)]
#plotSmear(lrt, de.tags=deGenes)
#abline(h=c(-1, 1), col=2)
save.image(snakemake@output[["r_data"]])
write.csv(res, snakemake@output[["de_table"]])
......@@ -37,6 +37,26 @@
value: "",
label: "Annotation file (gtf or gff) : ",
},
{
name: edger_normfact,
type: radio,
selected: "TMM",
choices: [TMM: TMM, RLE: RLE, upperquartile: upperquartile, none: none],
label: "Calculate normalization factors to scale the raw library sizes.",
},
{
name: "edger_dispersion",
type: "text",
value: "0",
label: "Dispersion: either a numeric vector of dispersions or a character string indicating that dispersions should be taken from the data.",
},
{
name: edger_testType,
type: radio,
selected: "exactTest",
choices: [exactTest: exactTest, glmLRT: glmLRT],
label: "Test type: exactTest: Computes p-values for differential abundance for each gene between two samples, conditioning on the total count for each gene. The counts in each group are assumed to follow a binomial distribution\nglmLRT: Fits a negative binomial generalized log-linear model to the read counts for each gene and conducts genewise statistical tests.",
},
],
},
],
......
fastqc_threads: 1
fastqc_threads: 4
# Specifies the number of files which can be processed
# simultaneously. Each thread will be allocated 250MB of
# memory so you shouldn't run more threads than your
......
......@@ -5,7 +5,7 @@ kallisto_index_output_dir: "kallisto/index"
# kallisto quant
kallisto_quant_output_dir: "kallisto/quant"
kallisto_quant_threads: 1
kallisto_quant_threads: 4
kallisto_quant_pe_stranded: "" # can be "", "--fr-stranded" or "--rf-stranded"
## estimated fragment length and standard deviation required with single end data
......
......@@ -4,7 +4,7 @@ salmon_index_output_dir: "salmon/index"
# salmon quant
salmon_quant_output_dir: "salmon/quant"
salmon_quant_threads: 1
salmon_quant_threads: 4
salmon_outputs:
counts: salmon/quant/{sample}/quant.sf
......
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