deseq2.script.R

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



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 = "", 
		design = ~ condition
		)

} else{ # counts where generated with one of these software "salmon", "sailfish", "alevin", "kallisto", "rsem", "stringtie"
    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"

    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)
percentVar <- round(100 * attr(dat, "percentVar"))

png(filename = snakemake@output[["PCA"]], height=800, width=800)
 p <- ggplot(dat,aes(x=PC1,y=PC2,col=group))
 p + geom_point(size=3) + xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed() + ggtitle("PCA on normalized, variance-stabilized transformed counts")
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)
#top 40 genes basé sur les p-values
topGenes = head(as.data.frame(res[order(res$pvalue),]) , 40)

cat("# id: Top_genes
# section_name: 'Top 40 genes'
# description: 'Tableau des 40 top genes avec les plus faibles pvalues'
# 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=TRUE, annotation_col=df)
dev.off()

speciesData = snakemake@params[["biomaRtDataset"]]

if (speciesData != "None") {
  # Add annotation with biomaRt
  library(biomaRt)
  tryCatch({
  mart <- useMart(biomart="ENSEMBL_MART_ENSEMBL", host="uswest.ensembl.org")

  species_mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", host = "uswest.ensembl.org", dataset = speciesData)

  #for the moment the top genes change to significant genes ?
  #what if we have tarnscript ID instead of gene ID ?
  ann <- getBM(filter="ensembl_gene_id",value=rownames(topGenes),attributes=c("ensembl_gene_id","description","transcript_length"),mart=species_mart)

  go_ann <- getBM(filter="ensembl_gene_id",value=rownames(topGenes),attributes=c("ensembl_gene_id","description","go_id","name_1006","namespace_1003"),mart=species_mart)
  write.table(go_ann,paste0(snakemake@params[["outDir"]],"/TopGenes_go_ann_mqc.tsv"),sep="\t",row.names =FALSE, quote = FALSE)
  },
  error = function(err) {
    print("Can't access Ensembl Biomart !")
  }
  )
  #test for gene enrichment with gprofiler2
  library(gprofiler2)
  #for the moment the top genes change to significant genes ?
  organisme = gsub("_gene_ensembl", "", speciesData)
  #limit to GO sources

  sources = c("GO:BP", "GO:MF", "GO:CC")
  query = rownames(topGenes) #setNames(as.list(rownames(topGenes)), rownames(topGenes))
  gostres <- gost(query = query, 
                  organism = organisme, ordered_query = FALSE, 
                  multi_query = FALSE, significant = FALSE, exclude_iea = TRUE, 
                  measure_underrepresentation = FALSE, evcodes = FALSE, 
                  user_threshold = 0.05, correction_method = "g_SCS", 
                  domain_scope = "annotated", custom_bg = NULL, 
                  numeric_ns = "", sources = sources, as_short_link = FALSE)

  png(filename = paste0(snakemake@params[["outDir"]],"/gprofiler2_plot_mqc.png"), height=800, width=800)
    p <- gostplot(gostres, capped = FALSE, interactive = FALSE)
    print(p  + ggtitle("gprofiler2 enrichment analysis for top genes to find over-representation of terms from Gene Ontology"))
  dev.off()

  
  #la colonne parents est une liste qu'il faut transformer avant write.table !!
  #write.table(apply(gostres$result,2,as.character), paste0(snakemake@params[["outDir"]],"/gprofiler2_mqc.tsv"), sep="\t",row.names =FALSE, quote = FALSE)
  #ou
  ttable = cbind(1:nrow(gostres$result), gostres$result[,-14])
  colnames(ttable)[1]="id"
  write.table(ttable, paste0(snakemake@params[["outDir"]],"/gprofiler2_mqc.tsv"), sep="\t",row.names =FALSE, quote = FALSE)

}