install.packages("remotes")
remotes::install_github("robinmeyers/niboR")
library(niboR)
下午有个客户用GEO的数据放到TCGA的脚本运行,查了半天才得知这个问题
最后帮他改通脚本,但跑出来的结果行名不对
再次查找问题,发现是.gct文件的行名没有改成gsm的
以下是更改后的代码:
install.packages("remotes")
remotes::install_github("robinmeyers/niboR")
library(niboR)
library(limma)
library(estimate)
#setwd("C:\\Users\\biowolf\\Desktop\\ssGSEA\\10.estimate") #设置工作目录
inputFile="symbol.txt" #输入文件名字
#读取文件
rt=read.table(inputFile,sep="\t",header=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)
#删除正常,只保留肿瘤样品
#group=sapply(strsplit(colnames(data),"\\-"),"[",4)
#group=sapply(strsplit(group,""),"[",1)
#group=gsub("2","1",group)
#data=data[,group==0]
#out=data[rowMeans(data)>0,]
#out=rbind(ID=colnames(out),out)
out <- data[,149:246]
#输出整理后的矩阵文件
write.table(out,file="uniq.symbol.txt",sep="\t",quote=F,col.names=F)
#运行estimate包
filterCommonGenes(input.f="uniq.symbol.txt",
output.f="commonGenes.gct",
id="GeneSymbol")
ds <- read.gct("commonGenes.gct")
colnames(ds)<- (colnames(exp)[149:246])
write.gct(ds,"commonGenes.gct")
estimateScore(input.ds = "commonGenes.gct",
output.ds="estimateScore.gct")
#输出每个样品的打分
scores=read.table("estimateScore.gct",skip = 2,header = T)
rownames(scores)=scores[,1]
scores=t(scores[,3:ncol(scores)])
rownames(scores)=gsub("\\.","\\-",rownames(scores))
out=rbind(ID=colnames(scores),scores)
write.table(out,file="scores.txt",sep="\t",quote=F,col.names=F)
网友评论