forked from adeslatt/deseq2-docker
-
Notifications
You must be signed in to change notification settings - Fork 1
/
DESeq2.R
36 lines (27 loc) · 1.04 KB
/
DESeq2.R
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
library(DESeq2)
library(tidyverse)
library(magrittr)
library(genefilter)
#take input from combined matrix
matrix <- commandArgs(trailingOnly=TRUE)
data <- read.csv(matrix, sep="\t")
de_input <- as.matrix(data[,-1])
row.names(de_input) <- data$gene_id
#object construction
dds <- DESeqDataSetFromMatrix(countData = round(de_input),
colData = data.frame(colnames(de_input)),
design = ~1)
#standard DESeq analysis
dds <- DESeq(dds)
#varianceStabilizingTransformation returns a DESeqTransform if a DESeqDataSet was provided,
#or returns a a matrix if a count matrix was provided
vsd <- varianceStabilizingTransformation(dds)
#returns a matrix
wpn_vsd <- getVarianceStabilizedData(dds)
#variance estimates for each row
rv_wpn <- rowVars(wpn_vsd)
summary(rv_wpn)
q75_wpn <- quantile( rowVars(wpn_vsd), .75) # <= original
q95_wpn <- quantile( rowVars(wpn_vsd), .95) # <= changed to 95 quantile to reduce dataset
expr_normalized <- wpn_vsd[ rv_wpn > q95_wpn, ]
write.csv(expr_normalized, "DESeq.csv")