CPTAC數(shù)據(jù)庫(https://cptac-data-portal.georgetown.edu/datasets)
image.png
可以查看各種研究,有Asmb的就是有處理好的表達(dá)矩陣
先點(diǎn)進(jìn)去對(duì)應(yīng)的study下載臨床數(shù)據(jù)看是否有自己想研究的信息相味,然后再下載對(duì)應(yīng)的蛋白數(shù)據(jù)拾积。
點(diǎn)進(jìn)來對(duì)臨床信息進(jìn)行下載
然后對(duì)需要的蛋白數(shù)據(jù)進(jìn)行下載
下載出來的只有tmt的是蛋白表達(dá)數(shù)據(jù)
然后對(duì)表達(dá)矩陣進(jìn)行整理,并比對(duì)需要的臨床信息。
表達(dá)矩陣整理
library(impute)#用來補(bǔ)缺
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列到最后一列取偶數(shù)列拓巧。
selectCol=seq(2,ncol(rt),2)
rt=rt[,selectCol]
rt<-as.data.frame(rt)
#計(jì)數(shù)每一行的NA值
f<-function(x) sum(is.na(x))
a<-as.data.frame(apply(rt,1,f)) #計(jì)數(shù)每行
colnames(a)<-'NA'
rt<-cbind(a,rt)
#按照NA列降序排序斯碌,然后輸出na>50%樣本數(shù)的基因。
rt<-rt[order(rt$"NA",decreasing = T),]
rt<-rt[rt$'NA'<76,] #需要修改
rt<-rt[,-1]
#數(shù)據(jù)補(bǔ)缺肛度,需要是矩陣
rt<-as.matrix(rt)
mat=impute.knn(rt)
rt=mat$data
#對(duì)數(shù)據(jù)進(jìn)行矯正
normalizeData=normalizeBetweenArrays(as.matrix(rt))
#輸出結(jié)果
normalizeData=rbind(geneNames=colnames(normalizeData),normalizeData)
data<-as.data.frame(normalizeData)
write.table(data,file="normalize.txt",sep="\t",quote=F,col.names=F)
差異分析
#先手動(dòng)整理得到E3-EXP-GROUP文件
pFilter=0.05
logFCfilter=0 #不對(duì)logfc進(jìn)行過濾傻唾,只按照P值進(jìn)行過濾
conNum=31 #正常樣品數(shù)目
treatNum=100 #腫瘤樣品數(shù)目
#讀取輸入文件
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) #蛋白質(zhì)組大多進(jìn)行T檢驗(yàn)
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)