返回生信工具
序列分析

利用生存分析进一步筛选基因

1.作用

差异分析筛选出基因之后再进行生存分析以进一步筛选出目标基因。 准备工作： 1.我们得到差异基因的表达量后还需准备TCGA的临床数据文件。从TCGA数据库下载之后进行整理。 2.将临床数据和差异基因表达数据合并

2.代码

[!NOTE]- 具体的实现代码，在源代码的基础上修改较多

#install.packages("survival")
library(tidyverse)
library(survival)
library(limma)
setwd("D:/生信代码复现/3.生存分析/3.生存分析")  #💚💚💚💚💚💚💚💚💚💚工作目录（需修改）
expFile="diffGeneExp.txt"  #❗❗❗❗❗❗❗❗❗❗❗❗这个文件就是差异基因表达的文件，但是很奇怪的是按照教程来弄，我们这个得到的文件的格式是excel的格式，而这里需要的是txt格式
cliFile="time.txt"  #❗❗❗❗❗❗❗❗❗❗❗❗这个临床的生存时间分析暂时也是不知道从哪里弄出来的
rt2=read.table(expFile, header=T, check.names=F, row.names=1)
exp_data_T = rt2%>% dplyr::select(str_which(colnames(.), "-01A")) # 匹配列名或用下示写法
nT = ncol(exp_data_T) 
#读取生存数据文件
cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)
cli$futime=cli$futime/365
group1=sapply(strsplit(colnames(rt2),"\\-"), "[", 4)
group1=sapply(strsplit(group1,""), "[", 1)
group1=gsub("2", "1", group1)
conNum=length(group1[group1==1])       #正常组样品数目
treatNum=length(group1[group1==0])     #肿瘤组样品数目
Type=c(rep(1,conNum), rep(2,treatNum))
#删掉正常样品
tumorData=exp_data_T
tumorData=as.matrix(tumorData)
tumorData=t(tumorData)
rownames(tumorData)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(tumorData))
data=avereps(tumorData)
#数据合并并输出结果
sameSample=intersect(row.names(data), row.names(cli))
data=data[sameSample,,drop=F]
cli=cli[sameSample,,drop=F]
rt2=cbind(cli, data)
#输出合并后的数据
outTab=cbind(ID=row.names(rt2), rt2)
outTab=outTab[,-ncol(outTab)]


#差异分析
pFilter=0.05                                                                 #显著性过滤条件
rt=outTab    #读取输入文件                                                   #如果以月为单位，除以30；以年为单位，除以365
outTab=data.frame()
sigGenes=c("futime","fustat")
for(gene in colnames(rt[,4:ncol(rt)]))
{
  if(sd(rt[,gene])<0.1){
    next}
  a=rt[,gene]<=median(rt[,gene])
  
  rt1=rt[a,]
  b=setdiff(rownames(rt),rownames(rt1))
  rt2=rt[b,]
  surTab1=summary(survfit(Surv(futime, fustat) ~ 1, data = rt1))
  surTab2=summary(survfit(Surv(futime, fustat) ~ 1, data = rt2))
  survivalTab1=cbind(time=surTab1$time, surv=surTab1$surv,lower=surTab1$lower,upper=surTab1$upper)
  survivalTab1=survivalTab1[survivalTab1[,"time"]<5,]

  survivalTab2=cbind(time=surTab2$time, surv=surTab2$surv,lower=surTab2$lower,upper=surTab2$upper)
  survivalTab2=survivalTab2[survivalTab2[,"time"]<5,]

  fiveYearsDiff=abs(survivalTab1["surv"]-survivalTab2["surv"])
  
  #km方法
  diff=survdiff(Surv(futime, fustat) ~a,data = rt)
  pValue=1-pchisq(diff$chisq,df=1)
  fit=survfit(Surv(futime, fustat) ~ a, data = rt)
  #cox方法
  cox=coxph(Surv(futime, fustat) ~ rt[,gene], data = rt)
  coxSummary = summary(cox)
  coxP=coxSummary$coefficients[,"Pr(>|z|)"]
  
  if(pValue<pFilter) {
    sigGenes=c(sigGenes,gene)
    outTab=rbind(outTab, 
                 cbind(gene=gene,
                       KM=pValue,
                       HR=coxSummary$conf.int[,"exp(coef)"],
                       HR.95L=coxSummary$conf.int[,"lower .95"],
                       HR.95H=coxSummary$conf.int[,"upper .95"],
                       coxPvalue=coxP) )
  }
}
write.table(outTab,file="survival.xls",sep="\t",row.names=F,quote=F)    #输出基因和p值表格文件
surSigExp=rt[,sigGenes]
surSigExp=cbind(id=row.names(surSigExp),surSigExp)
write.table(surSigExp,file="surSigExp.txt",sep="\t",row.names=F,quote=F)

3.效果

文件夹的效果 image.png

R语言中的效果 image.png

4.遇到的问题

1.首先遇到的问题就是这个time.txt文件我就找不到，查阅了网站以后了解到了这个

[!NOTE]- 如何获得time.txt文件 image.png

[!NOTE]- 矩阵长度超过1，还有就是最后的生存分析出来的表格啥也没有 报错了，无法读取我的文件 image.png 问了人工智能

# 如果没有安装'survival'包，取消下一行代码的注释来安装它
#install.packages("survival")

library(tidyverse)  # 加载'tidyverse'包，这是一组用于数据科学的R包的集合，包含数据操作、绘图等功能的包
library(survival)   # 加载'survival'包，提供了生存分析的函数和模型
library(limma)      # 加载'limma'包，它主要用于生物信息学领域的表达数据分析

# 设定工作目录，这一行代码被注释掉了，需要运行时根据实际路径取消注释并修改路径
#setwd("D:/生信代码复现/3.生存分析/3.生存分析")  # 设置工作目录（需根据实际情况修改）

expFile="diffGeneExp.txt"  # 指定差异基因表达数据文件的名称，注意需确保文件是文本格式，而不是Excel格式
cliFile="time.txt"         # 指定临床生存时间数据文件的名称，这里先做一个占位，具体文件从哪里获取还不明确

# 使用read.table函数读取名为diffGeneExp.txt的表达数据文件
# header=T表示文件的第一行是列名，check.names=F意味着不检查变量名是否为合法的R变量名，row.names=1表示文件的第一列是行名
rt2=read.table(expFile, header=T, check.names=F, row.names=1)

结果我发现：

因为还没有弄清楚，所以打算对代码的含义进行深层次理解。

# install.packages("survival") # 如果未安装survival包，取消注释该行以安装
library(tidyverse) # 载入tidyverse包，用于数据整理和可视化
library(survival) # 载入survival包，用于生存分析
library(limma) # 载入limma包，用于生物信息学和统计学分析
setwd("D:/生信代码复现/3.生存分析/3.生存分析")  # 设置工作目录，这里需要根据实际情况修改路径
expFile="diffGeneExp.txt"  # 定义差异基因表达文件的路径，需要的格式是文本文件（txt）
cliFile="time.txt"  # 定义临床生存时间分析数据文件的路径，来源暂时未知
rt2=read.table(expFile, header=T, check.names=F, row.names=1) # 读取差异基因表达数据
exp_data_T = rt2 %>% dplyr::select(str_which(colnames(.), "-01A")) # 选择包含"-01A"的列
nT = ncol(exp_data_T) # 计算所选列的数量

# 读取临床数据文件
cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1) # 读取临床数据
cli$futime=cli$futime/365 # 将临床数据中的随访时间从天转换为年
group1=sapply(strsplit(colnames(rt2), "\\-"), "[", 4) # 提取列名中的第四部分用于分组
group1=sapply(strsplit(group1, ""), "[", 1) # 进一步提取分组信息
group1=gsub("2", "1", group1) # 替换分组标识，将"2"替换为"1"
conNum=length(group1[group1==1]) # 计算正常组样本数量
treatNum=length(group1[group1==0]) # 计算肿瘤组样本数量
Type=c(rep(1,conNum), rep(2,treatNum)) # 创建一个向量记录正常组和肿瘤组的样本编号

# 去除正常样品数据
tumorData=exp_data_T # 将差异基因表达数据赋值给肿瘤数据变量
tumorData=as.matrix(tumorData) # 将数据框转换为矩阵
tumorData=t(tumorData) # 转置矩阵，使得行为样品，列为基因
rownames(tumorData)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(tumorData)) # 将矩阵的行名进行格式化处理
data=avereps(tumorData) # 对矩阵中的肿瘤数据执行求均值操作（limma包中的函数）

# 合并数据并输出结果
sameSample=intersect(row.names(data), row.names(cli)) # 找出既在差异基因表达数据中也在临床数据中的样本
data=data[sameSample,,drop=F] # 选取这些共有样本的差异基因表达数据
cli=cli[sameSample,,drop=F] # 选取这些共有样本的临床数据
rt2=cbind(cli, data) # 将临床数据和差异基因表达数据合并

# 输出合并后的数据
outTab=cbind(ID=row.names(rt2), rt2) # 合并数据后在最前面增加一个“ID”列用于标记样本ID
outTab=outTab[,-ncol(outTab)] # 删除最后一个冗余的列


# 差异分析部分
pFilter=0.05 # 设定显著性水平的阈值为0.05
rt=outTab # 将之前合并好的表格数据赋值给rt变量
outTab=data.frame() # 创建一个新的空数据框，用于存放差异分析结果
sigGenes=c("futime","fustat") # 初始化一个向量，用于存放显著的基因名，包括生存时间和状态

# 遍历表格中所有的基因
for(gene in colnames(rt[,4:ncol(rt)])) {
   if(sd(rt[,gene])<0.1) { # 如果基因的表达标准差小于0.1，则跳过这个基因
     next
   }
   
   a=rt[,gene]<=median(rt[,gene]) # 创建一个逻辑向量，基因表达低于中位数的为TRUE
  
   rt1=rt[a,] # 根据基因表达低于中位数的条件筛选数据
   b=setdiff(rownames(rt), rownames(rt1)) # 找出不在rt1中的样本
   rt2=rt[b,] # 根据上面找出的样本筛选数据
   surTab1=summary(survfit(Surv(futime, fustat) ~ 1, data = rt1)) # 对表达水平低的组进行生存分析并获取摘要
   surTab2=summary(survfit(Surv(futime, fustat) ~ 1, data = rt2)) # 对表达水平高的组进行生存分析并获取摘要
   survivalTab1=cbind(time=surTab1$time, surv=surTab1$surv,lower=surTab1$lower,upper=surTab1$upper) # 将分析结果整合到一个表中
   survivalTab1=survivalTab1[survivalTab1[, "time"]<5,] # 仅保留生存时间小于5年的数据 
   if(class(survivalTab1) == "matrix") {
     survivalTab1=survivalTab1[nrow(survivalTab1),] # 确保结果是一个矩阵
   }
   survivalTab2=cbind(time=surTab2$time, surv=surTab2$surv,lower=surTab2$lower,upper=surTab2$upper)
   survivalTab2=survivalTab2[survivalTab2[, "time"]<5,]
   if(class(survivalTab2) == "matrix") {
     survivalTab2=survivalTab2[nrow(survivalTab2),]
   }
   fiveYearsDiff=abs(survivalTab1["surv"]-survivalTab2["surv"]) # 计算两组之间在五年生存率上的差异 

   # 使用Kaplan-Meier方法进行差异测试
   diff=survdiff(Surv(futime, fustat) ~ a, data = rt) # 进行生存差异性检验
   pValue=1-pchisq(diff$chisq, df=1) # 计算卡方检验的p值
   fit=survfit(Surv(futime, fustat) ~ a, data = rt) # 使用Kaplan-Meier估计生存曲线
   # 使用Cox比例风险模型进行差异分析
   cox=coxph(Surv(futime, fustat) ~ rt[,gene], data = rt) # 拟合Cox比例风险模型
   coxSummary = summary(cox) # 获取Cox模型的摘要统计
   coxP=coxSummary$coefficients[, "Pr(>|z|)"] # 提取Cox模型的p值
 
   # 如果Kaplan-Meier方法和Cox模型的p值均小于阈值，且五年生存率差异大于0.15，则认为该基因的表达差异与生存时间有显著相关性
   if((pValue < pFilter) & (coxP < pFilter) & (fiveYearsDiff > 0.15)) {
       sigGenes=c(sigGenes, gene) # 将该基因添加到显著基因列表中
       outTab=rbind(outTab, cbind(gene=gene, KM=pValue, HR=coxSummary$conf.int[, "exp(coef)"], HR.95L=coxSummary$conf.int[, "lower .95"], HR.95H=coxSummary$conf.int[, "upper .95"], coxPvalue=coxP)) # 将结果添加到输出表中
   }
}
# 把最终的差异分析的结果写入一个Excel文件
write.table(outTab, file="survival.xls", sep="\t", row.names=F, quote=F)
# 把显著基因的表达数据和id一起写入文本文件
surSigExp=rt[,sigGenes]
surSigExp=cbind(id=row.names(surSigExp), surSigExp)
write.table(surSigExp, file="surSigExp.txt", sep="\t", row.names=F, quote=F)

最后问了原作者，发现需要将类似这种代码的删掉 if(class(survivalTab1) == "matrix") { survivalTab1=survivalTab1[nrow(survivalTab1),] # 确保结果是一个矩阵 } 因为最新的R已经不需要这样的确认方式了。

最后的生存分析啥也没有 需要将if((pValue < pFilter) & (coxP < pFilter) & (fiveYearsDiff > 0.15))仅保留if((pValue < pFilter)