1.导入数据

setwd('D:/陈思杰/大创/稳定对')
library('GEOquery')
rm(list = ls())
##数据导入
data=getGEO('GSE113486',destdir='.')[[1]]
info=pData(data) #标签信息：癌症分类信息在第36列
info=as.character(info[,36])
exp=exprs(data) #表达谱矩阵
miRNA=data@featureData@data$miRNA_ID_LIST #ID注释
exp=exp[-grep(',',miRNA),]
miRNA=miRNA[-grep(',',miRNA)]
# 这里没有做获得miRNA名字并转化，如MIMAT000062 变为 62，是因为下面的逆转对代码我加了一行进行修改

2.去离群值

#离群值处理函数定义
func=function(exp,name,info){
  cancer=exp[,info==name] #选取乳腺癌样本
  set.seed(2021) #随机种子，可更改任意值
  if(name=='Bladder Cancer'){
    cancer=cancer[,sample(1:ncol(cancer), 40, replace=FALSE)]#随机取40个样本
  }#这个if语句使392例膀胱癌随机抽取40个样本
  R=cor(cancer) #相关系数矩阵
  means=apply(R,1,mean) #各行均值
  means2=mean(means) #各行均值的均值
  sds=sd(means)  #各行均值的标准差
  a=c(means<means2-2*sds) #均值小于均值的均值减去2倍标准差的样本算离群
  cat(a+0,end='\n') #查看异常值数量结果：有x个样本异常值
  cat('共',sum(a),'个样本离群',sep='',end='\n')
  #因此对于Breast Cancer而言，cancer只选择40-sum(a)个样本来做分析
  cancer=cancer[,!a]
}

#每种数据去离群（392例膀胱癌随机挑40例）
exp1=func(exp,'Bladder Cancer',info)
exp2=func(exp,'Non-cancer control',info)
exp3=func(exp,'Sarcoma',info)
exp4=func(exp,'Prostate Cancer',info)
exp5=func(exp,'Pancreatic Cancer',info)
exp6=func(exp,'Ovarian Cancer',info)
exp7=func(exp,'Lung Cancer',info)
exp8=func(exp,'Hepatocellular Carcinoma',info)
exp9=func(exp,'Glioma',info)
exp10=func(exp,'Gastric Cancer',info)
exp11=func(exp,'Esophageal Cancer',info)
exp12=func(exp,'Colorectal Cancer',info)
exp13=func(exp,'Breast Cancer',info)
exp14=func(exp,'Biliary Tract Cancer',info)

3.找稳定对

#100个正常样本中稳定存在的秩序关系(a>b)的miRNA对称为稳定对
#稳定对自定义函数
wdd=function(exp,id,Thr){ 
  star_time <- Sys.time() 
  exp=exp-min(exp)+1
  num.simple=ncol(exp)
  num.mirna=nrow(exp)
  k=matrix(0,num.mirna,num.mirna)
  pb <- txtProgressBar(style=3)
  for(i in 1:num.simple){
    setTxtProgressBar(pb, i/num.simple)
    a=exp[,i]
    a=a%*%t(1/a)
    a[a<=1]=0
    a[a>1]=1
    k=k+a
  }
  k=k/num.simple
  index=which(k>=Thr)
  count=k[k>=Thr]
  miRNA_a=id[index%%num.mirna]   #取余
  miRNA_b=id[index%/%num.mirna+1]  #取整
  Sta_pair=data.frame(miRNA_a,miRNA_b,count)
  end_time <- Sys.time()
  cat('\n 程序用时:')
  print(end_time-star_time)
  Sta_pair
}
#函数调用
dataN=exp2
spair=wdd(dataN,1:2550,0.9)

4.找逆转对

#你研究的这种癌症中这种秩序关系变成(a<b)了，则这个稳定对叫做逆转对
#逆转对代码
Reverse_pair<-function(DataN,DataT,GID,Pairs,Step,savename){
  #首先挑选Pairs中的两个基因都在GID中的Pairs并将Pairs中的GeneID转化为索引.
  Pairs2Index<-na.exclude(cbind(match(Pairs[,1],GID),match(Pairs[,2],GID)));
  LN<-dim(DataN)[2]
  LT<-dim(DataT)[2]
  m<-dim(Pairs)[1]
  Pair.News<-matrix(nrow=m,ncol=5)
  
  if (m<=Step)
  {
    TempN<-apply((DataN[Pairs2Index[,1],]-DataN[Pairs2Index[,2],])>0,1,sum)
    TempT<-apply((DataT[Pairs2Index[,1],]-DataT[Pairs2Index[,2],])>0,1,sum)
    Pvalue<-phyper(TempN-1, TempN+TempT, LN+LT-TempN-TempT, LN, lower.tail=F)
    TempPair=cbind(Pairs2Index,TempN/LN,TempT/LT,Pvalue)
    Pair.News[1:m,]=TempPair
  }
  
  else
  {
    count<-0
    for (i in seq(Step,m,Step)) 
    {
      TempPair<-Pairs2Index[(count+1):i,]
      TempN<-apply((DataN[TempPair[,1],]-DataN[TempPair[,2],])>0,1,sum)
      TempT<-apply((DataT[TempPair[,1],]-DataT[TempPair[,2],])>0,1,sum)
      Pvalue<-phyper(TempN-1, TempN+TempT, LN+LT-TempN-TempT, LN, lower.tail=F)
      TempPair1<-cbind(TempPair,TempN/LN,TempT/LT,Pvalue)
      Pair.News[(count+1):i,]=TempPair1
      count=i;
    }
    TempPair<-Pairs2Index[(count+1):m,]
    TempN<-apply((DataN[TempPair[,1],]-DataN[TempPair[,2],])>0,1,sum)
    TempT<-apply((DataT[TempPair[,1],]-DataT[TempPair[,2],])>0,1,sum)
    Pvalue<-phyper(TempN-1, TempN+TempT, LN+LT-TempN-TempT, LN, lower.tail=F)
    TempPair1<-cbind(TempPair,TempN/LN,TempT/LT,Pvalue)
    Pair.News[(count+1):m,]=TempPair1
  }
  qValues= p.adjust(Pair.News[,5],method="BH")
  Pair.News<-cbind(Pair.News,qValues)
  Pair.News<-Pair.News[which(Pair.News[,5]<=0.05),]
  #save(Pair.News,file=paste(savename,'.Rdata'))
  Pair.News=as.data.frame(Pair.News) #这行是我加的，允许稳定对用miRNA名称表示，而非代号
  Pair.News[,1]=GID[Pair.News[,1]]
  Pair.News[,2]=GID[Pair.News[,2]]
  Pair.News
}
#逆转对解读
if(0){
    第三列：TempN/LN表示将正常样本判为正常的概率
    第四列：TempT/LT表示将疾病样本判为正常的概率
    故：第三例越大，第四列越小，标记越好
}

#函数调用
dataT=exp13
rpair=Reverse_pair(dataN,dataT,1:2550,spair,nrow(spair),' ')

5.挑选标记

View(rpair[order(-rpair[,3],rpair[,4]),])
#观察发现乳腺癌前135的逆转对完完全全将该癌症与正常分开
#排除干扰取top前200作为标记
 T200=rpair[order(-rpair[,3],rpair[,4]),][1:200,1:2]

6.验证标记

#验证函数
fun<-function(exp,tag,miRNA,miRNA1){
  index<-na.exclude(cbind(match(miRNA[tag[,1]],miRNA1),match(miRNA[tag[,2]],miRNA1)));
  Temp=apply((exp[index[,1],]-exp[index[,2],])>0,1,sum) #判为正常的数量
  Temp/ncol(exp)  #各个标记将疾病或正常样本判为正常的概率
}
#若exp选择的是疾病样本，则概率越低，标记越好
#若exp选择的是正常样本，则概率越高，标记越好

#选取含有该癌症的数据集
data=getGEO('GSE110651',destdir='.')[[1]]
exp1=exprs(data) #表达谱矩阵
miRNA1=data@featureData@data$miRNA_ID_LIST #ID注释
exp1=exp[-grep(',',miRNA1),]
miRNA1=miRNA[-grep(',',miRNA1)]

fun(exp1,T200,miRNA,miRNA1)

7.最后结果

这top200的逆转对中存在该癌症特异的逆转对，能将正常与癌症区分开来