deseq2.script.R

#log <- file(snakemake@log[[1]], open="wt")
#sink(log)
#sink(log, type="message")
library(DESeq2)
library(edgeR)
library(apeglm)
library(BiocParallel)
library(tximport)
library(rhdf5)
library(GenomicFeatures)
library(ggplot2)

register(MulticoreParam(snakemake@threads))

samples = read.table(snakemake@config[["popmap_file"]],stringsAsFactors=F)$V1
conditions = read.table(snakemake@config[["popmap_file"]],stringsAsFactors=F)$V2
#samples = c("DRR016125","DRR016126","DRR016127","DRR016128","DRR016129","DRR016130","DRR016131","DRR016132")
#files = file.path("/home/mbb/Documents/results2/kallisto/quant",samples,"abundance.h5")
files = sort(snakemake@input[["counts"]])
names(files) = samples

if (snakemake@config[["quantification"]] == "htseq_count")
{
   print("Quantification was done with htseq-count !")
   sampleTable <- data.frame(
		sampleName = samples,
		fileName   = files,
		condition  = conditions
		)
    # create the DESeq data object
    dds <- DESeqDataSetFromHTSeqCount(
		sampleTable = sampleTable, 
		directory = "", #paste0(snakemake@config[["results_dir"]],"/",snakemake@config[["<step_name>__deseq2_output_dir"]],"/"), 
		design = ~ condition
		)

} else{
    if (snakemake@params[["tx2gene"]] == TRUE){
      # conversion transcrits vers gènes
      TxDb <- makeTxDbFromGFF(file = snakemake@params[["annotations"]])
      k <- keys(TxDb, keytype = "TXNAME")
      tx2gene <- select(TxDb, k, "GENEID", "TXNAME")
      txi <- tximport(files, type = snakemake@config[["quantification"]], tx2gene = tx2gene) 
    } else {
      txi <- tximport(files, type = snakemake@config[["quantification"]], txOut = TRUE) 
    }

    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)

    dds <- DESeqDataSetFromTximport(txi, sampleTable, designMat)
}

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

dds <- DESeq(dds) #, fitType = snakemake@config[["deseq2_fitType"]], betaPrior = as.logical(snakemake@config[["deseq2_betaPrior"]]), test="Wald")

res <- results(dds)

save.image(snakemake@output[["r_data"]])
write.csv(res, snakemake@output[["de_table"]])

###### PREPARE REPORT ######
# normalized, variance-stabilized transformed counts for visualization
if (sum(keep) < 1000) {
  vsd = varianceStabilizingTransformation(dds)
  } else {
    vsd <- vst(dds, blind=FALSE)
    }

dat <- plotPCA(vsd, intgroup="condition",returnData=TRUE)


png(filename = snakemake@output[["PCA"]], height=800, width=800)
 p <- ggplot(dat,aes(x=PC1,y=PC2,col=group))
 p + geom_point() + ggtitle("PCA on normalized, variance-stabilized transformed counts")
dev.off()

# cpm_counts = cpm(txi$counts)
# pca = prcomp(t(cpm_counts))

# png(filename = snakemake@output[["PCA"]], height=800, width=800)
# ggplot(data.frame(pca$x),aes(x =pca$x[,1] , y = pca$x[,2],col = conditions)) + geom_point()
# dev.off()

options(digit=4)

resSig <- subset(res, padj < 0.1)
resSigdf = data.frame(resSig@listData,row.names = resSig@rownames)

#datatable(resSigdf,options = list(scrollX = '300px'))

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[["Top_genes"]])

write.table(topGenes,snakemake@output[["Top_genes"]],append=TRUE,sep="\t")

png(filename = snakemake@output[["MA_plot"]], height=800, width=800)
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()


# Heatmap
library(pheatmap)
select <- order(rowMeans(counts(dds,normalized=TRUE)), decreasing=TRUE)[1:20]
df <- as.data.frame(colData(dds)[,c("condition")])
ntd <- normTransform(dds)
rownames(df) = samples
colnames(df) = "Conditions"
assayNTD = assay(ntd)

png(filename = snakemake@output[["Heatmap"]], height=800, width=800)
pheatmap(assayNTD[select,], cluster_rows=FALSE, show_rownames=FALSE, cluster_cols=FALSE, annotation_col=df)

#cpmTop = cpm_counts[rownames(resSigdf),]

#heatmap(cpmTop,main = "Heatmap",margins = c(10,1))
dev.off()

#resLFC <- lfcShrink(dds, coef="groups", type="apeglm")
#resOrdered <- res[order(res$pvalue),]
 
#plotMA(res, ylim=c(-2,2))
#plotMA(resLFC, ylim=c(-2,2))
# 
# ntd <- normTransform(dds)
# 
# library("vsn")
# meanSdPlot(assay(ntd))
# 
# library("pheatmap")
# select <- order(rowMeans(counts(dds,normalized=TRUE)), decreasing=TRUE)[1:20]
# df <- as.data.frame(colData(dds)[,c("condition")])
# pheatmap(assay(ntd)[select,], cluster_rows=FALSE, show_rownames=FALSE, cluster_cols=FALSE, annotation_col=df)
# 
# 
# sampleDists <- dist(t(assay(vsd)))
# 
# library("RColorBrewer")
# sampleDistMatrix <- as.matrix(sampleDists)
# rownames(sampleDistMatrix) <- paste(vsd$condition, vsd$type, sep="-")
# colnames(sampleDistMatrix) <- NULL
# colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
# pheatmap(sampleDistMatrix,
#          clustering_distance_rows=sampleDists,
#          clustering_distance_cols=sampleDists,
#          col=colors)
# 
# plotPCA(vsd, intgroup=c("condition"))
#sink(NULL, type = "message")