-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathabi1_Interaction_effect.txt
119 lines (105 loc) · 5 KB
/
abi1_Interaction_effect.txt
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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
library("DESeq2")
library("gplots")
library(ggplot2)
library(pheatmap)
library("RColorBrewer")
library(scatterplot3d)
counts <- read.delim("E:/MicroRNA Data Analysis/SEQ Data analysis/3. Third Batch/2. Interaction Effect/WT_ABI1/abi_v3.txt", sep="\t", header=T, row.names=1)
counts <- as.matrix(counts)
design <- data.frame( condition=factor( c("ABA", "ABA", "ABA", "Control", "Control", "Control", "ABA", "ABA", "ABA", "Control", "Control", "Control") ) )
rownames(design) <- colnames(counts)
head(counts)
dds <- DESeqDataSetFromMatrix(countData = counts, colData = design, design = ~ condition)
dds$genotype <- factor( c("abi1", "abi1", "abi1", "abi1", "abi1", "abi1", "wt", "wt", "wt", "wt", "wt", "wt") )
dds$genotype = relevel( dds$genotype, "wt")
dds$condition = relevel( dds$condition, "Control")
design(dds) <- ~ genotype + condition + genotype:condition
dds <- DESeq(dds)
resultsNames(dds)
res1 <- results(dds, name="genotypeabi1.conditionABA")
res1 <- res1[ complete.cases(res1$padj), ]
res1 <- res1[order(res1$padj),]
new_columns <- data.frame(GeneID=rownames(res1))
res1 <- cbind( new_columns, res1)
write.csv(res1, file='abi1_interaction_effect.csv', quote=F, row.names=F)
rld <- rlog(dds, blind=FALSE)
rlogMat <- assay(rld)
pdf("PCA_1&2_abi1_interaction_effect.pdf")
pcaData <- plotPCA(rld, intgroup=c("condition", "genotype"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=condition, shape=genotype)) +
geom_point(size=3) +
xlab(paste0("PC1: ",percentVar[1],"% variance")) +
ylab(paste0("PC2: ",percentVar[2],"% variance")) +
coord_fixed()
dev.off()
select <- order(rowMeans(counts(dds,normalized=TRUE)), decreasing=TRUE)[1:20]
df <- as.data.frame(colData(dds)[,c("condition","genotype")])
pdf("pheatmap_abi1_interaction_effect.pdf")
pheatmap(assay(rld)[select,], cluster_rows=FALSE, show_rownames=FALSE,
cluster_cols=FALSE, annotation_col=df)
dev.off()
sampleDists <- dist(t(assay(rld)))
sampleDistMatrix <- as.matrix(sampleDists)
pdf("SampletoSampleDistance_abi1_interaction_effect.pdf")
rownames(sampleDistMatrix) <- paste(rld$condition, rld$genotype, sep="-")
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
clustering_distance_rows=sampleDists,
clustering_distance_cols=sampleDists,
col=colors)
dev.off()
pdf("MA_plot_abi1_interaction_effect.pdf")
plotMA(dds, main="MA plot", ylim=c(-2, 2))
#res <- lfcShrink(dds, coef="genotypeabi1.conditionABA")
#plotMA(res)
dev.off()
pdf("Dispersionplot_abi1_interaction_effect.pdf")
plotDispEsts(dds)
dev.off()
rld <- rlog(dds, blind=FALSE)
rlogMat <- assay(rld)
PCA <- prcomp(t(rlogMat))
percentVar <- round(100*PCA$sdev^2/sum(PCA$sdev^2),1)
dataGG = data.frame(PC1 = PCA$x[,1], PC2 = PCA$x[,2], PC3 = PCA$x[,3], PC4 = PCA$x[,4], condition = colData(rld)$condition, genotype=colData(rld)$genotype)
pdf("3D plot.pdf")
colors <- c("#bf0f0f", "#56B4E9")
colors <- colors[as.numeric(dataGG$condition)]
shapes = c(15, 16, 17,22)
shapes <- shapes[as.numeric(dataGG$genotype)]
s3d <- scatterplot3d(dataGG[,1:3], pch = shapes, color=colors, main = "3D Plot", angle=50)
legend("bottomright", legend = levels(dataGG$condition), pch = c(15, 16, 17, 22),col = c("#bf0f0f", "#56B4E9"), inset=.05, xpd = TRUE, horiz = TRUE)
legend("topleft", legend = levels(dataGG$genotype), pch = c(15, 16, 17,22), inset=.05, xpd = TRUE, horiz = TRUE)
dev.off()
library("pheatmap")
mat = assay(rld)[ head(order(res$padj),30), ]
mat = mat - rowMeans(mat)
df = as.data.frame(colData(rld)[,c("condition")])
colnames(df) = "condition"
rownames(df) = colnames(mat)
pheatmap(mat, annotation_col=df)
pdf("plotcounts.pdf")
res <- results(dds)
topGene <- rownames(res)[which.min(res$padj)]
data <- plotCounts(dds, gene=topGene, intgroup=c("condition","genotype"), returnData=TRUE)
ggplot(data, aes(x=condition, y=count, color=genotype, group=genotype)) +
scale_y_log10() + geom_point(size=3) + geom_line()+labs(title = topGene)+theme(plot.title = element_text(hjust = 0.5))
dev.off()
pdf("Screeplot.pdf")
barplot(percentVar, cex.names=1, xlab=paste("Principal component (PC), 1-", length(PCA$sdev)), ylab="Proportion of variation (%)", main="Scree plot", ylim=c(0,100))
dev.off()
pdf("volcanoplot.pdf")
tab = data.frame(logFC = res$log2FoldChange, negLogPval = -log10(res$pvalue))
head(tab)
par(mar = c(5, 4, 4, 4))
plot(tab, pch = 16, cex = 0.6, xlab = expression(log[2]~fold~change), ylab = expression(-log[10]~pvalue))
lfc = 2
pval = 0.05
signGenes = (abs(tab$logFC) > lfc & tab$negLogPval > -log10(pval))
points(tab[signGenes, ], pch = 16, cex = 0.8, col = "red")
abline(h = -log10(pval), col = "green3", lty = 2)
abline(v = c(-lfc, lfc), col = "blue", lty = 2)
mtext(paste("pval =", pval), side = 4, at = -log10(pval), cex = 0.8, line = 0.5, las = 1)
mtext(c(paste("-", lfc, "fold"), paste("+", lfc, "fold")), side = 3, at = c(-lfc, lfc), cex = 0.8, line = 0.5)
dev.off()