LUAD文章复现代码汇总

最新推荐文章于 2024-07-08 00:01:13 发布
JOYE AO
最新推荐文章于 2024-07-08 00:01:13 发布
阅读量579
点赞数 1
文章标签： python 开发语言
本文链接：https://blog.csdn.net/zhuoyi520/article/details/132240961
版权
1+1#1+1代表普通的运算
#install.packages("tidyverse")
#也可以通过右下角的操作栏通过点击来安装R包
#chooseBioCmirror()
#a<-6
#a
#以上是试验赋值的用法，还可以用=来表示，不过==代表一种逻辑判断
#setwd("DAY1")#此代码用来设置工作目录（临时工作目录，退出后回到初始目录），以免在使用过程中将储存文件弄混
#dir.create()可以用来直接创建目录，不过麻烦
#getwd()#查看当前目录
#setwd("DAY1.1")
#ctrl+F可以用来查找代码框中的代码，ctrl+Z是撤销一次的含义
#read.table()#读取文本文档
#write.CSV()#是输出为表格
write.table(counts01A,"counts01A.txt",sep="\t",row.names=T,col.names=NA,quote=F)#是输出文本文档#对文件的读取或输出可以直接点击右上角的读取来导入#通常txt格式会更不容易出错
#t()是用于行列反转的函数
#当一个可以判断为数据框的格式但R语言并没有将其判断为数据框时，可以通过函数as.data.frame()来改变

#class()#用来判断数据框类型，若是数据框内的某一列数据通过输入$符号来确定
#[,]代表对数据框进行提取。其中：表示多少到多少。并且，前是行后是列
#%>%表示传导符号

#install.packages("tidyverse")
#library(tidyverse)
#a<-c("a","b","a","b","c")
#duplicated(a)#用来判断数据里的元素是否重复
#a<-!duplicated(a)#感叹号是指将数据集内判断结果反过来，[]中括号可以把判断中是TRUE的元素提取出来
#a[a]#以上操作中将判断结果赋予了a故导致原本数据的丢失
#a<-c("a","b","a","b","c")
#duplicated(a)
#a<-a[!duplicated(a)]
#a
#这才是正确操作结果

#inner_join()#是将两数据框中数据相同的进行合并的函数，包括（表一，表二，by='合并标准'）
#left_join()#则是代表以函数内左的数据框为标准进行合并


###接下来进行数据的下载实战#####
##建立好文件##
#setwd("TCGA-LUAD")
#setwd('TCGAdata')
library(tidyverse)
install.packages('BiocManager')
library(BiocManager)
BiocManager::install("BioinformaticsFMRP/TCGAbiolinksGUI.data")
BiocManager::install("remotes")
BiocManager::install("ExperimentHub")
BiocManager::install("BioinformaticsFMRP/TCGAbiolinks")
#library(TCGAbiolinks)
cancer_type<-"TCGA-LUAD"
expquery<-GDCquery(project=cancer_type,
                   data.category="Transcriptome Profiling",
                   data.type="Gene Expression Quantification",
                   workflow.type="STAR - Counts")
GDCdownload(expquery,directory="GDCdata")
expquery2<-GDCprepare(expquery,directory="GDCdata",summarizedExperiment = T)
save(expquery2,file="luad.gdc_2023.rad")#对数据进行保存，rad格式是属于R的格式

setwd('TCGA-LUAD')
setwd("TCGAdata")
load("luad.gdc_2023.rad")
load("gene_annotation_2023.rad")
gene_annotation_2023<-data.txt


library(tidyverse)
counts<-expquery2@assays@data@listData[["unstranded"]]
colnames(counts)<-expquery2@colData@rownames
rownames(counts)<-expquery2@rowRanges@ranges@NAMES
counts<-as.data.frame(counts)
counts<-rownames_to_column(counts,var = "ENSEMBL")
counts<-inner_join(counts,gene_annotation_2023,"ENSEMBL")
counts<-counts[!duplicated(counts$symbol),]


rownames(counts)<-NULL#去除行名
counts<-column_to_rownames(counts,var = "symbol")
table(counts$type)#数一下有多少的type
counts<-counts[counts$type == "protein_coding",]
counts<-counts[,-c(1,ncol(counts))]
colnames(counts)<-substring(colnames(counts),1,16)#将列名中1~16字符提取出来 
counts<-counts[,!duplicated(colnames(counts))]



table(substring(colnames(counts),14,16))#通常01A代表肿瘤样本，11A代表正常样本
counts01A<-counts[,substring(colnames(counts),14,16)=="01A"]
counts11A<-counts[,substring(colnames(counts),14,16)=="11A"]#将样本中01A和11A的样本提出来做分析



write.table(counts01A,"counts01A.txt",sep="\t",row.names=T,col.names=NA,quote=F)
write.table(counts11A,"counts11A.txt",sep="\t",row.names=T,col.names=NA,quote=F)
###接下来提取tpms###前面提的counts是用来做差异分析的，而tpms是用来做后面分析的

tpms1<-expquery2@assays@data@listData[["tpm_unstrand"]]
colnames(tpms1)<-expquery2@colData@rownames
rownames(tpms1)<-expquery2@rowRanges@ranges@NAMES
tpms1<-tpms1%>%
  as.data.frame()%>%
  rownames_to_column("ENSEMBL")%>%
  inner_join(gene_annotation_2023,"ENSEMBL")%>%
  .[!duplicated(.$symbol),]
rownames(tpms1)<-NULL
tpms1<-tpms1 %>%column_to_rownames("symbol")
  
table(tpms1$type)#查看分类计数

tpms1<-tpms1[tpms1$type == "protein_coding",]
tpms1<-tpms1[,-c(1,ncol(tpms1))]
colnames(tpms1)<-substring(colnames(tpms1),1,16)#将列名中1~16字符提取出来 
tpms1<-tpms1[,!duplicated(colnames(tpms1))]
table(substring(colnames(tpms1),14,16))

tpms01A<-tpms1[,substring(colnames(tpms1),14,16)=="01A"]
tpms11A<-tpms1[,substring(colnames(tpms1),14,16)=="11A"]


identical(rownames(counts01A),rownames(tpms01A))
identical(rownames(counts11A),rownames(tpms11A))
identical(colnames(counts01A),colnames(tpms01A))
identical(colnames(counts11A),colnames(tpms11A))#验证几个样本是否一致

write.table(tpms01A,"tpms01A.txt",sep="\t",row.names=T,col.names=NA,quote=F)
write.table(tpms11A,"tpms11A.txt",sep="\t",row.names=T,col.names=NA,quote=F)


counts<-cbind(counts01A,counts11A)#cbind是以列colnames来合并，rbind是以行rownames来合并
tpms<-cbind(tpms01A,tpms11A)
write.table(tpms,"tpms.txt",sep="\t",row.names=T,col.names=NA,quote=F)
write.table(counts,"counts.txt",sep="\t",row.names=T,col.names=NA,quote=F)


range(tpms)#查看数据范围
range(tpms01A)
range(tpms11A)
tpms_log2<-log2(tpms+1)
range(tpms_log2)
tpms01A_log2<-log2(tpms01A+1)
range(tpms01A_log2)
tpms11A_log2<-log2(tpms11A+1)
range(tpms11A_log2)
write.table(tpms_log2,"tpms_log2.txt",sep="\t",row.names=T,col.names=NA,quote=F)
write.table(tpms01A_log2,"tpms01A_log2.txt",sep="\t",row.names=T,col.names=NA,quote=F)
write.table(tpms11A_log2,"tpms11A_log2.txt",sep="\t",row.names=T,col.names=NA,quote=F)




###计算肿瘤患者的免疫评分#####
#（基质，免疫，肿瘤）肿瘤纯度=肿瘤所占比例
setwd("TCGA-LUAD")
setwd("ESTIMATE")
#rfoge<-"http://r-forge.r-project.org"#RFORGE就是R语言自己的R包平台
#install.packages("estimate",repos=rfoge,dependencies=TRUE)
library(estimate)#estimate是用来计算免疫评分的
library(tidyverse)
exp<-read.table("tpms01A_log2.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                ,header=T)

filterCommonGenes(input.f="tpms01A_log2.txt",
                  output.f="tpms01A_log2.gct",
                  id="GeneSymbol")
estimateScore("tpms01A_log2.gct",
              "tpms01A_log2_estimate_score.txt",
              platform="affymetrix")#affymetrix是一个基因芯片平台，相比于其它平台结果更好

ESTIMATE_result<-read.table("tpms01A_log2_estimate_score.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                            ,header=T)
ESTIMATE_result<-ESTIMATE_result[,-1]
column_to_rownames(ESTIMATE_result,var="NAME")
colnames(ESTIMATE_result)<-ESTIMATE_result[1,]
ESTIMATE_result<-t(ESTIMATE_result[-1,])#通常行列转换后会变得不是数据框形式
ESTIMATE_result2<-as.data.frame(ESTIMATE_result)#可以直接用ESTIMATE_result<-as.data.frame(t(ESTIMATE_result[-1,]))


#dat1 <- as.data.frame(ESTIMATE_result2)
#dat1 <-rownames_to_column(dat1 )
#dat1 <- str_replace(ESTIMATE_result2[,1], ".", "_") #错误代码，练习更换字符

rownames(ESTIMATE_result)<-colnames(exp)
write.table(ESTIMATE_result,"ESTIMATE_result.txt",sep="\t",row.names=T,col.names=NA,quote=F)


##进行生存分析####
setwd("Survival_data")
library(tidyverse)
survival_data<-read.table("OS.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
           ,header=T)
survival_data<-survival_data[,2:3]
survival<-rownames_to_column(survival_data,"sample")
survival$name<-paste0(survival$sample,"A")#当粘贴到的数据框中没有name这列时，R会自动生成新的一列
#paste0表示粘贴连接一字符“”
table(substring(survival$name,14,16))
rownames(survival)<-survival$name
survival<-survival[,2:3]

#接下将生存信息与基因表达谱合并起来#
tpms01A_log2<-read.table("tpms01A_log2.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                         ,header=T)
a<-intersect(colnames(tpms01A_log2),rownames(survival))#对生存数据和基因表达谱
#~的sample名取交集，主要是为了得到一个额外的sample顺序，用来接下来将两者sample对齐
table(substring(a,14,16))
exp_01A<-tpms01A_log2[,a]
surv_01A<-survival[a,]#以顺序或以集合a作为模板来将两个数据中以行或列进行数据提取
exp_01A<-exp_01A%>%t()%>%as.data.frame()
identical(rownames(exp_01A),rownames(surv_01A))
exp_surv_01A<-cbind(surv_01A,exp_01A)#cbind是将两数据列和并起来（前提是行相同）
#~rbind是将两数据行合并起来（前提是列相同）
#将合并后的生存信息和基因表达谱合并的文件保存
write.table(exp_surv_01A,"exp_surv_01A.txt",sep="\t",row.names=T,col.names=NA,quote=F)


#合并生存信息与免疫评分数据（ESTIMATE）####
ESTIMATE_result<-read.table("ESTIMATE_result.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                            ,header=T)
identical(rownames(ESTIMATE_result),rownames(surv_01A))
ESTIMATE_result_surv_01A<-cbind(surv_01A,ESTIMATE_result)
write.table(ESTIMATE_result_surv_01A,"ESTIMATE_result_surv_01A.txt",sep="\t",row.names=T,col.names=NA,quote=F)


setwd("TCGA-LUAD")
setwd("survival")
surv01A<-read.table("ESTIMATE_result_surv_01A.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                    ,header=T)
surv01A$OS.time<-surv01A$OS.time/365#将生存时间从天变为年，利于后续分析





surv01A$group<-ifelse(surv01A$ImmuneScore>median(surv01A$ImmuneScore),"High","Low")
#该函数添加的high和low都是字符形式的，而非因子形式，故无法做二分类变量判断
#此次的评分分组为关键，分组不同所做的生存分析标准不同


class(surv01A$group)
surv01A$group<-factor(surv01A$group,levels=c("Low","High"))
class(surv01A$group)
table(surv01A$group)

#install.packages("survival")#安装做生存分析的包
library(survival)
fitd<-survdiff(Surv(OS.time,OS)~group,
               data=surv01A,
               na.action = na.exclude)#na.omit：这个操作会从数据中删除任何包含缺失值的行。
#~na.fail：这个操作会在数据中存在缺失值时导致模型失败。
#~na.exclude：这个操作会将包含缺失值的行从分析中排除，但不从数据中删除。
pValue<-1-pchisq(fitd$chisq,length(fitd$n)-1)#该步骤是用来做卡方检验P值
#~与统计学上卡方检验一样，小于0.01才具有统计学意义


#拟合生存曲线
fit<-survfit(Surv(OS.time,OS)~group,data=surv01A)
summary(fit)
p.lab<-paste0("P",ifelse(pValue<0.001,"<0.001",paste0("=",round(pValue,3))))
#该步只是将到的P值规范地写下来
#install.packages("survminer")
library(survminer)
ggsurvplot(
  fit,
  data = surv01A,
  palette = "jco",#配色采用Jco，lancet柳叶刀配色，jama是jama期刊的配色
  pval = p.lab,
  conf.int = TRUE,#显示置信区间
  risk.table = TRUE,#显示风险表
  risk.table.col="strata",
  legend.labs=c("Low","High"),#图例
  size=1,
  xlim=c(0,20),#x轴长度，根据自己数据中OS.time的时间来决定，此数据中最大为19年多，
  #故用20作X轴长度，可用range函数判断自己数据范围
  break.time.by=5,#X轴步长为5
  legend.title="ImmuneScore",#TITLE
  surv.median.line = "hv",#限制垂直和水平的中位生存
  ylab="Survival probability(%)",#修改Y轴标签
  xlab="Time(Years)",#修改x轴标签
  ncensor.plot=TRUE,#显示缺失图块
  ncensor.plot.height=0.25,
  risk.table.y.text=FALSE)
dev.off()#不知道干什么用暂时


#StromalScore#只需要将上述分组时的High和Low通过机制评分来进行分组便可以
#~进行机制的生存分析



surv01A$group<-ifelse(surv01A$StromalScore>median(surv01A$StromalScore),"High","Low")
#该函数添加的high和low都是字符形式的，而非因子形式，故无法做二分类变量判断
class(surv01A$group)
surv01A$group<-factor(surv01A$group,levels=c("Low","High"))
class(surv01A$group)
table(surv01A$group)

#install.packages("survival")#安装做生存分析的包
library(survival)
fitd<-survdiff(Surv(OS.time,OS)~group,
               data=surv01A,
               na.action = na.exclude)#na.omit：这个操作会从数据中删除任何包含缺失值的行。
#~na.fail：这个操作会在数据中存在缺失值时导致模型失败。
#~na.exclude：这个操作会将包含缺失值的行从分析中排除，但不从数据中删除。
pValue<-1-pchisq(fitd$chisq,length(fitd$n)-1)#该步骤是用来做卡方检验P值
#~与统计学上卡方检验一样，小于0.01才具有统计学意义


#拟合生存曲线
fit<-survfit(Surv(OS.time,OS)~group,data=surv01A)
summary(fit)
p.lab<-paste0("P",ifelse(pValue<0.001,"<0.001",paste0("=",round(pValue,3))))
#该步只是将到的P值规范地写下来
#install.packages("survminer")
library(survminer)
ggsurvplot(
  fit,
  data = surv01A,
  palette = "jco",#配色采用Jco，lancet柳叶刀配色，jama是jama期刊的配色
  pval = p.lab,
  conf.int = TRUE,#显示置信区间
  risk.table = TRUE,#显示风险表
  risk.table.col="strata",
  legend.labs=c("Low","High"),#图例
  size=1,
  xlim=c(0,20),#x轴长度，根据自己数据中OS.time的时间来决定，此数据中最大为19年多，
  #故用20作X轴长度，可用range函数判断自己数据范围
  break.time.by=5,#X轴步长为5
  legend.title="StromalScore",#TITLE
  surv.median.line = "hv",#限制垂直和水平的中位生存
  ylab="Survival probability(%)",#修改Y轴标签
  xlab="Time(Years)",#修改x轴标签
  ncensor.plot=TRUE,#显示缺失图块
  ncensor.plot.height=0.25,
  risk.table.y.text=FALSE)




#ESTIMATEScore



surv01A$group<-ifelse(surv01A$ESTIMATEScore>median(surv01A$ESTIMATEScore),"High","Low")
#该函数添加的high和low都是字符形式的，而非因子形式，故无法做二分类变量判断
class(surv01A$group)
surv01A$group<-factor(surv01A$group,levels=c("Low","High"))
class(surv01A$group)
table(surv01A$group)

#install.packages("survival")#安装做生存分析的包
library(survival)
fitd<-survdiff(Surv(OS.time,OS)~group,
               data=surv01A,
               na.action = na.exclude)#na.omit：这个操作会从数据中删除任何包含缺失值的行。
#~na.fail：这个操作会在数据中存在缺失值时导致模型失败。
#~na.exclude：这个操作会将包含缺失值的行从分析中排除，但不从数据中删除。
pValue<-1-pchisq(fitd$chisq,length(fitd$n)-1)#该步骤是用来做卡方检验P值
#~与统计学上卡方检验一样，小于0.01才具有统计学意义


#拟合生存曲线
fit<-survfit(Surv(OS.time,OS)~group,data=surv01A)
summary(fit)
p.lab<-paste0("P",ifelse(pValue<0.001,"<0.001",paste0("=",round(pValue,3))))
#该步只是将到的P值规范地写下来
#install.packages("survminer")
library(survminer)
ggsurvplot(
  fit,
  data = surv01A,
  palette = "jco",#配色采用Jco，lancet柳叶刀配色，jama是jama期刊的配色
  pval = p.lab,
  conf.int = TRUE,#显示置信区间
  risk.table = TRUE,#显示风险表
  risk.table.col="strata",
  legend.labs=c("Low","High"),#图例
  size=1,
  xlim=c(0,20),#x轴长度，根据自己数据中OS.time的时间来决定，此数据中最大为19年多，
  #故用20作X轴长度，可用range函数判断自己数据范围
  break.time.by=5,#X轴步长为5
  legend.title="ESTIMATEScore",#TITLE
  surv.median.line = "hv",#限制垂直和水平的中位生存
  ylab="Survival probability(%)",#修改Y轴标签
  xlab="Time(Years)",#修改x轴标签
  ncensor.plot=TRUE,#显示缺失图块
  ncensor.plot.height=0.25,
  risk.table.y.text=FALSE)





#分组画图，肿瘤阶段，转移，####


setwd("TCGA-LUAD")
setwd("Clinical")
library(tidyverse)
load("luad.gdc_2023.rad")
clinical<-as.data.frame(expquery2@colData)%>%.[!duplicated(.$sample),]
#提取文章中所需要的临床信息，性别，年龄。。。。
clinical<-clinical[,c("gender","age_at_index","ajcc_pathologic_stage",
                      "ajcc_pathologic_t","ajcc_pathologic_n","ajcc_pathologic_m")]
#判断一下每一列的数据类型，做到心中有数
class(clinical$gender)
class(clinical$age_at_index)#结果是"integer"代表整数
class(clinical$ajcc_pathologic_stage)
class(clinical$ajcc_pathologic_t)
class(clinical$ajcc_pathologic_n)
class(clinical$ajcc_pathologic_m)

#分组计数一下，可以看出结果
table(clinical$gender)
table(clinical$age_at_index)
table(clinical$ajcc_pathologic_stage)#只需要里面的一二三期便可以不需要分压型，故要去掉亚型分类
table(clinical$ajcc_pathologic_t)
table(clinical$ajcc_pathologic_n)
table(clinical$ajcc_pathologic_m)

#去掉亚型分类
clinical$ajcc_pathologic_stage<-gsub("A","",clinical$ajcc_pathologic_stage)
clinical$ajcc_pathologic_stage<-gsub("B","",clinical$ajcc_pathologic_stage)
clinical$ajcc_pathologic_t<-gsub("a","",clinical$ajcc_pathologic_t)
clinical$ajcc_pathologic_t<-gsub("b","",clinical$ajcc_pathologic_t)
clinical$ajcc_pathologic_m<-gsub("a","",clinical$ajcc_pathologic_m)
clinical$ajcc_pathologic_m<-gsub("b","",clinical$ajcc_pathologic_m)
#数据里的X表示未知并不是亚型
rownames(clinical)<-substring(rownames(clinical),1,16)

#接下来要将基因表达谱导入进来并且将与临床数据合并#和先前的生存数据一样
exp01A<-read.table("tpms01A_log2.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                   ,header=T)
clinical01A<-clinical[colnames(exp01A),]#以exp01A中的列名为模板对clinical中数据进行提取
exp01A<-exp01A%>%t()%>%as.data.frame()
clinical.expr01A<-cbind(clinical01A,exp01A)
write.table(clinical.expr01A,"clinical.expr01A.txt",sep="\t",row.names=T,col.names=NA,quote=F)


#再将临床数据与免疫评分进行合并####
ESTIMATE_result<-read.table("ESTIMATE_result.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                            ,header=T)
clinical.ESTIMATE_result01A<-cbind(clinical01A,ESTIMATE_result)
write.csv(clinical.ESTIMATE_result01A,file="clinical.ESTIMATE_result01A.csv")




#自己尝试用R作分组箱式图####
setwd("TCGA-LUAD")
setwd("Clinical")
df<-read.csv(file="clinical.ESTIMATE_result01A.csv")
library(tidyverse)
df1<-df[complete.cases(df$ajcc_pathologic_stage),]
table(df1$ajcc_pathologic_stage)
#df4$ajcc_pathologic_m<-gsub("MX",NA,df4$ajcc_pathologic_m)#将TX值换位缺失值
#df1<-df1[complete.cases(df3$ajcc_pathologic_n),]#过滤掉了stage中的缺失值
df1<- df1[order(df1$ajcc_pathologic_stage), ]#使数据按顺序排列，以便于画图
library(ggplot2)
library(RColorBrewer)#这是R语言的调色板，可以选择绘图颜色
display.brewer.all()#查看调色板
mycol<-brewer.pal(n=2,name="Set2")#设置了四个颜色，用Set2为色调
mycol#查看我的调色板的代码名字
mycol<-c("#66C2A5" ,"#FC8D62", "#8DA0CB", "#E78AC3")#设置为集合
library(ggpubr)
stage_ESTIMATEScore<-ggboxplot( df1,
                              x="ajcc_pathologic_stage",
                              y="ESTIMATEScore",
                              #color="ajcc_pathologic_stage",#外边框颜色
                              palette=mycol,
                              legend="none",
                              #add="jitter",#添加样本点
                              outlier.shape=NA,
                              fill="ajcc_pathologic_stage",#颜色填充
                              title="stage",
                              xlab="stage",
                              ylab="ESTIMATEScore",
           size=0.5,)+stat_boxplot(geom="errorbar",width=0.5,size=0.1)+
  geom_boxplot(fill=mycol,size=0.1)+#将中轴线覆盖
  stat_compare_means(method="kruskal.test",label="p.signif")+
  stat_compare_means(comparisons=list(c("Stage I","Stage IV")),method="wilcox",label="p.format")
stage_ESTIMATEScore

#小提琴图
stage_ImmuneScore2<-ggviolin( df1,
                              x="ajcc_pathologic_stage",
                              y="ImmuneScore",
                              #color="ajcc_pathologic_stage",#外边框颜色
                              palette=mycol,
                              legend="none",
                              #add="jitter",#添加样本点
                              add="boxplot",#添加箱式图
                              alpha=0.5,#调整透明度
                              outlier.shape=NA,
                              fill="ajcc_pathologic_stage",#颜色填充
                              title="stage_ImmuneScore",
                              xlab="stage",
                              ylab="ImmuneScore",
                              size=0.1,)+stat_boxplot(geom="errorbar",width=0.3,size=0.1)+
  geom_boxplot(fill=mycol,size=0.1,width=0.3)+#将中轴线覆盖
  stat_compare_means(method="kruskal.test",label="p.signif")+
  stat_compare_means(comparisons=list(c("Stage I","Stage IV")),method="wilcox",label="p.format")
stage_ImmuneScore2



####差异分析####
setwd("TCGA-LUAD")
setwd("Immune_DEG")
library("BiocManager")
BiocManager::install("DESeq2")#DESeq2主要包括以下功能：
##数据归一化：DESeq2提供了多种数据归一化方法，如TMM（Trimmed Mean of M-values）归一化。
##差异分析：DESeq2使用负二项分布模型及贝叶斯估计方法，基于基因或转录本的计数数据，计算基因/转录本在不同条件之间的差异表达。差异分析结果包括差异表达基因/转录本的统计显著性和折变倍数等信息。
##数据可视化：DESeq2提供了丰富的数据可视化功能，包括MA图、PCA图、热图等，帮助用户更好地理解数据的特征和差异表达结果
library(DESeq2)
library(tidyverse)
counts_01A<-read.table("counts01A.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                       ,header=T)
###TCGA做差异分析时是使用的counts数据中肿瘤组
estimate<-read.table("ESTIMATE_result.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                     ,header=T)
x<-"ImmuneScore"
med<-as.numeric(median(estimate[,x]))
counts_01A2<-as.data.frame(t(counts_01A))
estimate2<-as.data.frame(t(estimate))
identical(rownames(counts_01A2),rownames(estimate))

conditions=data.frame(sample=rownames(counts_01A2),group=factor(ifelse(estimate[,x]>med,"high","low"),levels=c("low","high")))
conditions<-column_to_rownames(conditions,"sample")   

#接下来是差异分析的准备工作
dds<-DESeqDataSetFromMatrix(
  countData =counts_01A,#计数矩阵，列名为ID必须与colData数据行名一致
  colData = conditions,#要获取的数据的内容，即分组文件(如处理组和对照组)
  desig=~group)#以什么来将所获取的基因分组

#开始差异分析
dds<-DESeq(dds)
#DESeq()内部做了统计分析，我们得到差异基因。所以这也是为什么要用的计数矩阵是原始的计数矩阵，未被标准化处理过
resultsNames(dds)#查看一下结果是high比上low还是反之
res<-results(dds)
save(res,file="DEG_ImmuneScore.rda")                    


####绘制热图####
DEG<-as.data.frame(res)##结果中的log2Foldchange的结果就表示log2(某个基因在免疫评分高的样本的表达量/该基因在免疫评分低的样本中的表达量)
#读取基因表达谱
exp<-read.table("tpms01A_log2.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
           ,header=T)#Counts数据是原始的读取次数数据，而TPMs数据是通过归一化和标准化Count数据得到的相对表达量数据。TPMs数据更适合用于不同样本之间的比较和基因表达分析。
#添加上下调信息，因为我们需要真正具有高表达或低表达的才作为上调或下调，也就是要除去那些表达差异不明显的基因，
#~通常文章以log2(high/low)大于1作UP组，小于-1作DOWN组。即至少差异有2.7倍
a<-1
b<--1
type1<-(DEG$padj<0.05)&(DEG$log2FoldChange<b )
type2<-(DEG$padj<0.05)&(DEG$log2FoldChange>a)
DEG$change=ifelse(type1,"DOWN",ifelse(type2,"UP","NOT"))
table(DEG$change)

#install.packages("pheatmap")#画的图的R包
library(pheatmap)
DEG$change<-gsub("NOT",NA,DEG$change)#将我不需要的数据变为空值
table(DEG$change)
DEG<-DEG[complete.cases(DEG$change),]#将DEG中完整的行提取出来
DEG<-DEG[order(DEG$change),]#进行排序,实际不需要排序#以上步骤也可以用filter()函数筛选来做
d<-rownames(DEG)
exp_diff<-exp[d,]#这就是差异表达的基因
groups<-conditions#分组信息
a<-filter(groups,group=="high")
b<-filter(groups,group=="low")
exp_diff_high<-exp_diff[,rownames(a)]
exp_diff_low<-exp_diff[,rownames(b)]
exp_diff<-cbind(exp_diff_high,exp_diff_low)#将样本对称分为左右，左边样本为免疫评分高，右边为低
#开始画图
pheatmap(exp_diff,
         annotation_col=groups,
         scale="row",
         show_rownames=F,
         shou_colnames=F,
        color=colorRampPalette(c("#000080FF", "white", "#CD262680"))(100),# #热图色块颜色是从蓝到红分为100个等级
         cluster_cols=F,
         cluster_rows=T,
         fontsize=10,
         fontsize_row=3,
         fontsize_col=3)
original_color <- "#000080FF"  # 原始颜色为红色
adjusted_color <- adjustcolor(original_color,  1)  # 将颜色亮度调暗一半
adjusted_color






####差异分析基质评分####
setwd("TCGA-LUAD")
setwd("Stromal_DEG")
library("BiocManager")
#BiocManager::install("DESeq2")#DESeq2主要包括以下功能：
##数据归一化：DESeq2提供了多种数据归一化方法，如TMM（Trimmed Mean of M-values）归一化。
##差异分析：DESeq2使用负二项分布模型及贝叶斯估计方法，基于基因或转录本的计数数据，计算基因/转录本在不同条件之间的差异表达。差异分析结果包括差异表达基因/转录本的统计显著性和折变倍数等信息。
##数据可视化：DESeq2提供了丰富的数据可视化功能，包括MA图、PCA图、热图等，帮助用户更好地理解数据的特征和差异表达结果
library(DESeq2)
library(tidyverse)
counts_01A<-read.table("counts01A.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                       ,header=T)
###TCGA做差异分析时是使用的counts数据中肿瘤组
estimate<-read.table("ESTIMATE_result.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                     ,header=T)
x<-"StromalScore"
med<-as.numeric(median(estimate[,x]))
counts_01A2<-as.data.frame(t(counts_01A))
estimate2<-as.data.frame(t(estimate))
identical(rownames(counts_01A2),rownames(estimate))

conditions=data.frame(sample=rownames(counts_01A2),group=factor(ifelse(estimate[,x]>med,"high","low"),levels=c("low","high")))
conditions<-column_to_rownames(conditions,"sample")   

#接下来是差异分析的准备工作
dds<-DESeqDataSetFromMatrix(
  countData =counts_01A,#计数矩阵，列名为ID必须与colData数据行名一致
  colData = conditions,#要获取的数据的内容，即分组文件(如处理组和对照组)
  desig=~group)#以什么来将所获取的基因分组

#开始差异分析
dds<-DESeq(dds)
#DESeq()内部做了统计分析，我们得到差异基因。所以这也是为什么要用的计数矩阵是原始的计数矩阵，未被标准化处理过
resultsNames(dds)#查看一下结果是high比上low还是反之
res<-results(dds)
save(res,file="DEG_StromalScore.rda")                    


####绘制热图####
DEG<-as.data.frame(res)##结果中的log2Foldchange的结果就表示log2(某个基因在免疫评分高的样本的表达量/该基因在免疫评分低的样本中的表达量)
#读取基因表达谱
exp<-read.table("tpms01A_log2.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                ,header=T)#Counts数据是原始的读取次数数据，而TPMs数据是通过归一化和标准化Count数据得到的相对表达量数据。TPMs数据更适合用于不同样本之间的比较和基因表达分析。
#添加上下调信息，因为我们需要真正具有高表达或低表达的才作为上调或下调，也就是要除去那些表达差异不明显的基因，
#~通常文章以log2(high/low)大于1作UP组，小于-1作DOWN组。即至少差异有2.7倍
a<-1
b<--1
type1<-(DEG$padj<0.05)&(DEG$log2FoldChange<b )
type2<-(DEG$padj<0.05)&(DEG$log2FoldChange>a)
DEG$change=ifelse(type1,"DOWN",ifelse(type2,"UP","NOT"))
table(DEG$change)

#install.packages("pheatmap")#画的图的R包
library(pheatmap)
DEG$change<-gsub("NOT",NA,DEG$change)#将我不需要的数据变为空值
table(DEG$change)
DEG<-DEG[complete.cases(DEG$change),]#将DEG中完整的行提取出来
DEG<-DEG[order(DEG$change),]#进行排序,实际不需要排序#以上步骤也可以用filter()函数筛选来做
d<-rownames(DEG)
exp_diff<-exp[d,]#这就是差异表达的基因
groups<-conditions#分组信息
a<-filter(groups,group=="high")
b<-filter(groups,group=="low")
exp_diff_high<-exp_diff[,rownames(a)]
exp_diff_low<-exp_diff[,rownames(b)]
exp_diff<-cbind(exp_diff_high,exp_diff_low)#将样本对称分为左右，左边样本为免疫评分高，右边为低
#开始画图
pheatmap(exp_diff,
         annotation_col=groups,
         scale="row",
         show_rownames=F,
         shou_colnames=F,
         color=colorRampPalette(c("#000080FF", "white", "#CD262680"))(100),# #热图色块颜色是从蓝到红分为100个等级
         cluster_cols=F,
         cluster_rows=T,
         fontsize=10,
         fontsize_row=3,
         fontsize_col=3)
original_color <- "#000080FF"  # 原始颜色为红色
adjusted_color <- adjustcolor(original_color,  1)  # 将颜色亮度调暗一半
adjusted_color



####将两次差异分析的结果取交集####
setwd("TCGA-LUAD")
setwd("Immune_Stromal_DEG")
#提取免疫的上下调基因
load("DEG_ImmuneScore.rda")
DEG<-as.data.frame(res)
#添加上下调基因#和差异分析步骤一样
a<-1
b<--1
type1<-(DEG$padj<0.05)&(DEG$log2FoldChange<b )
type2<-(DEG$padj<0.05)&(DEG$log2FoldChange>a)
DEG$change=ifelse(type1,"DOWN",ifelse(type2,"UP","NOT"))
table(DEG$change)
#用filter()函数来删选提取出上下调基因的信息
library(tidyverse)
e<-filter(DEG,change=="UP")
f<-filter(DEG,change=="DOWN")
write.csv(e,file="Immune_up.csv")
write.csv(f,file="Immune_down.csv")

#基质的上下调基因也以同样的方式提取出来
load("DEG_StromalScore.rda")
DEG<-as.data.frame(res)
#添加上下调基因#和差异分析步骤一样
a<-1
b<--1
type1<-(DEG$padj<0.05)&(DEG$log2FoldChange<b )
type2<-(DEG$padj<0.05)&(DEG$log2FoldChange>a)
DEG$change=ifelse(type1,"DOWN",ifelse(type2,"UP","NOT"))
table(DEG$change)
#用filter()函数来删选提取出上下调基因的信息
library(tidyverse)
j<-filter(DEG,change=="UP")
h<-filter(DEG,change=="DOWN")
write.csv(j,file="Stromal_up.csv")
write.csv(h,file="Stromal_down.csv")


##用交集差异基因做(GO  KEGG)富集分析
setwd("TCGA-LUAD")
setwd("FUJI_Immune_Stromal_DEG")
library(tidyverse)
library(BiocManager)
#安装需要的包
BiocManager::install('clusterProfiler')
BiocManager::install('org.Hs.eg.db')
library(org.Hs.eg.db)
library(clusterProfiler)
##导入交集基因
#导入immune或stromal差异分析的结果，主要是为了获取各基因信息，对后面的评分高低不在乎
load("DEG_ImmuneScore.rda")
DEG<-as.data.frame(res)
DEG<-DEG[DEG_final$SYMBOL,]
DEG<-rownames_to_column(DEG,"SYMBOL")
genelist<-bitr(DEG$SYMBOL,fromType = "SYMBOL",
               toType = "ENTREZID",OrgDb = 'org.Hs.eg.db')#对基因名进行了ID转换
DEG<-inner_join(DEG,genelist,by="SYMBOL")


##GO分析####
ego<-enrichGO(gene =DEG$ENTREZID,
              OrgDb = org.Hs.eg.db,
              ont = "all",#包括，BP，MF ,CC
              pAdjustMethod = 'BH',#矫正方法
              minGSSize = 1,
              pvalueCutoff = 0.05,
              qvalueCutoff = 0.05,
              readable = TRUE)
ego_res<-ego@result
save(ego,ego_res,file = "GO_DEG_final.Rda")

####KEGG####
kk<-enrichKEGG(gene = DEG$ENTREZID,
               organism = "hsa",
               keyType = "kegg",
               pvalueCutoff = 0.1,
               qvalueCutoff = 0.1,)
library(DOSE)
kk<-DOSE::setReadable(kk,
                          OrgDb = org.Hs.eg.db,
                          keyType = "ENTREZID")
head(kk)
kk_res<-kk@result
save(kk,kk_res,file="KEGG_DEG_final.Rda")

#网络图
library(ggnewscale)
List=DEG$log2FoldChange#画图时需要的是离散型的向量，而不是数据框
names(List)=DEG$ENTREZID
head(List)
List=sort(List,decreasing = T)#降序排列
#GO
cnetplot(ego,foldChange=List,circular=TRUE,colorEdge = FALSE,node_label = "all")
#KEGG
cnetplot(kk,foldChange=List,circular=TRUE,colorEdge = FALSE,node_label = "all")


####自己尝试直方图####
#setwd("TCGA-LUAD")
#setwd("PPI")
#PPI<-read.csv(file="network_ppi.csv")
#library(ggplot2)
#ggplot(data=PPI,aes(x=name,y=Degree))+geom_point()+
#  labs(title="ppi",x="name",y="Degree")##直方图或条形图需要的是单变量样本，以统计数量为做Y轴


#COX回归####
setwd("TCGA-LUAD")
setwd("COX")
library(survival)
library(tidyverse)
#install.packages("forestplot")
library(forestplot)
#读取文件
exp_surv_01A<-read.table("exp_surv_01A.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                         ,header=T)
DEG_final<-read.table("DEG_final_venny.txt",sep="\t",check.names=F,stringsAsFactors = F
                      ,header=T)
a<-exp_surv_01A[,DEG_final$SYMBOL]
b<-exp_surv_01A[,1:2]
DEG_exp_surv<-cbind(b,a)
#开始cox单因素分析
Coxoutput<-NULL
for(i in 3:ncol(DEG_exp_surv)){
  g<-colnames(DEG_exp_surv)[i]
  cox<-coxph(Surv(OS.time,OS)~DEG_exp_surv[,i],data=DEG_exp_surv)
  coxSummary=summary(cox)
  
  #修改某列列名，colnames(数据框)[第几列]<-"修改后的列名"
  Coxoutput<-rbind.data.frame(Coxoutput,
                              data.frame(gene=g,
                                         HR=as.numeric(coxSummary$coefficients[,"exp(coef)"])[1],
                                         z=as.numeric(coxSummary$coefficients[,"z"])[1],
                                         pvalue=as.numeric(coxSummary$coefficients[,"Pr(>|z|)"])[1],
                                         lower=as.numeric(coxSummary$conf.int[,3][1]),
                                         upper=as.numeric(coxSummary$conf.int[,4][1]),
                                         stringsAsFactors = F),
                              stringsAsFactors = F)
}

Coxoutput<-arrange(Coxoutput,pvalue)
gene_sig<-Coxoutput[Coxoutput$pvalue<0.005,]#根据文章的要求来进行P值小于0.005进行提取
write.csv(gene_sig,file="gene_sig.csv")
topgene<-gene_sig
##接下来绘制森林图
tabletext<-cbind(c("Gene",topgene$gene),
                 c("HR",format(round(as.numeric(topgene$HR),3),nsmall=3)),
                 c("lower 95%CI",format(round(as.numeric(topgene$lower),3),nsmall=3)),
                 c("upper 95%CI",format(round(as.numeric(topgene$upper),3),nsmall=3)),
                 c("pvalue",format(round(as.numeric(topgene$p),3),nsmall=3)))

#绘制森林图
forest<-forestplot(labeltext=tabletext,
           mean=c(NA,as.numeric(topgene$HR)),
           lower=c(NA,as.numeric(topgene$lower)),
           upper=c(NA,as.numeric(topgene$upper)),
           graph.pos=5,
           graphwidth=unit(.25,"npc"),
           fn.ci_norm="fpDrawDiamondCI",
           col=fpColors(box="#A1D99B",lines="black",zero="#41AB5D"),
           boxsize=0.4,
           zero=1,
           lwd.xaxis=1,
           xlab="a",
           txt_gp=fpTxtGp(label = gpar(cex=1),
                          ticks=gpar(cex=1),
                          xlab=gpar(cex=1),
                          title=gpar(cex=1)),
           hrzl_lines=list("1"=gpar(lwd=1.5,col="black"),
                           "2"=gpar(lwd=1.5,col="black"),
                           "53"=gpar(lwd=1.5,col="black")),
           lineheight=unit(.75,"cm"),
           colgap=unit(0.3,"cm"),
           mar=unit(rep(1.5,times=4),"cm"),
           new_page=F)
forest
#install.packages("forestploter")
#library(forestploter)
#tabletext2<-topgene
#tabletext2$` `<-paste(rep("",7),collapse="")
#p<-forest(data=tabletext2[,c(1,5,6,7,4)],
 #         lower=tabletext2$lower,
  #        upper=tabletext2$upper,
   #       est=tabletext2$HR,
    #      ci_column=4)
#print(p)
#forest1<-edit_plot(forest,
 #                  row=topgene$HR>1,
  #                 col=4,
   #                gp=gpar(col="red"))


setwd("TCGA-LUAD")
setwd("BTK")
library(tidyverse)
tpms01A<-read.table("tpms01A_log2.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
           ,header=T)
tpms11A<-read.table("tpms11A_log2.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
           ,header=T)
gene<-"BTK"
tpms01A<-tpms01A[gene,]
tpms11A<-tpms11A[gene,]
tpms01A<-tpms01A%>%t()%>%as.data.frame()
tpms11A<-tpms11A%>%t()%>%as.data.frame()
write.csv(tpms01A,file="BTK_01A.csv")
write.csv(tpms11A,file="BTK_11A.csv")

#
btk_01_11<-read.csv(file="BTKbar.csv")
library(ggplot2)
library(dplyr)
btk_01_11<-arrange(btk_01_11,desc(X01A_11A))
a<-ggplot(data=btk_01_11,aes(x=X01A_11A,y=BTK))+geom_boxplot()+
  geom_point(position="jitter",
             aes(color=X01A_11A))
a
##
BTK_01A<-read.csv(file="BTK_01A.csv")
BTK_11A<-read.csv(file="BTK_11A.csv")
rownames(BTK_01A)<-NULL
BTK_01A<-BTK_01A%>%as.data.frame%>%column_to_rownames("X")
rownames(BTK_11A)<-NULL
BTK_11A<-BTK_11A%>%as.data.frame%>%column_to_rownames("X")
rownames(BTK_01A)<-substring(rownames(BTK_01A),1,12)
rownames(BTK_11A)<-substring(rownames(BTK_11A),1,12)
BTK_peidui<-cbind(BTK_01A,BTK_11A)

a<-intersect(rownames(BTK_01A),rownames(BTK_11A))
btk_01a<-BTK_01A[a,,drop=F]
btk_11a<-BTK_11A[a,,drop=F]#加上drop=F保证提取出来的数据任然为数据框
BTK_peidui<-cbind(btk_01a,btk_11a)
write.csv(BTK_peidui,file="BTK_peidui.csv")#用EXCEL进行修改
BTK_peidui<-read.csv(file="BTK_peidui.csv")

#画配对图
ggplot(data=BTK_peidui,aes(x=SAMPLE,y=BTK))+
  geom_point(aes(color=SAMPLE))+
  geom_line(aes(group=group))
  



#BTK的生存分析####
setwd("TCGA-LUAD")
setwd("BTK_surv")
surv01A<-read.table("exp_surv_01A.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                    ,header=T)
surv01A$OS.time<-surv01A$OS.time/365#将生存时间从天变为年，利于后续分析





surv01A$group<-ifelse(surv01A$BTK>median(surv01A$BTK),"High","Low")
#该函数添加的high和low都是字符形式的，而非因子形式，故无法做二分类变量判断
#此次的评分分组为关键，分组不同所做的生存分析标准不同


class(surv01A$group)
surv01A$group<-factor(surv01A$group,levels=c("Low","High"))
class(surv01A$group)
table(surv01A$group)

#install.packages("survival")#安装做生存分析的包
library(survival)
fitd<-survdiff(Surv(OS.time,OS)~group,
               data=surv01A,
               na.action = na.exclude)#na.omit：这个操作会从数据中删除任何包含缺失值的行。
#~na.fail：这个操作会在数据中存在缺失值时导致模型失败。
#~na.exclude：这个操作会将包含缺失值的行从分析中排除，但不从数据中删除。
pValue<-1-pchisq(fitd$chisq,length(fitd$n)-1)#该步骤是用来做卡方检验P值
#~与统计学上卡方检验一样，小于0.01才具有统计学意义


#拟合生存曲线
fit<-survfit(Surv(OS.time,OS)~group,data=surv01A)
summary(fit)
p.lab<-paste0("P",ifelse(pValue<0.001,"<0.001",paste0("=",round(pValue,3))))
#该步只是将到的P值规范地写下来
#install.packages("survminer")
library(survminer)
ggsurvplot(
  fit,
  data = surv01A,
  palette = "jco",#配色采用Jco，lancet柳叶刀配色，jama是jama期刊的配色
  pval = p.lab,
  conf.int = TRUE,#显示置信区间
  risk.table = TRUE,#显示风险表
  risk.table.col="strata",
  legend.labs=c("Low","High"),#图例
  size=1,
  xlim=c(0,20),#x轴长度，根据自己数据中OS.time的时间来决定，此数据中最大为19年多，
  #故用20作X轴长度，可用range函数判断自己数据范围
  break.time.by=5,#X轴步长为5
  legend.title="BTK",#TITLE
  surv.median.line = "hv",#限制垂直和水平的中位生存
  ylab="Survival probability(%)",#修改Y轴标签
  xlab="Time(Years)",#修改x轴标签
  ncensor.plot=TRUE,#显示缺失图块
  ncensor.plot.height=0.25,
  risk.table.y.text=FALSE)

##自己代码做分组比较

df<-read.table("clinical.expr01A.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
           ,header=T)
library(tidyverse)
df1<-df[complete.cases(df$ajcc_pathologic_stage),]
table(df1$ajcc_pathologic_stage)
#df4$ajcc_pathologic_m<-gsub("MX",NA,df4$ajcc_pathologic_m)#将TX值换位缺失值
#df1<-df1[complete.cases(df3$ajcc_pathologic_n),]#过滤掉了stage中的缺失值
df1<- df1[order(df1$ajcc_pathologic_stage), ]#使数据按顺序排列，以便于画图
library(ggplot2)
library(RColorBrewer)#这是R语言的调色板，可以选择绘图颜色
display.brewer.all()#查看调色板
mycol<-brewer.pal(n=4,name="Set2")#设置了四个颜色，用Set2为色调
mycol#查看我的调色板的代码名字
mycol<-c("#66C2A5" ,"#FC8D62", "#8DA0CB", "#E78AC3")#设置为集合
library(ggpubr)
stage<-ggboxplot( df1,
                                x="ajcc_pathologic_stage",
                                y="BTK",
                                #color="ajcc_pathologic_stage",#外边框颜色
                                palette=mycol,
                                legend="none",
                                #add="jitter",#添加样本点
                                outlier.shape=NA,
                                fill="ajcc_pathologic_stage",#颜色填充
                                title="stage",
                                xlab="stage",
                                ylab="BTK",
                                size=0.5,)+stat_boxplot(geom="errorbar",width=0.5,size=0.1)+
  geom_boxplot(fill=mycol,size=0.1)+#将中轴线覆盖
  stat_compare_means(method="kruskal.test",label="p.signif")+
  stat_compare_means(comparisons=list(c("Stage I","Stage IV")),method="wilcox",label="p.format")
stage

#t

df1<-df[complete.cases(df$ajcc_pathologic_t),]
table(df1$ajcc_pathologic_t)
df1$ajcc_pathologic_t<-gsub("TX",NA,df1$ajcc_pathologic_t)#将TX值换位缺失值
df1<-df1[complete.cases(df1$ajcc_pathologic_t),]#过滤掉了stage中的缺失值
df1<- df1[order(df1$ajcc_pathologic_t), ]#使数据按顺序排列，以便于画图
table(df1$ajcc_pathologic_t)
library(RColorBrewer)#这是R语言的调色板，可以选择绘图颜色
display.brewer.all()#查看调色板
mycol<-brewer.pal(n=4,name="Set2")#设置了四个颜色，用Set2为色调
mycol#查看我的调色板的代码名字
mycol<-c("#66C2A5" ,"#FC8D62", "#8DA0CB", "#E78AC3")#设置为集合
library(ggpubr)
tclassification <-ggboxplot( df1,
                                x="ajcc_pathologic_t",
                                y="BTK",
                                #color="ajcc_pathologic_stage",#外边框颜色
                                palette=mycol,
                                legend="none",
                                #add="jitter",#添加样本点
                                outlier.shape=NA,
                                fill="ajcc_pathologic_t",#颜色填充
                                title="T Classification",
                                xlab="T",
                                ylab="BTK",
                                size=0.5,)+stat_boxplot(geom="errorbar",width=0.5,size=0.1)+
  geom_boxplot(fill=mycol,size=0.1)+#将中轴线覆盖
  stat_compare_means(method="kruskal.test",label="p.signif")+
  stat_compare_means(comparisons=list(c("T1","T4")),method="wilcox",label="p.format")
tclassification


#n
df1<-df[complete.cases(df$ajcc_pathologic_n),]
table(df1$ajcc_pathologic_n)
df1$ajcc_pathologic_n<-gsub("NX",NA,df1$ajcc_pathologic_n)#将TX值换位缺失值
df1<-df1[complete.cases(df1$ajcc_pathologic_n),]#过滤掉了stage中的缺失值
df1<- df1[order(df1$ajcc_pathologic_n), ]#使数据按顺序排列，以便于画图
table(df1$ajcc_pathologic_n)
library(RColorBrewer)#这是R语言的调色板，可以选择绘图颜色
display.brewer.all()#查看调色板
mycol<-brewer.pal(n=4,name="Set2")#设置了四个颜色，用Set2为色调
mycol#查看我的调色板的代码名字
mycol<-c("#66C2A5" ,"#FC8D62", "#8DA0CB", "#E78AC3")#设置为集合
library(ggpubr)
nclassification <-ggboxplot( df1,
                             x="ajcc_pathologic_n",
                             y="BTK",
                             #color="ajcc_pathologic_stage",#外边框颜色
                             palette=mycol,
                             legend="none",
                             #add="jitter",#添加样本点
                             outlier.shape=NA,
                             fill="ajcc_pathologic_n",#颜色填充
                             title="N Classification",
                             xlab="N",
                             ylab="BTK",
                             size=0.5,)+stat_boxplot(geom="errorbar",width=0.5,size=0.1)+
  geom_boxplot(fill=mycol,size=0.1)+#将中轴线覆盖
  stat_compare_means(method="kruskal.test",label="p.signif")+
  stat_compare_means(comparisons=list(c("N0","N3")),method="wilcox",label="p.format")
nclassification

#m

df1<-df[complete.cases(df$ajcc_pathologic_m),]
table(df1$ajcc_pathologic_m)
df1$ajcc_pathologic_m<-gsub("MX",NA,df1$ajcc_pathologic_m)#将TX值换位缺失值
df1<-df1[complete.cases(df1$ajcc_pathologic_m),]#过滤掉了stage中的缺失值
df1<- df1[order(df1$ajcc_pathologic_m), ]#使数据按顺序排列，以便于画图
table(df1$ajcc_pathologic_m)
library(RColorBrewer)#这是R语言的调色板，可以选择绘图颜色
display.brewer.all()#查看调色板
mycol<-brewer.pal(n=2,name="Set2")#设置了四个颜色，用Set2为色调
mycol#查看我的调色板的代码名字
mycol<-c("#66C2A5" ,"#FC8D62")#设置为集合
library(ggpubr)
mclassification <-ggboxplot( df1,
                             x="ajcc_pathologic_m",
                             y="BTK",
                             #color="ajcc_pathologic_stage",#外边框颜色
                             palette=mycol,
                             legend="none",
                             #add="jitter",#添加样本点
                             outlier.shape=NA,
                             fill="ajcc_pathologic_m",#颜色填充
                             title="M Classification",
                             xlab="M",
                             ylab="BTK",
                             size=0.5,)+stat_boxplot(geom="errorbar",width=0.5,size=0.1)+
  geom_boxplot(fill=mycol,size=0.1)+#将中轴线覆盖
  stat_compare_means(method="kruskal.test",label="p.signif")+
  stat_compare_means(comparisons=list(c("M0","M1")),method="wilcox",label="p.format")
mclassification

#查看四组图
stage
tclassification
nclassification
mclassification



#GESA分析####
setwd("TCGA-LUAD")
setwd("BTK_DEG")
library("BiocManager")
#BiocManager::install("DESeq2")#DESeq2主要包括以下功能：
##数据归一化：DESeq2提供了多种数据归一化方法，如TMM（Trimmed Mean of M-values）归一化。
##差异分析：DESeq2使用负二项分布模型及贝叶斯估计方法，基于基因或转录本的计数数据，计算基因/转录本在不同条件之间的差异表达。差异分析结果包括差异表达基因/转录本的统计显著性和折变倍数等信息。
##数据可视化：DESeq2提供了丰富的数据可视化功能，包括MA图、PCA图、热图等，帮助用户更好地理解数据的特征和差异表达结果
library(DESeq2)
library(tidyverse)
counts_01A<-read.table("counts01A.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                       ,header=T)
###TCGA做差异分析时是使用的counts数据中肿瘤组
exp<-read.table("tpms01A_log2.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                     ,header=T)
x<-"BTK"
med<-median(as.numeric(exp[x,]))
counts_01A2<-as.data.frame(t(counts_01A))
exp2<-as.data.frame(t(exp))
identical(rownames(counts_01A2),rownames(exp2))

conditions=data.frame(sample=rownames(counts_01A2),group=factor(ifelse(exp2[,x]>med,"high","low"),levels=c("low","high")))
conditions<-column_to_rownames(conditions,"sample")   

#接下来是差异分析的准备工作
dds<-DESeqDataSetFromMatrix(
  countData =counts_01A,#计数矩阵，列名为ID必须与colData数据行名一致
  colData = conditions,#要获取的数据的内容，即分组文件(如处理组和对照组)
  desig=~group)#以什么来将所获取的基因分组

#开始差异分析
dds<-DESeq(dds)
#DESeq()内部做了统计分析，我们得到差异基因。所以这也是为什么要用的计数矩阵是原始的计数矩阵，未被标准化处理过
resultsNames(dds)#查看一下结果是high比上low还是反之
res<-results(dds)
save(res,file="BTK_DEG.rda")                    


#数据整理，只要log2fc的分组和基因名即可
library(org.Hs.eg.db)
library(clusterProfiler)


DEG<-as.data.frame(res)
DEG<-arrange(DEG,padj)

DEG<-rownames_to_column(DEG,"SYMBOL")
geneList=DEG[,3]#log2fc进行提取，下面是排序
names(geneList)=as.character(DEG[,"SYMBOL"])
head(geneList)
geneList=sort(geneList,decreasing = TRUE)#降序排列
head(geneList)



#msigdbr包对富集数据提取
#H基因集进行提取分析
#BiocManager::install("msigdbr")
library(msigdbr)
h_gene_sets<-msigdbr(species="Homo sapiens",category="H")%>%
  dplyr::select(gs_name,gene_symbol)
head(h_gene_sets)

gsea<-GSEA(geneList,TERM2GENE =h_gene_sets )
gsea_result_df<-as.data.frame(gsea)
save(gsea,gsea_result_df,file="GSEA_BTK_h.all.rda")

#绘图
library(enrichplot)
gseaplot2(gsea,1,color="red",pvalue_table = T)#数据内的NES是峰值，根据NES分上下弯曲
H<-gseaplot2(gsea,geneSetID = c(1,2,3,4,5,7,8,11),subplots=1:3)
H
H2<-gseaplot2(gsea,geneSetID = c(6,9,10,12,13),subplots=1:3)
H2

#C7（免疫类）基因集进行提取分析

c7_gene_sets<-msigdbr(species="Homo sapiens",category="C7")%>%
  dplyr::select(gs_name,gene_symbol)
head(c7_gene_sets)

gsea<-GSEA(geneList,TERM2GENE =c7_gene_sets )
gsea_result_df<-as.data.frame(gsea)
save(gsea,gsea_result_df,file="GSEA_BTK_c7.all.rda")

#绘图
library(enrichplot)
c7_1<-gseaplot2(gsea,geneSetID =1:7,subplots=1:3)
c7_2<-gseaplot2(gsea,825,color="blue",pvalue_table = T)
c7_1
c7_2



免疫浸润分析#####

#install.packages("devtools")
#library(devtools)
#if(!require(CIBERSORT))devtools::install_github("Moonerss/CIBERSORT")
setwd("TCGA-LUAD")
setwd("CIBERSORT")
#library(CIBERSORT)
#source("CIBERSORT.R")<-CIBERSORT("LM22.txt","tpms01A.txt",perm=1000,QN=T)
install.packages("e1071")
install.packages("parallel")
BiocManager::install("preprocessCore")
library(e1071)
library(parallel)
library(preprocessCore)
library(tidyverse)
source("CIBERSORT.R")
sig_matrix<-"LM22.txt"
mixture_file="tpms01A.txt"
res_cibersort<-CIBERSORT(sig_matrix,mixture_file,perm = 100,QN=TRUE)
res_cibersort<-res_cibersort[,1:22]#将最后三列不需要的统计结果删去
ciber.res<-res_cibersort[,colSums(res_cibersort)>0]
ciber.res<-as.data.frame(ciber.res)
write.table(ciber.res,"ciber.res.txt",sep="\t",row.names=T,col.names=NA,quote=F)


#cibersort彩虹图
library(ggplot2)
mycol<-ggplot2::alpha(rainbow(ncol(ciber.res)),0.7)
par(bty="o",mgp=c(2.5,0.3,0),mar=c(2.1,4.1,2.1,10.1),tcl=-.25,las=1,xpd=F)
barplot(as.matrix(t(ciber.res)),
        border = NA,#柱子边框
        names.arg=rep("",nrow(ciber.res)),
        yaxt="n",
        ylab="Relative percentage",
        col=mycol)+
  axis(side=2,at=c(0,0.2,0.4,0.6,0.8,1),
       labels=c("0%","20%","40%","60%","80%","100%"))+
  legend(par("usr")[2]-20,
         par("usr")[4],
         legend = colnames(ciber.res),
         xpd=T,
         fill=mycol,
         cex=0.6,
         border=NA,
         y.intersp = 1,
         x.intersp = 0.2,
         bty="n")
#分组比较图
a<-ciber.res
exp<-read.table("tpms01A_log2.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                ,header=T)
med=median(as.numeric(exp["BTK",]))
exp<-exp%>%t()%>%as.data.frame()
exp<-exp%>%mutate(group=factor(ifelse(exp$BTK>med,"high","low"),levels=c("low","high")))#mutate是直接给数据框新增一列的函数
class(exp$group)
identical(rownames(a),rownames(exp))
a$group<-exp$group
a<-rownames_to_column(a,"sample")
library(ggsci)
library(tidyr)
library(ggpubr)
b<-gather(a,key=CIBERSORT,value = Fraction,-c(group,sample))
ggboxplot(b,x="CIBERSORT",y="Fraction",
          fill="group",palette="lancet")+
  stat_compare_means(aes(group=group),
                     method="wilcox.test",
                     label="p.signif",
                     symnum.args = list(cutpoints=c(0,0.001,0.01,0.05,1),
                                        symbols=c("***","**","*","ns")))+
  theme(text=element_text(size=10),
        axis.text.x = element_text(angle = 45,hjust=1))


#相关性热土
install.packages("ggstatsplot")
install.packages("ggcorrpolt")
install.packages("corrplot")
library(ggstatsplot)
library(ggcorrplot)
library(corrplot)
cor<-sapply(ciber.res,function(x,y) cor(x,y,method = "spearman"),ciber.res)
rownames(cor)<-colnames(ciber.res)

ggcorrplot(cor,
           hc.order=TRUE,#使用hc.order排序
           type="lower",
           outline.color="white",
           lab=TRUE,
           color=c("#01468b","white","#ee0000"))




library(tidyverse)
exp<-read.table("tpms01A_log2.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                ,header=T)
exp<-exp["BTK",]
ciber<-read.table("ciber.res.txt",sep="\t",row.names=1,check.names=F,stringsAsFactors = F
                ,header=T)
ciber<-ciber%>%t()%>%as.data.frame()
identical(colnames(ciber),colnames(exp))
exp_ciber<-rbind(exp,ciber)
exp_ciber<-exp_ciber%>%t()%>%as.data.frame()
colnames(exp_ciber)<-gsub(" ",".",colnames(exp_ciber))#将空格变为点才能画图
install.packages("ggside")
library(ggstatsplot)
library(ggside)
ggscatterstats(data=exp_ciber,
               y=BTK,
               x=B.cells.naive,
               type="nonparametric",
               margins="both",
               xfill="blue",
               yfill="red",
               marginal.type="densigram")