#install.packages('boot')
library(boot)
###########寻找预后相关的G蛋白偶联受体##############
library(survival)      #引用包
pFilter= 0.05         #显著性过滤条件

library(limma)               #引用包
setwd("G:\\肾细胞癌")
# 数据处理
## 载入TCGA胃癌表达谱
load('KIRC_tpm.Rdata')
gpr=read.table('G_protein.txt')
gpr=as.data.frame(t(gpr))
gpr=gpr$V1
# 获取GPR基因
gpr=gpr[3:872]

# GPR表达谱的差异分析
data=exprSet_tcga_mRNA[rownames(exprSet_tcga_mRNA) %in% gpr,]
## 肿瘤和正常样品
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group_list=ifelse(group=="0",'tumor','normal')
group_list=factor(group_list,levels = c('normal','tumor'))
library(limma)
design=model.matrix(~ group_list)

fit=lmFit(data,design)
fit=eBayes(fit) 
allDiff=topTable(fit,adjust='fdr',coef=2,number=Inf,p.value=0.05) 


write.csv(allDiff,file ='allDiff.csv',quote = F)


## 更新data
data=data[rownames(allDiff),]
#删掉正常样品
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[,group==0]

## 去除在大多数样本中都表达为0的基因
keep <- rowSums(data>0) >= floor(0.75*ncol(data))
table(keep)
## 351
data<- data[keep,]

data=as.data.frame(t(data))


# 单因素寻找预后影响的GPR
## 读取生存数据
suv=read.table('TCGA-KIRC.survival.tsv',row.names = 1,header = T,check.names = F)
cli=dplyr::select(suv,'OS.time','OS')
colnames(cli)=c("futime", "fustat")

## 数据合并并输出结果
sameSample=intersect(row.names(data),row.names(cli))
data=data[sameSample,]
cli=cli[sameSample,]
#你需要牢记这样的结构，用于cox
out=cbind(cli,data)
out=cbind(id=row.names(out),out)
## 写出GPR的基础矩阵
write.table(out,file="expTime.txt",sep="\t",row.names=F,quote=F)

rt=read.table("expTime.txt", header=T, sep="\t", check.names=F, row.names=1)     #读取输入文件


##单因素cox分析
outTab=data.frame()
sigGenes=c("futime","fustat")
for(gene in colnames(rt[,3:ncol(rt)])){
  set.seed(123456)
  cox=coxph(Surv(futime, fustat) ~ rt[,gene], data = rt)
  coxSummary = summary(cox)
  coxP=coxSummary$coefficients[,"Pr(>|z|)"]
  if(coxP<pFilter){
    sigGenes=c(sigGenes,gene)
    outTab=rbind(outTab,
                 cbind(gene=gene,
                       HR=coxSummary$conf.int[,"exp(coef)"],
                       HR.95L=coxSummary$conf.int[,"lower .95"],
                       HR.95H=coxSummary$conf.int[,"upper .95"],
                       pvalue=coxP))
    print(coxP)
  }
}


#输出单因素结果
write.table(outTab,file="uniCox_gpr.txt",sep="\t",row.names=F,quote=F)
surSigExp=rt[,sigGenes]
surSigExp=cbind(id=row.names(surSigExp),surSigExp)
write.table(surSigExp,file="uniSigExp_gpr.txt",sep="\t",row.names=F,quote=F)


#### 热图
##########热图
#正常和肿瘤数目、
load('KIRC_tpm.Rdata')
gene=read.table('uniCox_gpr.txt',header = T)
gene=gene$gene
data=exprSet_tcga_mRNA[gene,]
group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2", "1", group)

conNum=length(group[group==1])       #正常组样品数目
treatNum=length(group[group==0])     #肿瘤组样品数目

sampleType=ifelse(group=='1',1,2)

identical(colnames(data),colnames(exprSet_tcga_mRNA))

#基因差异分析
sigVec=c()
allDiff=read.csv('allDiff.csv',header = T,row.names = 1)
alldiff_cox=allDiff[gene,]
pvalue=alldiff_cox$adj.P.Val
Sig=ifelse(pvalue<0.001,"***",ifelse(pvalue<0.01,"**",ifelse(pvalue<0.05,"*","")))
sigVec=paste0(gene, Sig)


## 另起一个compare矩阵，避免破坏原来的
compare=data
# 修饰一下行名
row.names(compare)=sigVec

#调整顺序，保证出图美观
normal=compare[,sampleType==1]
tumor=compare[,sampleType==2]
compare=cbind(normal,tumor)

#热图可视化
Type=c(rep("Normal",conNum), rep("Tumor",treatNum))
names(Type)=colnames(compare)
Type=as.data.frame(Type)
library(pheatmap)
a=pheatmap::pheatmap(compare,
                   annotation=Type,
                   breaks = c(seq(-3,3,length.out = 100)),
                   cluster_cols =F,
                   cluster_rows =T,
                   scale="row",
                   show_colnames=F,
                   show_rownames=T,
                   fontsize=6,
                   fontsize_row=7,
                   fontsize_col=6)

ggsave("pheatmap.png", plot = a, width = 15, height = 12, dpi = 500)


################先lasso,进行bootstrap_multicox回归#####################
gene=read.table('uniCox_gpr.txt',header = T)
gene=gene$gene
# 没有包先安装
#install.packages('survival')
library(survival)
library(survminer)
rt=read.table("expTime.txt", header=T, sep="\t", check.names=F, row.names=1)     #读取输入文件
rt=rt[,c('futime','fustat',gene)]




## 先lasso筛基因!!!!!!!
set.seed(2)   #设定随机种子 
x=as.matrix(rt[,c(3:ncol(rt))]) 
y=data.matrix(Surv(rt$futime,rt$fustat)) 

# 没包先安装
# install.packages('glmnet')
library(glmnet)
fit=glmnet(x, y, family = "cox", alpha = 1) 
plot(fit, xvar = "lambda", label = F)

cvfit = cv.glmnet(x, y, family="cox",nfolds = 10,alpha=1) 
plot(cvfit) 
#其中两条虚线分别指示了两个特殊的λ值 
abline(v = log(c(cvfit$lambda.min,cvfit$lambda.1se)),lty="dashed")

coef =coef(fit,s = cvfit$lambda.1se)
index = which(coef !=0)
actCoef = coef[index] 
lassoGene = row.names(coef)[index] 
geneCoef = cbind(Gene=lassoGene,Coef=actCoef) 
geneCoef   #查看模型的相关系数

## 剩下基因
gene=read.table('uniCox_gpr.txt',header = T)
gene=gene[gene$gene %in% geneCoef[,1],]

write.table(gene,file = 'uniCox_lasso_gpr.txt',quote = F,sep = '\t',row.names = F)



########进行bootstrap_multicox回归####
#install.packages('boot')
library(boot)
gene=read.table('uniCox_lasso_gpr.txt',header = T)
# boot_coef=coef/Boot_sd
gene=gene$gene
library(survival)
library(survminer)
rt=read.table("expTime.txt", header=T, sep="\t", check.names=F, row.names=1)     #读取输入文件
rt=rt[,c('futime','fustat',gene)]

# 初始HR
cox=coxph(Surv(futime, fustat) ~.,data = rt)
ggforest(cox)

#install.packages('boot')
library(boot)
rsq <- function(formula, data, indices) { 
  d <- data[indices,] 
  fit <- coxph(formula, data=d) 
  return(fit$coefficients) 
} 

## bootstrap，稍等待
set.seed(123456)
boot_results <- boot(data=rt, statistic=rsq, 
                     R=1000, formula=Surv(futime, fustat) ~ .)

## 单因素不良预后的这边可能变好，但是因为是多因素，不管
print(boot_results)

## 获取参数
coef=boot_results$t0
sd=as.data.frame(boot_results$t)
sd=apply(sd, 2, sd)

## 定义coef/sd为新的参数
ratio=coef/sd

GPR=data.frame('Coef'=coef,'boot_SD'=sd,'Coef\\/boot_SD'=ratio)
write.csv(GPR,file= 'GPR_coef.csv',quote = F)

# 构建GPRscore
load('KIRC_tpm.Rdata')
gene=read.table('uniCox_lasso_gpr.txt',header = T)
gene=gene$gene
data=exprSet_tcga_mRNA[gene,]

#删掉正常样品
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[,group==0]

# 读取系数
coef=read.csv('GPR_coef.csv',header = T,check.names = F)
coef=coef$Coef..boot_SD
gpr_score=c()
for (i in 1:ncol(data)) {
  score=sum(as.numeric(data[,i])*coef)
  gpr_score=c(gpr_score,score)
}

data=as.data.frame(t(data))
data$gpr_score=gpr_score

#读取生存数据
suv=read.table('TCGA-KIRC.survival.tsv',row.names = 1,header = T,check.names = F)
cli=dplyr::select(suv,'OS.time','OS')
colnames(cli)=c("futime", "fustat")

##数据合并并输出结果
sameSample=intersect(row.names(data),row.names(cli))
data=data[sameSample,]
cli=cli[sameSample,]


## K-M生存分析
rt=cbind(cli,data)
rt$futime=rt$futime/30

library(survival)
library(survminer)

### 中位值划分
Type=ifelse(data[,'gpr_score']<= median(rt$gpr_score), "Low", "High")
data=rt
data$group=Type
data$group=factor(data$group, levels=c("Low", "High"))
diff=survdiff(Surv(futime, fustat) ~ group, data = data)
length=length(levels(factor(data[,"group"])))
pValue=1-pchisq(diff$chisq, df=length-1)

fit <- survfit(Surv(futime, fustat) ~ group, data = data)

bioCol <- c("#FFA500", "#8B4513", "#1E90FF")
p=ggsurvplot(fit, 
             data=data,
             conf.int=F,
             pval=pValue,
             pval.size=6,
             legend.title='GPR_score',
             legend.labs=levels(factor(data[,"group"])),
             legend = c(0.88, 0.9),
             font.legend=12,
             xlab="Time(Months)",
             palette = bioCol,
             surv.median.line = "hv",
             risk.table=T,
             cumevents=F,
             risk.table.height=.25)

p    

save(data,file ='GPRscore_and_group.Rdata')
