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", "#CD"))(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", "#CD"))(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")