2020-02-27 R读取.gct文件

2020-02-27  本文已影响0人  my_derek
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)
上一篇 下一篇

猜你喜欢

热点阅读