-
Notifications
You must be signed in to change notification settings - Fork 1
/
tatus
87 lines (87 loc) · 4.17 KB
/
tatus
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
[1mdiff --git a/STARlight.smk b/STARlight.smk[m
[1mindex 1b0e06e..afc6f6e 100644[m
[1m--- a/STARlight.smk[m
[1m+++ b/STARlight.smk[m
[36m@@ -189,6 +189,7 @@[m [mrule edgeR:[m
cpm = config['cpm'],[m
nsamp = config['nsamples'],[m
logFC = config['logFC'],[m
[32m+[m[32m FDR = config['FDR'],[m
pval = config['pval'],[m
path = RESULTS,[m
edgeR = config['edgeR'],[m
[36m@@ -196,5 +197,5 @@[m [mrule edgeR:[m
24[m
shell:[m
"""[m
[31m- Rscript {params.edgeR} {input.genes} {input.sampleInfo} {input.compsTab} {params.species} {params.cpm} {params.nsamp} {params.logFC} {params.pval} {params.path}[m
[32m+[m[32m Rscript {params.edgeR} {input.genes} {input.sampleInfo} {input.compsTab} {params.species} {params.cpm} {params.nsamp} {params.logFC} {params.FDR} {params.pval} {params.path}[m
"""[m
[1mdiff --git a/config.json b/config.json[m
[1mindex 64d5c53..8cce5ed 100644[m
[1m--- a/config.json[m
[1m+++ b/config.json[m
[36m@@ -11,10 +11,11 @@[m
"sampleInfo": "/mnt/WD2/BIN22P009_CE/input/sampleInfo.txt",[m
"compsTab": "/mnt/WD2/BIN22P009_CE/input/compsTab.txt",[m
"species": "hsapiens_gene_ensembl",[m
[31m- "cpm": "1",[m
[32m+[m[32m #filter on poorly expressed genes[m
[32m+[m[32m "cpm": "10",[m
"nsamples": "3",[m
"logFC": "1.5",[m
[31m- "FDR": "0.1",[m
[32m+[m[32m "FDR": "0.05",[m
#pval could be only 0 or 1.[m
#when set to 1 filter on PValue instead of FDR.[m
"pval": "0",[m
[1mdiff --git a/scripts/edger.R b/scripts/edger.R[m
[1mindex c04b5bf..04790e6 100644[m
[1m--- a/scripts/edger.R[m
[1m+++ b/scripts/edger.R[m
[36m@@ -1,6 +1,6 @@[m
#!/usr/bin/env Rscript[m
[m
[31m-#Rscript edger.R ./Results/featureCounts_genes_mod.txt ./Results/sampleInfo.txt ./Results/compsTab.txt mmusculus_gene_ensembl 1 2 1.5 0.05 /path/to/results[m
[32m+[m[32m#Rscript edger.R ./Results/featureCounts_genes_mod.txt ./Results/sampleInfo.txt ./Results/compsTab.txt mmusculus_gene_ensembl 1 2 1.5 0.1 0 /path/to/results[m
args = commandArgs(trailingOnly=TRUE)[m
[m
# Load the libraries[m
[36m@@ -69,20 +69,21 @@[m [medgeR2anno <- function(etObj, annoObj, fileName) {[m
}[m
[m
print_tab <- function(comps, compT) {[m
[32m+[m [32mcol <- ifelse(fpval == 0, quote(FDR), quote(PValue))[m
for (i in comps$comp){[m
dt <- compT[[i]][m
write.table(dt, file=file.path(out, paste(i, "_unfiltered.tsv", sep="")), row.names=FALSE, col.names=TRUE, quote=FALSE, sep="\t")[m
[31m- write.xlsx(dt, file=file.path(out, paste(i, "_unfiltered.xlsx", sep="")), overwrite=TRUE, asTable=TRUE)[m
[31m- dt <- dt[(logFC >= flogFC & PValue <= fpval) | (logFC <= -flogFC & PValue <= fpval) ,][m
[31m- write.table(dt, file=file.path(out, paste(i, "_filtered.tsv", sep="")), row.names=FALSE, col.names=TRUE, quote=FALSE, sep="\t")[m
[31m- write.xlsx(dt, file=file.path(out, paste(i, "_filtered.xlsx", sep="")), overwrite=TRUE, asTable=TRUE)[m
[32m+[m [32mwrite.xlsx(dt, file=file.path(out, paste(i, "_unfiltered.xlsx", sep="")), overwrite=TRUE, asTable=TRUE)[m
[32m+[m[32m dt <- dt[(logFC >= flogFC & eval(col) <= fFDR) | (logFC <= -flogFC & eval(col) <= fFDR) ,][m
[32m+[m [32mwrite.table(dt, file=file.path(out, paste(i, "_filtered.tsv", sep="")), row.names=FALSE, col.names=TRUE, quote=FALSE, sep="\t")[m
[32m+[m [32mwrite.xlsx(dt, file=file.path(out, paste(i, "_filtered.xlsx", sep="")), overwrite=TRUE, asTable=TRUE)[m
}[m
}[m
[m
[m
# Create the results folder[m
#out <- file.path(getwd(), "edgeR_results")[m
[31m-out <- file.path(args[9], "edgeR_results")[m
[32m+[m[32mout <- file.path(args[10], "edgeR_results")[m
test_output(out)[m
[m
[m
[36m@@ -97,7 +98,9 @@[m [mspecies <- args[4][m
nCPM <- as.numeric(args[5])[m
nSamples <- as.numeric(args[6])[m
flogFC <- as.numeric(args[7])[m
[31m-fpval <- as.numeric(args[8])[m
[32m+[m[32mfFDR <- as.numeric(args[8])[m
[32m+[m[32mfpval <- as.numeric(args[9])[m
[32m+[m
[m
# Create the DGE list[m
dgeFull <- DGEList(data_raw, genes=rownames(data_raw), group=sampleInfo$condition)[m