蛋白分析相关

2019-07-17  本文已影响0人  存存baby

找差异蛋白

rm(list=ls()) #清空环境

pcut=0.05

fccut=1.5

MaxNA=1  #生物学重复最大允许几个缺失

#第一列是蛋白名 2、3、4列是生物学重复

#把EXCEL保存为txt文件(另存为-制表符分割)

oridata=read.table("data/demo.pro.foldChange.txt",header=T,row.name=1)

#将0和10的蛋白表达量均替换为NA

oridata[oridata==0]=NA

oridata[oridata==10]=NA

#按行(蛋白)计算NA的数目,取0或者1个NA的蛋白

na_count=-apply(oridata,1,function(x)sum(is.na(x)))

matrix=oridata[na_count<=MaxNA,]

matrix=as.data.frame(matrix)

trData=log2(matrix)  #对FC取对数

pvalue=NULL

#循环 将每个蛋白的log2(FC)与均值为0(即FC=1)做t检验,取p.value

for(i in (1:nrow(trData))){pvalue[i]=t.test(trData[i,],mu=0)$p.value}

FCmean=rowMeans(matrix,na.rm = T) #计算每个蛋白的FC均值

out=cbind(matrix,FCmean,pvalue)  #输出结果包含原始矩阵、FC值、t检验的p值

write.csv(out,file="demo/out.all.csv",quote = F)

#abs是绝对值 log2X x大于1则是正值 X小于1则是负值

FCindex=abs(log2(FCmean))>abs(log2(fccut))  #取FC均值大于阈值的蛋白

pindex=pvalue<pcut

difindex=FCindex & pindex #dif表示对于FC和P值同时满足的

diff=out(difindex,)  #选出满足的,付给新变量diff

write.csv(diff,file="demo.out.diff.csv",quote=F)

上一篇下一篇

猜你喜欢

热点阅读