This repository has been archived by the owner on Jul 30, 2019. It is now read-only.
forked from jpverta/verta_jones_elife_2019
-
Notifications
You must be signed in to change notification settings - Fork 0
/
DEanalysisPureEcotypes.r
83 lines (66 loc) · 3.22 KB
/
DEanalysisPureEcotypes.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
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
library(gtools)
library(gplots)
#--------------------------------
qThreshold = 0.2
#--------------------------------
# DE test results - TYNE
tyneDEtest = read.table('.../expression/cuffdiff/STAR/fw_vs_mar_Tyne_UCSCgtf/gene_exp.diff',header=T,stringsAsFactors=F,sep='\t')
tyneDEsig = tyneDEtest[which(tyneDEtest$q_value<qThreshold),]
tyneDE = tyneDEsig$gene_id
tyneFC = tyneDEtest$log2.fold_change.
names(tyneFC) = tyneDEtest$gene_id
tyneLR = foldchange2logratio(tyneFC)
#--------------------------------
# DE test results - LITC
litcDEtest = read.table('.../expression/cuffdiff/STAR/fw_vs_mar_Litc_UCSCgtf/gene_exp.diff',header=T,stringsAsFactors=F,sep='\t')
litcDEsig = litcDEtest[which(litcDEtest$q_value<qThreshold),]
litcDE = litcDEsig$gene_id
litcFC = litcDEtest$log2.fold_change.
names(litcFC) = litcDEtest$gene_id
litcLR = foldchange2logratio(litcFC)
#--------------------------------
# compare Tyne and Litc
v = venn(list(tyneDE,litcDE))
sharedSig = attr(v,"intersections")$"A:B"
litcSpecific = attr(v,"intersections")$"B"
tyneSpecific = attr(v,"intersections")$"A"
# plot
plot(litcFC,tyneFC,pch=20,cex=0.2)
points(litcFC[litcSpecific],tyneFC[litcSpecific],col='blue')
points(litcFC[tyneSpecific],tyneFC[tyneSpecific],col='green')
points(litcFC[sharedSig],tyneFC[sharedSig],col='red')
# parallel
parallel = sharedSig[which(sign(litcFC[sharedSig]) == sign(tyneFC[sharedSig]))]
# anti-parallel
antiParallel = sharedSig[which(sign(litcFC[sharedSig]) != sign(tyneFC[sharedSig]))]
#--------------------------------
# DE test results - both together
parallelDEtest = read.table('.../expression/cuffdiff/STAR/fw_vs_mar_Litc_Tyne_UCSCgtf/gene_exp.diff',header=T,stringsAsFactors=F,sep='\t')
rownames(parallelDEtest) = parallelDEtest$gene_id
parallelDEsig = parallelDEtest[which(parallelDEtest$q_value<qThreshold),]
parallelDE = parallelDEsig$gene_id
# parallel
parallel = parallelDE[which(sign(litcFC[parallelDE]) == sign(tyneFC[parallelDE]))]
# plot
plot(litcFC,tyneFC,pch=20,cex=0.2)
points(litcFC[litcDE],tyneFC[litcDE],col='blue')
points(litcFC[tyneDE],tyneFC[tyneDE],col='green')
points(litcFC[parallel],tyneFC[parallel],col='red')
# compPC
sameSign = names(litcLR)[which(sign(litcLR) == sign(tyneLR))]
PC2 = read.table('.../ASE/UCSCgtf/Tyne_Litc_PC2_loadings.txt')
PC5 = read.table('.../ASE/UCSCgtf/Tyne_Litc_PC5_loadings.txt')
PC = data.frame(x=0.145*PC2$x+0.063*-1*PC5$x)
rownames(PC) = rownames(PC2)
extremePCglobal = rownames(PC)[which(PC$x<quantile(PC$x,probs=c(0.01)) | PC$x>quantile(PC$x,probs=c(0.99)))]
PC$zPC = scale(PC$x)
extremePCglobalSameSign = extremePCglobal[which(extremePCglobal %in% sameSign)]
# plot
plotData = data.frame(litc=litcFC,tyne=tyneFC,parallel=names(litcLR)%in%parallel,PC=names(litcFC)%in%extremePCglobal)
g = ggplot(plotData[which(plotData$parallel == F),],aes(litc,tyne))
g = g + geom_point(color='grey')
g = g + geom_point(data=plotData[which(plotData$parallel == T),],aes(litc,tyne),color='black')
g = g + theme(axis.text=element_text(hjust=1,size=24,colour='black'),panel.background=element_rect(fill='white',colour='black'),text=element_text(size=22),legend.position="right")
#g = g + geom_point(data=plotData[which(plotData$PC == T),],aes(litc,tyne),color='blue')
g
ggsave(".../ASE/UCSCgtf/LitcTyneParallelDE.ps", device=cairo_ps)