-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy paths3_Heatmap.R
164 lines (130 loc) · 5.11 KB
/
s3_Heatmap.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
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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
library(ggplot2)
library(pheatmap)
library(dplyr)
height <- function(df){
if(nrow(df) >= 1 & nrow(df) <= 20 ){
height=7
} else if(nrow(df) > 20 & nrow(df) <= 60 ) {
height=10
} else if(nrow(df) > 60 & nrow(df) <= 200 ) {
height=15
} else if (nrow(df) > 200 ) {
height=140
}
return(height)
}
save_pheatmap_pdf <- function(x, filename, width=7, height=7) {
stopifnot(!missing(x))
stopifnot(!missing(filename))
pdf(filename, width=width, height=height)
grid::grid.newpage()
grid::grid.draw(x$gtable)
dev.off()
}
EXP <- read.csv('/home/disk/ZF/mTOR_0720/res/FPKM/combined_fpkm.csv')
#op="/home/disk/ZF/mTOR_0720/res/Heatmap"
op="/home/disk/ZF/mTOR_0720/res/Heatmap_noIg"
FC <- 1
Qvalue <- 0.1
CUT=paste('FC',FC,'_Q',Qvalue,sep='')
#DEG_file='DEG_new'
DEG_file='DEG_new_noIg'
## 1. 两两比
cond=c('WT-DMSO','WT-200nMAZD','WT-500nMAZD','KI-DMSO','KI-200nMAZD','KI-500nMAZD')
combn_lst <- t(combn(cond,2))
if(!file.exists(file.path(op,CUT))){
dir.create(file.path(op,CUT),recursive = TRUE)
}
## offtarget
for(i in 13:nrow(combn_lst)){
cond1 <- combn_lst[i,1]
cond2 <- combn_lst[i,2]
TAG <- paste(cond1,cond2,sep='_')
sigUP_lst <- read.csv(paste('/home/disk/ZF/mTOR_0720/res/',DEG_file,'/',CUT,'/',cond1,'_',cond2,'/sigUP_gene.csv',sep=''))[,1]
sigDN_lst <- read.csv(paste('/home/disk/ZF/mTOR_0720/res/',DEG_file,'/',CUT,'/',cond1,'_',cond2,'/sigDN_gene.csv',sep=''))[,1]
df <- EXP[EXP$gene_short_name %in% c(sigUP_lst,sigDN_lst) ,
c(colnames(EXP)=='gene_short_name',grep(sub('-','_',cond1),colnames(EXP)),grep(sub('-','_',cond2),colnames(EXP)))]
print(dim(df))
rownames(df) <- df$gene_short_name
df <- df[,-1]
h=height(df)
tmp=df
stmp=apply((tmp+1),2,log,2)
stmp=t(apply(tmp,1,scale))
rownames(stmp)=rownames(df)
colnames(stmp)=colnames(df)
this_diff=apply(stmp[,1:3],1,mean)-apply(stmp[,4:6],1,mean)
stmp=stmp[order(this_diff),]
bk <- c(seq(-3,-0.1,by=0.01),seq(0,3,by=0.01))
p <- pheatmap(stmp,scale = "none",cexRow=0.5,fontsize=8,fontsize_row = 8,fontsize_col=8,cluster_col = F,
cluster_row = F,treeheight_row=0, treeheight_col=0,border_color=NA,
color = colorRampPalette(colors = c("yellow2","blue2"))(400),
breaks=seq(-2,2,0.01),
legend_labels = c( "-2", "-1", "0", "1", "2"))
save_pheatmap_pdf(p, paste(op,'/',CUT, '/',TAG,'.pdf',sep=''),4,height = h)
}
### 2. 3列比
# FC <- 2
# Qvalue <- 0.05
# CUT=paste('FC',FC,'_Q',Qvalue,sep='')
sigUP_lst <- read.csv(paste('/home/disk/ZF/mTOR_0720/res/',DEG_file,'/',CUT,'/KI200_500_overlap/sigUP_gene.csv',sep=''))[,1]
sigDN_lst <- read.csv(paste('/home/disk/ZF/mTOR_0720/res/',DEG_file,'/',CUT,'/KI200_500_overlap/sigDN_gene.csv',sep=''))[,1]
df <- EXP[EXP$gene_short_name %in% c(sigUP_lst,sigDN_lst) ,
c(colnames(EXP)=='gene_short_name',
grep('KI_DMSO',colnames(EXP)),grep('KI_200nMAZD',colnames(EXP)),grep('KI_500nMAZD',colnames(EXP)))]
rownames(df) <- df$gene_short_name
df <- df[,-1]
df <- as.data.frame(t(df))
genes <- c(sigUP_lst,sigDN_lst)
df <- as.data.frame(t(select(df,genes)))
dim(df)
h=10
tmp=df
stmp=apply((tmp+1),2,log,2)
stmp=t(apply(tmp,1,scale))
rownames(stmp)=rownames(df)
colnames(stmp)=colnames(df)
p <- pheatmap(stmp,scale = "none",cexRow=0.5,fontsize=8,fontsize_row = 8,fontsize_col=8,cluster_col = F,
cluster_row = F,treeheight_row=0, treeheight_col=0,border_color=NA,
color = colorRampPalette(colors = c("yellow2","blue2"))(400),
breaks=seq(-2,2,0.01),
legend_labels = c( "-2", "-1", "0", "1", "2"))
save_pheatmap_pdf(p, paste(op,'/',CUT, '/KI200_500_overlap.pdf',sep=''),6,height = h)
### 4. on-target gene
op_DEG=paste('/home/disk/ZF/mTOR_0720/res/',DEG_file,'/',sep='')
#op_heatmap='/home/disk/ZF/mTOR_0720/res/Heatmap/'
op_heatmap='/home/disk/ZF/mTOR_0720/res/Heatmap_noIg/'
AZD=500
# FC=1.3
# Qvalue=0.05
# CUT=paste('FC',FC,'_Q',Qvalue,sep='')
if(AZD==200){
TAG='WT_KI_200AZD_overlap'
}else if(AZD==500){
TAG='WT_KI_500AZD_overlap'
}
TYPE <- c('ALL','UP','DN')
for(i in 1:length(TYPE)){
type=TYPE[i] #DN/UP
sigG_lst <- read.csv(paste(op_DEG,'FC',FC,'_Q',Qvalue,'/',TAG,'/sig',type,'_gene.csv',sep=''))[,1]
df <- EXP[EXP$gene_short_name %in% sigG_lst ,
c(colnames(EXP)=='gene_short_name',
grep('KI_DMSO',colnames(EXP)),grep('KI_200nMAZD',colnames(EXP)),grep('KI_500nMAZD',colnames(EXP)))]
rownames(df) <- df$gene_short_name
df <- df[,-1]
df <- as.data.frame(t(df))
df <- as.data.frame(t(select(df,sigG_lst)))
dim(df)
h=height(df)
tmp=df
stmp=apply((tmp+1),2,log,2)
stmp=t(apply(tmp,1,scale))
rownames(stmp)=rownames(df)
colnames(stmp)=colnames(df)
p <- pheatmap(stmp,scale = "none",cexRow=0.5,fontsize=8,fontsize_row = 8,fontsize_col=8,cluster_col = F,
cluster_row = F,treeheight_row=0, treeheight_col=0,border_color=NA,
color = colorRampPalette(colors = c("yellow2","blue2"))(400),
breaks=seq(-2,2,0.01),
legend_labels = c("-2", "-1", "0", "1", "2"))
save_pheatmap_pdf(p, paste(op_heatmap,'/',CUT, '/',TAG,'_',type,'.pdf',sep=''),6,height = h)
}