diff --git a/files/Snakefile b/files/Snakefile index ddb07525692313ba961fe778f06a107fb320687f..d2fed04c06292d50bebb2033a3c5b4e5e26d9673 100644 --- a/files/Snakefile +++ b/files/Snakefile @@ -360,6 +360,7 @@ rule edger: PCA = config["results_dir"]+'/'+config["edger_output_dir"]+'/PCA_mqc.png', MD = config["results_dir"]+'/'+config["edger_output_dir"]+'/MD_mqc.png', Top_genes = config["results_dir"]+'/'+config["edger_output_dir"]+'/Top_genes_mqc.tsv', + Volcano_plot = config["results_dir"]+'/'+config["edger_output_dir"]+'/Volcano_plot_mqc.png', Heatmap = config["results_dir"]+'/'+config["edger_output_dir"]+'/Heatmap_mqc.png', log: config["results_dir"]+'/logs/edger/edger_log.txt' script: @@ -372,8 +373,9 @@ rule deseq2: de_table = config["results_dir"]+'/'+config["deseq2_output_dir"]+'/de_table.csv', r_data = config["results_dir"]+'/'+config["deseq2_output_dir"]+'/data.RData', PCA = config["results_dir"]+'/'+config["deseq2_output_dir"]+'/PCA_mqc.png', - DE_Table = config["results_dir"]+'/'+config["deseq2_output_dir"]+'/DE_Table_mqc.tsv', + Top_genes = config["results_dir"]+'/'+config["deseq2_output_dir"]+'/Top_genes_mqc.tsv', MA_plot = config["results_dir"]+'/'+config["deseq2_output_dir"]+'/MA_plot_mqc.png', + Volcano_plot = config["results_dir"]+'/'+config["deseq2_output_dir"]+'/Volcano_plot_mqc.png', Heatmap = config["results_dir"]+'/'+config["deseq2_output_dir"]+'/Heatmap_mqc.png', log: config["results_dir"]+'/logs/deseq2/deseq2_log.txt' script: diff --git a/files/scripts/deseq2.script.R b/files/scripts/deseq2.script.R index 705bbf0b4e362e31400f01251b9a81e3686a7354..e003397ca58d13c103a65a00c0dbe7b476b226cd 100644 --- a/files/scripts/deseq2.script.R +++ b/files/scripts/deseq2.script.R @@ -29,8 +29,10 @@ if (snakemake@config[["deseq2_tx2gene"]] == TRUE){ txi <- tximport(files, type = snakemake@config[["quantification"]], txOut = TRUE) } -designMat <- model.matrix(~conditions) +groups = as.factor(conditions) +designMat <- model.matrix(~groups) colnames(designMat)[2] = "condition" +#dds$condition <- factor(dds$condition, levels = unique(conditions)) sampleTable <- data.frame(condition = factor(conditions)) rownames(sampleTable) <- colnames(txi$counts) @@ -40,8 +42,6 @@ dds <- DESeqDataSetFromTximport(txi, sampleTable, designMat) keep <- rowSums(counts(dds)) >= 10 dds <- dds[keep,] -dds$condition <- factor(dds$condition, levels = unique(conditions)) - dds <- DESeq(dds) #, fitType = snakemake@config[["deseq2_fitType"]], betaPrior = as.logical(snakemake@config[["deseq2_betaPrior"]]), test="Wald") res <- results(dds) @@ -65,22 +65,26 @@ resSigdf = data.frame(resSig@listData,row.names = resSig@rownames) #datatable(resSigdf,options = list(scrollX = '300px')) -cat("# id: DE_Table -# section_name: 'DE Table' -# description: 'Tableau d'expression différentielle avec p-value ajustée < 0.1' +topGenes = head(resSigdf[order(resSigdf$padj),], 40) + +cat("# id: Top_genes +# section_name: 'Top 40 genes' +# description: 'Tableau des 40 top genes' # format: 'tsv' -# plot_type: 'table'\ngene\t", file=snakemake@output[["DE_Table"]]) +# plot_type: 'table'\ngene\t", file=snakemake@output[["Top_genes"]]) -write.table(resSigdf,snakemake@output[["DE_Table"]],append=TRUE,sep="\t") +write.table(topGenes,snakemake@output[["Top_genes"]],append=TRUE,sep="\t") png(filename = snakemake@output[["MA_plot"]], height=800, width=800) -plotMA(res, ylim=c(-2,2)) +DESeq2::plotMA(res, ylim=c(-5,5)) dev.off() +png(filename = snakemake@output[["Volcano_plot"]], height=800, width=800) with(res, plot(log2FoldChange, -log10(pvalue), pch=20, main="Volcano plot", xlim=c(-2.5,2))) with(subset(res, padj<.05 ), points(log2FoldChange, -log10(pvalue), pch=20, col="red")) with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(pvalue), pch=20, col="orange")) with(subset(res, padj<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(pvalue), pch=20, col="green")) +dev.off() cpmTop = cpm_counts[rownames(resSigdf),] @@ -93,7 +97,6 @@ rownames(df) = samples colnames(df) = "Conditions" assayNTD = assay(ntd) -dev.off() png(filename = snakemake@output[["Heatmap"]], height=800, width=800) pheatmap(assayNTD[select,], cluster_rows=FALSE, show_rownames=FALSE, cluster_cols=FALSE, annotation_col=df) #heatmap(cpmTop,main = "Heatmap",margins = c(10,1)) diff --git a/files/scripts/edger.script.R b/files/scripts/edger.script.R index 0aa3eec20c8985e02b4ce96f1777dff88ca981d6..46a7c2656b4b4a5f68ca109d68ca54f58beff089 100644 --- a/files/scripts/edger.script.R +++ b/files/scripts/edger.script.R @@ -91,10 +91,12 @@ cat("# id: Top_genes write.table(as.data.frame(topGenes),snakemake@output[["Top_genes"]],append=TRUE,sep="\t") +png(filename = snakemake@output[["Volcano_plot"]], height=800, width=800) with(res, plot(log2FoldChange, -log10(pvalue), pch=20, main="Volcano plot", xlim=c(-2.5,2))) with(subset(res, padj<.05 ), points(log2FoldChange, -log10(pvalue), pch=20, col="red")) with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(pvalue), pch=20, col="orange")) with(subset(res, padj<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(pvalue), pch=20, col="green")) +dev.off() cpmTop = cpm_counts[rownames(topGenes$table),] @@ -103,7 +105,6 @@ df <- as.data.frame(conditions) rownames(df) = samples colnames(df) = "Conditions" -dev.off() png(filename = paste(snakemake@output[["Heatmap"]],sep = "/"), height=800, width=800) pheatmap(cpmTop, cluster_rows=FALSE, show_rownames=FALSE, cluster_cols=FALSE, annotation_col=df) #heatmap(cpmTop)