Commit 69d2a16c authored by khalid's avatar khalid
Browse files

Add GO and functional enrichment to edgeR

parent 9e373320
......@@ -14,5 +14,6 @@ rule <step_name>__deseq2:
script = config["<step_name>__deseq2_command"],
tx2gene = config["<step_name>__deseq2_tx2gene"],
annotations = config["<step_name>__deseq2_annotations"],
biomaRtDataset = = config["<step_name>__deseq2_biomaRtDataset"],
script:
config["<step_name>__deseq2_script"]
\ No newline at end of file
......@@ -133,13 +133,51 @@ 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)
pheatmap(assayNTD[select,], cluster_rows=FALSE, show_rownames=FALSE, cluster_cols=TRUE, annotation_col=df)
#heatmap(cpmTop,main = "Heatmap",margins = c(10,1))
dev.off()
#cpmTop = cpm_counts[rownames(resSigdf),]
speciesData = snakemake@config[["biomaRtDataset"]]
#heatmap(cpmTop,main = "Heatmap",margins = c(10,1))
if(speciesData != "None")
{
#Add annotation with biomaRt
library(biomaRt)
mart <- useMart(biomart="ENSEMBL_MART_ENSEMBL", host="ensembl.org")
species_mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", host = "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,paste(snakemake@output[["Top_genes"]],"_mqc.tsv",sep="\t")
#test for gene enrichment with gprofiler2
library(gprofiler2)
#for the moment the top genes change to significant genes ?
organisme = gsub("_gene_ensembl", "", speciesData)
gostres <- gost(query = rownames(topGenes),
organism = organisme, ordered_query = FALSE,
multi_query = FALSE, significant = FALSE, exclude_iea = FALSE,
measure_underrepresentation = FALSE, evcodes = FALSE,
user_threshold = 0.05, correction_method = "g_SCS",
domain_scope = "annotated", custom_bg = NULL,
numeric_ns = "", sources = NULL, as_short_link = FALSE)
png(filename = paste0(snakemake@config[["results_dir"]],"/",snakemake@config[["<step_name>__deseq2_output_dir"]],"/gprofiler2_mqc.png"), height=800, width=800)
gostplot(gostres, capped = TRUE, interactive = TRUE)
dev.off()
write.table(gostres$result, paste0(snakemake@config[["results_dir"]],"/",snakemake@config[["<step_name>__deseq2_output_dir"]],"/gprofiler2_mqc.tsv"), sep="\t")
}
#cpmTop = cpm_counts[rownames(resSigdf),]
#resLFC <- lfcShrink(dds, coef="groups", type="apeglm")
#resOrdered <- res[order(res$pvalue),]
......
......@@ -54,6 +54,19 @@
value: "",
label: "Annotation file (gtf or gff) : ",
},
{
name: "deseq2_biomaRtDataset",
type: select,
choices: [
None: "None",
Mus musculus: "mm_gene_ensembl",
Homo sapiens: "hsapiens_gene_ensembl",
Nile tilapia: "oniloticus_gene_ensembl",
oniloticus_gene_ensembl: "dmelanogaster_gene_ensembl"
],
value: "None",
label: "Species"
},
# {
# name: deseq2_fitType,
# type: radio,
......
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