CPTAC数据库(https://cptac-data-portal.georgetown.edu/datasets)
先点进去对应的study下载临床数据看是否有自己想研究的信息,然后再下载对应的蛋白数据。
然后对表达矩阵进行整理,并比对需要的临床信息。
表达矩阵整理
library(impute)#用来补缺
library(limma)#用来矫正
rt=read.table('UCEC-PROTEIN.txt',sep="\t",header=T,check.names=F,row.names=NULL)
rt <- rt[!duplicated(rt$Gene),]
rownames(rt)<-rt[,1]
rt<-rt[,-1]
rt=as.matrix(rt)
#保留unshared列,从第2列到最后一列取偶数列。
selectCol=seq(2,ncol(rt),2)
rt=rt[,selectCol]
rt<-as.data.frame(rt)
#计数每一行的NA值
f<-function(x) sum(is.na(x))
a<-as.data.frame(apply(rt,1,f)) #计数每行
colnames(a)<-'NA'
rt<-cbind(a,rt)
#按照NA列降序排序,然后输出na>50%样本数的基因。
rt<-rt[order(rt$"NA",decreasing = T),]
rt<-rt[rt$'NA'<76,] #需要修改
rt<-rt[,-1]
#数据补缺,需要是矩阵
rt<-as.matrix(rt)
mat=impute.knn(rt)
rt=mat$data
#对数据进行矫正
normalizeData=normalizeBetweenArrays(as.matrix(rt))
#输出结果
normalizeData=rbind(geneNames=colnames(normalizeData),normalizeData)
data<-as.data.frame(normalizeData)
write.table(data,file="normalize.txt",sep="\t",quote=F,col.names=F)
差异分析
#先手动整理得到E3-EXP-GROUP文件
pFilter=0.05
logFCfilter=0 #不对logfc进行过滤,只按照P值进行过滤
conNum=31 #正常样品数目
treatNum=100 #肿瘤样品数目
#读取输入文件
outTab=data.frame()
group=c(rep(1,conNum),rep(2,treatNum))
data=read.table('E3-EXP-GROUP.txt',sep="\t",header=T,check.names=F,row.names=1)
data=as.matrix(data)
#差异分析
for(i in row.names(data)){
geneName=i
rt=rbind(expression=data[i,],group=group)
rt=as.matrix(t(rt))
tTest<-t.test(expression ~ group, data=rt) #蛋白质组大多进行T检验
pvalue=tTest$p.value
conGeneMeans=mean(data[i,1:conNum])
treatGeneMeans=mean(data[i,(conNum+1):ncol(data)])
logFC=treatGeneMeans-conGeneMeans
conMed=median(data[i,1:conNum])
treatMed=median(data[i,(conNum+1):ncol(data)])
diffMed=treatMed-conMed
if( ((logFC>0) & (diffMed>0)) | ((logFC<0) & (diffMed<0)) ){
outTab=rbind(outTab,cbind(gene=i,conMean=conGeneMeans,treatMean=treatGeneMeans,logFC=logFC,pValue=pvalue))
}
}
#输出所有蛋白的差异情况
write.table(outTab,file="all-UCEC.xls",sep="\t",row.names=F,quote=F)