DESeq2_example.R 1.32 KB
Newer Older
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
# Sample DESeq2 script using data from vignette examples in Docker container
# A real data analysis script should have much more than what is shown here :)

library(DESeq2)
library(pasilla)

pasCts <- system.file("extdata",
                      "pasilla_gene_counts.tsv",
                      package="pasilla", mustWork=TRUE)

pasAnno <- system.file("extdata",
                       "pasilla_sample_annotation.csv",
                       package="pasilla", mustWork=TRUE)

cts <- as.matrix(read.csv(pasCts,sep="\t",row.names="gene_id"))
head(cts)

coldata <- read.csv(pasAnno, row.names=1)
coldata <- coldata[,c("condition","type")]
coldata

coldata$condition <- factor(coldata$condition)
coldata$type <- factor(coldata$type)

# re-ordering sample names to match b/t coldata and counts

rownames(coldata) <- sub("fb", "", rownames(coldata))
all(rownames(coldata) %in% colnames(cts))

all(rownames(coldata) == colnames(cts))

cts <- cts[, rownames(coldata)]
all(rownames(coldata) == colnames(cts)) # now TRUE

# make DESeqDataSet object

dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ condition)
dds

# PCA

rld <- rlog(dds, blind=FALSE)

plotPCA(rld, intgroup=c("condition", "type"))


# simple diff. expression

dds <- DESeq(dds)
res <- results(dds)
res