CPTAC蛋白质差异分析详细步骤

CPTAC蛋白质差异分析

1. CPTAC数据下载

根据自己的需求选择对应的癌症，然后下载，这里就不展开讲了。以下为下载的数据目录

2. 数据补缺校正

从CPTAC下载的蛋白质数据很多都是有缺失值的，这样不利于我们后期分析，所以需要对这些缺失值进行补缺，然后校正（归一化）。

• 首先新建空文件夹名为1_数据补缺校正，将下载的第二个文件（肿瘤文件）复制到空文件夹中

• 然后新建的txt文件，名为input.txt，将复制到空文件夹的表格内容粘贴到input.txt中，然后保存

• 之后在新建的空文件夹新建脚本，脚本代码如下所示：

  if (!requireNamespace("BiocManager", quietly = TRUE))
     install.packages("BiocManager")
  BiocManager::install("impute")
  
  if (!requireNamespace("BiocManager", quietly = TRUE))
     install.packages("BiocManager")
  BiocManager::install("limma")
  
  library(impute)
  library(limma)
  setwd("C:\\Users\\Administrator\\Desktop\\cptac\\1_数据补缺校正")           #设置工作目录
  file="input.txt"                                                  #输入文件名字       
  
  #读取file并对基因整理
  rt=read.table(file,sep="\t",header=T,check.names=F,row.names=1)
  rt=as.matrix(rt)
  
  #选取数值,需要修改
  #selectCol=seq(2,ncol(rt),2)        #保留列号为偶数的列，如果需要保留列号为奇数的列，把第一个2改成1
  #rt=rt[,selectCol]
  
  #数据补缺
  mat=impute.knn(rt)
  rt=mat$data
  
  #对数据进行矫正
  normalizeData=normalizeBetweenArrays(as.matrix(rt))
  
  #输出结果
  normalizeData=rbind(geneNames=colnames(normalizeData),normalizeData)
  write.table(normalizeData,file="normalize.txt",sep="\t",quote=F,col.names=F)
  

• 最后将上述代码复制到R中运行，运行完在名为1_数据补缺校正的文件夹中就会有normalize.txt了。至此就实现了数据补缺和校正

3. 数据分组

差异分析肯定是对某个类别进行的分析，比如对性别、癌症的诊断方式等，这里我们以癌症患者的性别进行分组演示。

• 首先新建文件夹为2_分组分析，将数据补缺校正得到的normalize.txt文件和数据库下载的临床文件复制到该文件夹中

• 打开临床文件，对性别这一列进行排序，得到的结果如下图所示：

• 接着新建文本文件sample1.txt和sample2.txt，将排序之后的临床文件的female对应的id复制粘贴到sample1.txt中，将male对应的id复制粘贴到sample2.txt中

• 然后建立perl脚本文件cptac_group.pl，在此文件夹目录下打开powershell，输入perl cptac_group.pl，（此脚本太大了，代码就先不放了）程序运行完即可显示sample1和sample2对应的样本数，并且在文件夹下按照sample1和sample2的顺序生成了新的样本表达文件sampleExp.txt
  

• 至此，样本分组完毕

4. 蛋白质组差异分析

在进行完上面三步之后，接下来就可以进行差异分析了，步骤如下：

• 新建空白文件夹3_差异分析，并将数据分组之后的结果sampleExp.txt复制到该文件夹下

• 新建R语言脚本cpta_diff.R，其代码如下：

  setwd("C:\\Users\\Administrator\\Desktop\\cptac\\3_差异分析")            #设置工作目录
  inputFile="sampleExp.txt"                                     #输入文件
  pFilter=0.05                                                  #p值临界值
  logFCfilter=0                                                 #logFC临界值
  conNum=56                                                    #sample1组样品数目
  treatNum=41                                                  #sample2组样品数目
  
  #读取输入文件
  outTab=data.frame()
  group=c(rep(1,conNum),rep(2,treatNum))
  data=read.table(inputFile,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)
    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.xls",sep="\t",row.names=F,quote=F)
  
  #输出差异表格
  outDiff=outTab[( abs(as.numeric(as.vector(outTab$logFC)))>logFCfilter & as.numeric(as.vector(outTab$pValue))<pFilter),]
  write.table(outDiff,file="diff.xls",sep="\t",row.names=F,quote=F)
  
  #输出差异蛋白表达
  heatmap=rbind(ID=colnames(data[as.vector(outDiff[,1]),]),data[as.vector(outDiff[,1]),])
  write.table(heatmap,file="diffProteinExp.txt",sep="\t",col.names=F,quote=F)
  

• 然后打开R软件，粘贴此代码运行，最终在文件夹下即可得到结果，如下图所示：
  

  其中：

  • all.xls为所有基因的差异数据，如下图所示：

  

  • diff.xls为满足p值小于0.05的蛋白质的差异数据，如下图所示：

    

    • diffProteinExp.txt为差异蛋白的表达量数据，如下图所示：

      

至此，蛋白质组的差异分析就结束了，如果需要作图，用表格的差异数据在graphpad可自行作图。