一、生信分析的目的

TCGA的表达矩阵基因多、样本多，筛选的时候难免会受到极端基因的影响（标准化数据后可能也会有影响）。并且筛选之后的基因又非常的多。
今天就给同学们提供一个解决的思路——利用热点：铜死亡
我给同学们搜集了4篇论文的铜死亡基因（同学们也可以选择铁死亡基因或线粒体基因等进行操作）

二、具体的代码

remove(list = ls()) ##清空当前环境
setwd("D:/生信代码复现/27.铜死亡相关基因筛选/27.铜死亡相关基因筛选") ##💚💚💚设置路径

  
library(limma)            
expFile="combined_RNAseq_FPKM.txt"      #💛💛💛表达数据文件
geneFile="gene.txt"       #💛💛💛基因列表文件，这些是和铜死亡相关的文献


#读取输入文件，并对数据进行处理
rt=read.table(expFile, header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[rowMeans(data)>0,]

#读取基因列表文件,提取铜死亡相关基因的表达量
gene=read.table(geneFile, header=F, sep="\t", check.names=F)
sameGene=intersect(as.vector(gene[,1]), rownames(data))
geneExp=data[sameGene,]

#输出结果
outTab=rbind(ID=colnames(geneExp),geneExp)
write.table(outTab, file="铜死亡基因表达量.txt", sep="\t", quote=F, col.names=F)

#####共表达分析
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""), "[", 1)
group=gsub("2","1",group)
RNA=data[,group==0]
conNum=length(group[group==1])       #正常组样品数目
treatNum=length(group[group==0])     #肿瘤组样品数目
sampleType=c(rep(1,conNum), rep(2,treatNum))


rt1=rt1=read.table("铜死亡基因表达量.txt", header=T, sep="\t", check.names=F)
rt1=as.matrix(rt1)
rownames(rt1)=rt1[,1]
exp1=rt1[,2:ncol(rt1)]
dimnames1=list(rownames(exp1),colnames(exp1))
cuproptosis=matrix(as.numeric(as.matrix(exp1)), nrow=nrow(exp1), dimnames=dimnames1)
cuproptosis=avereps(cuproptosis)
cuproptosis=cuproptosis[rowMeans(cuproptosis)>0,]

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

#相关性检验
corFilter=0.3           
pvalueFilter=0.05  
outTab=data.frame()
for(i in row.names(RNA)){
  if(sd(RNA[i,])>0.1){
    test=wilcox.test(data[i,] ~ sampleType)
    if(test$p.value<0.05){
      for(j in row.names(cuproptosis)){
        x=as.numeric(RNA[i,])
        y=as.numeric(cuproptosis[j,])
        corT=cor.test(x,y)
        cor=corT$estimate
        pvalue=corT$p.value
        if((cor>corFilter) & (pvalue<pvalueFilter)){
          outTab=rbind(outTab,cbind(Cuproptosis=j,RNA=i,cor,pvalue,Regulation="postive"))
        }
        if((cor< -corFilter) & (pvalue<pvalueFilter)){
          outTab=rbind(outTab,cbind(Cuproptosis=j,RNA=i,cor,pvalue,Regulation="negative"))
        }
      }
    }
  }
}

#输出相关性的结果
write.table(file="Result.txt",outTab,sep="\t",quote=F,row.names=F)

#提取铜死亡相关RNA的表达量
cuproptosisRNA=unique(as.vector(outTab[,"RNA"]))
cuproptosisRNAexp=data[cuproptosisRNA,]
cuproptosisRNAexp=rbind(ID=colnames(cuproptosisRNAexp), cuproptosisRNAexp)
write.table(cuproptosisRNAexp,file="cuproptosisExp.txt",sep="\t",quote=F,col.names=F)

三、最终的效果

①R语言中的效果：

②文件效果：

四、遇到的问题

无

五、参考

①具体代码参考——公众号《叉叉滴同学的生信笔记》