rm(list = ls())
options(stringsAsFactors = F)

#step1 get genes_expr
library(AnnoProbe)
#suppressPackageStartupMessages(library(GEOquery))
#gset=AnnoProbe::geoChina('GSE6364')

a=read.table('rifpatient_tissue.txt',sep='\t',quote = "",fill = T,
comment.char = "!",header = T) # 提取表达矩阵

rownames(a)=a[,8]
a <- a[,-1]

expMatrix <- a[,1:6]
fpkmToTpm <- function(fpkm)
{
exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
tpms <- apply(expMatrix,2,fpkmToTpm)
tpms[1:3,]
colSums(tpms)

genes_expr<-tpms
group_list=c(rep('CTL',3),rep('RIF',3))
save(genes_expr,group_list,file = 'step1-genes_expr.Rdata')

load('step1-genes_expr.Rdata')

temps=as.data.frame(genes_expr)
temps$id=row.names(temps)
b=a[,7:10]
temps=merge.data.frame(temps,b,by.x="id",by.y = 'symbol')
head(temps)
rownames(temps)=temps[,1]
temps <- temps[,-1]
degexpr=temps
save(degexpr,group_list,file = 'deg-genes_expr.Rdata')
deg=degexpr[,7:9]

head(deg)
save(deg,file = 'deg.Rdata')

rm(list = ls())
options(stringsAsFactors = F)

load('deg.Rdata')
load('step1-genes_expr.Rdata')
head(deg)

visualization

need_deg=data.frame(symbols=rownames(deg), logFC=deg$logFC, p=deg$Pvalue)
deg_volcano(need_deg,1)
deg_volcano(need_deg,2)

deg_heatmap(deg,genes_expr,group_list)
deg_heatmap(deg,genes_expr,group_list,30)

check_diff_genes('GPX4',genes_expr,group_list)
check_diff_genes('FOXO1',genes_expr,group_list)
check_diff_genes('ACSL4',genes_expr,group_list)
check_diff_genes('MAP3K4',genes_expr,group_list)
check_diff_genes(c('GPX4','ACSL4','LPCAT3'),genes_expr,group_list)

#tonglujiyinbianhua
if(F){
#ferroptosis
library(KEGGREST)
exprSet<-genes_expr
#cg <- KEGGREST::keggGet("hsa04216")[[1]]$GENE
cg <- KEGGREST::keggGet("hsa04350")[[1]]$GENE
cg=as.character(sapply(cg[seq(2,length(cg),by=2)], function(x) strsplit(x,';')[[1]][1]))
check_diff_genes( cg ,exprSet,group_list)
p1<-check_diff_genes( cg ,exprSet,group_list)
p1

ggsave("RIFCTL TGFb.png",p1,height = 10, width = 5)

cg=cg[1:88]
ferroptosisrelategene=cg
ferrodeg=a[ferroptosisrelategene,]
save(ferrodeg,file = 'ferrodeg.Rddta')
check<-genes_expr[cg,]
check<-na.omit(genes_expr[cg,])
ac=data.frame(g=group_list)
rownames(ac)=colnames(check)
pheatmap(check,show_colnames =F,show_rownames = T) #对dat按照cg取行，所得到的矩阵来画热图

n=t(scale(t(check)))#通过“scale”对log-ratio数值进行归一化，现在的dat是行名为探针，列名为样本名，由于scale这个函数应用在不同组数据间存在差异时，需要行名为样本，因此需要用t(dat[cg,])来转换，最后再转换回来
n=t(scale(t(check)))#通过“scale”对log-ratio数值进行归一化，现在的dat是行名为探针，列名为样本名，由于scale这个函数应用在不同组数据间存在差异时，需要行名为样本，因此需要用t(dat[cg,])来转换，最后再转换回来

n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = T)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n) #将ac的行名也就分组信息（是‘no TNBC’还是‘TNBC’）给到n的列名，即热图中位于上方的分组信息
pheatmap(n,show_colnames =F,
show_rownames = T,
cluster_cols = F,
annotation_col=ac,
filename = 'Pre ESE ferroptosis deg.png',width = 8, height= 9) #列名注释信息为ac即分组信息

#autophagy
library(KEGGREST)
cg <- KEGGREST::keggGet("hsa04140")[[1]]$GENE
cg=as.character(sapply(cg[seq(2,length(cg),by=2)], function(x) strsplit(x,';')[[1]][1]))
check_diff_genes( cg ,gene_exCPEMS,group_list)

#apotosis
library(KEGGREST)
cg <- KEGGREST::keggGet("hsa04210")[[1]]$GENE
cg=as.character(sapply(cg[seq(2,length(cg),by=2)], function(x) strsplit(x,';')[[1]][1]))
check_diff_genes( cg ,gene_exCPEMS,group_list)

细胞衰老

map04218
library(KEGGREST)
cg <- KEGGREST::keggGet("map04933")[[1]]$GENE
cg=as.character(sapply(cg[seq(2,length(cg),by=2)], function(x) strsplit(x,';')[[1]][1]))
check_diff_genes( cg ,genes_expr,group_list)
}

for volcano

if(F){
nrDEG=deg
head(nrDEG)
attach(nrDEG)
plot(logFC,-log10(Pvalue))
library(ggpubr)
df=nrDEG
df$v= -log10(Pvalue) #df新增加一列'v',值为-log10(P.Value)
ggscatter(df, x = "logFC", y = "v",size=0.5)

df$g=ifelse(df$Pvalue>0.05,'stable', #if 判断：如果这一基因的P.Value>0.01，则为stable基因
ifelse( df$logFC >1,'up', #接上句else 否则：接下来开始判断那些P.Value<0.01的基因，再if 判断：如果logFC >1.5,则为up（上调）基因
ifelse( df$logFC < -1,'down','stable') )#接上句else 否则：接下来开始判断那些logFC <1.5 的基因，再if 判断：如果logFC <1.5，则为down（下调）基因，否则为stable基因
)
table(df$g)
df$name=rownames(df)
head(df)
ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')
ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
label = "name", repel = T,
#label.select = rownames(df)[df$g != 'stable'] ,
label.select = c('ATF3'), #挑选一些基因在图中显示出来
palette = c("#00AFBB", "#E7B800", "#FC4E07") )
ggsave('volcano2.png')

}

for heatmap

if(F){
group_list
table(group_list)
x=deg$logFC #deg取logFC这列并将其重新赋值给x
names(x)=rownames(deg) #deg取probe_id这列，并将其作为名字给x
cg=c(names(head(sort(x),100)),#对x进行从小到大排列，取前100及后100，并取其对应的探针名，作为向量赋值给cg
names(tail(sort(x),100)))
library(pheatmap)
dat=genes_expr
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对dat按照cg取行，所得到的矩阵来画热图
n=t(scale(t(dat[cg,])))#通过“scale”对log-ratio数值进行归一化，现在的dat是行名为探针，列名为样本名，由于scale这个函数应用在不同组数据间存在差异时，需要行名为样本，因此需要用t(dat[cg,])来转换，最后再转换回来

n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n) #将ac的行名也就分组信息（是‘no TNBC’还是‘TNBC’）给到n的列名，即热图中位于上方的分组信息
pheatmap(n,show_colnames =T,
show_rownames = F,
cluster_cols = T,
annotation_col=ac,filename = 'heatmap_top200_DEG2.png') #列名注释信息为ac即分组信息
}

write.csv(DEG,file = 'DEG.csv')

rm(list = ls())
load('RIFanno_DEG.Rdata')
a=DEG
row.names(a)=DEG$ENTREZID
b=as.numeric(a$logFC)
names(b)=row.names(a)
head(b)
library('pathview')
p <- pathview(gene.data = b,
pathway.id = "04216",
species = "hsa",
out.suffix = "RIFCTL",
#limit = list(gene=as.integer(max(abs(a))))
)
p <- pathview(gene.data = b,
pathway.id = "04110",
species = "hsa",
out.suffix = "RIFCTL2",
limit = 2
#limit = list(gene=as.integer(max(abs(a))))
)
p <- pathview(gene.data = b,
pathway.id = "04350",
species = "hsa",
out.suffix = "RIFCTL1",
limit = 1
#limit = list(gene=as.integer(max(abs(a))))
)

check_diff_genes <- function(gene,genes_expr,group_list ){
if(length(gene)==1){
# gene='NKILA'
if(! gene %in% rownames(genes_expr)){
stop(paste0(gene,' in not in your expression matrix'))
}
df=data.frame(value=as.numeric(genes_expr[gene,]),
group=group_list)
library(ggpubr)
ggboxplot(df, "group", "value",
color = "group", palette =c("#00AFBB", "#E7B800"),
add = "jitter", shape = "group")
}else{
library(pheatmap)
cg=gene
cg=cg[cg %in% rownames(genes_expr) ]
warning(paste0('Only ',length(cg),' in ',length(gene),' genes are in your expression matrix'))
if(length(cg)<1){
stop('None of the gene in your expression matrix')
}
n=t(scale(t(genes_expr[cg,])))
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
ac=data.frame(group_list=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,show_rownames = T,cluster_cols= F,
annotation_col=ac)
}

}