Commit cf84f6d1 authored by khalid's avatar khalid
Browse files

debug deseq2 R script

parent 709cc54b
......@@ -2,11 +2,9 @@
#sink(log)
#sink(log, type="message")
library(DESeq2)
library(edgeR)
library(apeglm)
library(BiocParallel)
library(tximport)
library(rhdf5)
library(GenomicFeatures)
library(ggplot2)
......@@ -137,41 +135,40 @@ pheatmap(assayNTD[select,], cluster_rows=FALSE, show_rownames=FALSE, cluster_col
#heatmap(cpmTop,main = "Heatmap",margins = c(10,1))
dev.off()
speciesData = snakemake@config[["biomaRtDataset"]]
speciesData = snakemake@params[["biomaRtDataset"]]
print(speciesData)
if (speciesData != "None") {
# Add annotation with biomaRt
library(biomaRt)
mart <- useMart(biomart="ENSEMBL_MART_ENSEMBL", host="ensembl.org")
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()
species_mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", host = "ensembl.org", dataset = )
#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")
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),]
......
......@@ -59,7 +59,7 @@
type: select,
choices: [
None: "None",
Mus musculus: "mm_gene_ensembl",
Mus musculus: "mmusculus_gene_ensembl",
Homo sapiens: "hsapiens_gene_ensembl",
Nile tilapia: "oniloticus_gene_ensembl",
oniloticus_gene_ensembl: "dmelanogaster_gene_ensembl"
......@@ -92,7 +92,8 @@
tximport: ['Rscript -e ''BiocManager::install("tximport", version = "3.10",Ncpus=8, clean=TRUE);library("tximport")'''],
rhdf5: ['Rscript -e ''BiocManager::install("rhdf5", version = "3.10",Ncpus=8, clean=TRUE);library("rhdf5")'''],
GenomicFeatures: ['Rscript -e ''BiocManager::install("GenomicFeatures", version = "3.10",Ncpus=8, clean=TRUE);library("GenomicFeatures")'''],
pheatmap: ['Rscript -e ''install.packages("pheatmap",Ncpus=8, clean=TRUE);library("pheatmap")''']
pheatmap: ['Rscript -e ''install.packages("pheatmap",Ncpus=8, clean=TRUE);library("pheatmap")'''],
certificates: [apt-get install -y ca-certificates]
},
script: deseq2.script.R,
citations: {
......
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