-
Notifications
You must be signed in to change notification settings - Fork 9
/
RUVSeq.R
39 lines (29 loc) · 1.7 KB
/
RUVSeq.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
37
38
39
library(DESeq2)
library(RUVSeq) # http://bioconductor.org/packages/release/bioc/vignettes/RUVSeq/inst/doc/RUVSeq.pdf
library(EDASeq)
setwd('/Volumes/ncatssctl/NGS_related/BulkRNA/Vukasin_SRA/')
sampleTable <- fread('sampletable_v5set.csv')
load('Vukasin_SRAv5.DDS.RData') # dds is loaded
x <- as.factor(sampleTable$condition)
set <- newSeqExpressionSet(as.matrix(counts(dds, normalized=F)),
phenoData = data.frame(x, row.names= c(sampleTable$sampleFiles) ))
set
library(RColorBrewer)
colors <- brewer.pal(9, "Set2")
# before correction:
plotRLE(set, outline=FALSE, ylim=c(-4, 4), col=colors[x])
plotPCA(set, col=colors[x], cex=1.2)
spikes <- unique(c('ANAPC5', 'ANAPC15', 'ARID3B', 'ARL10', 'ATXN2', 'C16orf62', 'C3orf49', 'CCAR1', 'CCDC125', 'CCDC90B', 'CHFR', 'DHRSX', 'FRMD8', 'GGA1', 'HERC4', 'MKNK1', 'NASP', 'NME4', 'OTUB1', 'PMF1', 'POLR2B', 'POLR3A', 'POMK', 'PSMA3-AS1', 'PTPN14', 'RAPGEF6', 'REL', 'RRP1', 'RUNDC1', 'SAMD4B', 'SLC4A1AP', 'SLMAP', 'SMARCAL1', 'SNAP29', 'SNRNP200', 'SUPT4H1', 'TBC1D22A', 'THUMPD3-AS1', 'TSPOAP1-AS1', 'TUBGCP2', 'WDTC1', 'ZNF544','C1orf43','CHMP2A','EMC7','GPI','PSMB2','PSMB4','RAB7A','REEP5','SNRPD3','VCP','VPS29')) #https://bmcmedgenomics.biomedcentral.com/articles/10.1186/s12920-019-0538-z#Tab3
spikes <- spikes[spikes %in% row.names(dds)]
set2 <- RUVg(set, spikes, k=1)
pData(set2)
# after correction:
plotRLE(set2, outline=FALSE, ylim=c(-4, 4), col=colors[x])
plotPCA(set2, col=colors[x], cex=1.2)
condition <- x
names(pData(set2)) <- c('condition', 'W_1')
dds <- DESeqDataSetFromMatrix(countData = counts(set2),
colData = pData(set2),
design = ~ W_1 + condition)
dds <- DESeq(dds)
# save(...)