写在前面：我认为此次的任务重点在理清楚实验原理和步骤，代码部分比较基础
也可用来研究共突变基因

一．实验内容

共突变通路识别
具体原理步骤如下

image.png
image.png
image.png
image.png

二．实验目的

（1）掌握通路富集方法
（2）掌握共突变通路计算方法

三．实验数据工具及步骤

实验数据

1.突变基因列表var_exp.csv

image.png

2.通路名字pathway_Name.csv

image.png
3.通路涉及到的基因Path_Gene.csv
image.png

4.通路基因对应矩阵path_info_mtx.csv

image.png

实验步骤（version:R.3.6.3）

1.获得每个样本的突变基因列表，来自数据1，共986个样本，18062个基因

2. 突变基因列表与KEGG通路中所涉及基因作交集
3. 利用超几何富集模型，识别每个样本富集到显著富集突变基因的通路
   1-phyper(k-1,M,N-M,n)
4. 利用超几何模型识别共突变通路
5. 得到共突变通路对。

四．实验代码

setwd("E:\\实验\\转录组学\\实验四")
exp<-read.table("varexp.txt",header=T) #18062*986,突变表达谱
gene<-read.table("path_Gene.csv",as.is=T)#6200行，通路里的所有基因
pathway<-read.csv("pathway_Name.csv",as.is=T,header=F)#247行，涉及到的所有通路
path_gene<-read.csv("path_info_mtx.csv",as.is=T,header=F)#247*6200，每个通路富集的基因情况
#加载要用的包
install.packages("plyr")
library(plyr)
install.packages("dplyr")
library(dplyr)
#求基因交集表达谱，我们只研究既在突变表达谱里的基因又在通路里的基因的表达情况，所以对gene和exp的行取交集，得到5659（基因）*986（样本）交集表达谱
a<-as.numeric(rownames(exp))
jiao<-intersect(a,gene[,1])#5659，取通路里的所有基因和突变表达谱基因的交集
mu_exp<-exp %>% filter(rownames(exp)%in% jiao)  #5659*986，交集表达谱


#分别对突变表达谱和path_gene表达谱统计表达情况
mu<-apply(mu_exp,2,sum)#按列，求和，统计个数，统计每个样本突变的基因数
fuji<-apply(path_gene,1,sum)#按行，统计每个通路的基因数

#把两个矩阵 突变的数据和富集到的基因数据（=1）的提出来
muExp1<-as.data.frame(list())#建立两个空数据框
pathExp1<-as.data.frame(list())

for (i in 1:length(colnames(mu_exp)))
{
muExp<-as.data.frame(t(rownames(mu_exp)[which(mu_exp[,i]==1)]))
muExp1<-rbind.fill(muExp,muExp1)
}
#986*1419，突变表达谱

pathExp1<-as.data.frame(list())
for (j in 1:length(rownames(path_gene)))
{
pathExp<-as.data.frame(t(gene[which(path_gene[j,]==1),1]))
pathExp1<-rbind.fill(pathExp,pathExp1)
}
#247*1146，每个通路富集的基因情况

#取交集，也就是得到每个通路在每个样本突变基因数矩阵，247*986
mu_fu<-matrix(NA,length((pathway)[,1]),length(colnames(exp)))

for (i in 1:length(rownames(pathExp1)))
{
for (j in 1:length(rownames(muExp1)))
{
mu_fu[i,j]<-length(na.omit(intersect(muExp1[,j],pathExp1[,i])))
}

}
#247*986


#累积超几何分布，算p值,得到p值矩阵，看显著性
pvalue<-matrix(NA,length((pathway)[,1]),length(colnames(exp)))
for (i in 1:length(fuji))
{
n<-fuji[i]
for (j in 1:length(mu))
{
M<-mu[j]
k<-mu_fu[i,j]
pvalue[i,j]<-1-phyper(k-1,M,5659-M,n)
}
}
#显著的p找出来，0为<0.05,1为>0.05，识别每个样本富集到显著富集突变基因的通路

aa<-(pvalue<0.05)
for(i  in 1:length((pathway)[,1]))
{
for(j in 1:length(colnames(exp)))
{
if(aa[i,j]==TRUE)
aa[i,j]=1
else
aa[i,j]=0
}
}
rownames(aa)<-pathway[,1]
colnames(aa)<-colnames(exp)
write.table(aa,"pathway_sample.txt",sep="\t")#写出

#aa
     0      1 
230536  13006
#p值较正，显著的p找出来，0为<0.05,1为>0.05
qvalue<-p.adjust(pvalue,method='fdr')
write.table(qvalue,"qvalue.txt",sep="\t")
bb<-(qvalue<0.05)
bb<-matrix(bb,247,986,byrow=T)
for(i  in 1:length((pathway)[,1]))
{
for(j in 1:length(colnames(exp)))
{
if(bb[i,j]==TRUE)
bb[i,j]=1
else
bb[i,j]=0
}
}
#bb
 0     1
237097   6445
rownames(bb)<-pathway[,1]
colnames(bb)<-colnames(exp)

#共突变（用的是没有矫正的p）
#构建一个三列矩阵，第一二列为通路名字，组合数为30381行，第三列为p值
pathway_sample<-aa
path_path<-matrix(0,choose(length(rownames(pathway_sample)),2) ,3) #30381*3
path_path[,1]<-t(combn(rownames(pathway_sample),2))[,1]  #组合数
path_path[,2]<-t(combn(rownames(pathway_sample),2))[,2]
colnames(path_path)<-c("pathway1","pathway2","p")
path_path[,3]<-0
path_path[,3]<-as.numeric(path_path[,3])
path_sa<-apply(pathway_sample,1,sum)#按行求和,247个值，代表每个通路富集的基因数

#按通路找=1的样本,找交集，把通路上显著富集的样本名列出，247行
path_s<-list() #构建空list
for (i in 1:length(rownames(pathway_sample)))
{
path_s[i]<-list(which(pathway_sample[i,]==1))
}

path_s[[1]]结构如下

image.png

qq<-matrix(NA,247,247) #建立矩阵，把对角线设为0
for (i in 1:247)
{
for (j in 1:247)
{   
qq[i,j]<-length(intersect(path_s[[i]],path_s[[j]]))
if(i==j)
qq[i,j]<-0
}}

K<-qq[upper.tri(qq)] #取上三角，通路交集的得到组合数30381个值（两个通路之间共同富集的基因数）

n<-t(combn(path_sa,2))[,1]
M<-t(combn(path_sa,2))[,2]

#累积超几何分布，找p值显著的通路对
path_path[,3]<-1-phyper(K-1,M,30381-M,n)
write.table(path_path,"path_path.txt",sep="\t")
path_pair<-path_path[which(path_path[,3]<0.05),]
write.table(path_pair,"path_pair.txt",sep="\t") #20246
#####################################################################
同理，矫正后的p值
setwd("E:\\实验\\转录组学\\实验四\\adjust")
qvalue<-p.adjust(pvalue,method='fdr')
write.table(qvalue,"qvalue.txt",sep="\t")
bb<-(qvalue<0.05)
bb<-matrix(bb,247,986,byrow=T)
for(i  in 1:length((pathway)[,1]))
{
for(j in 1:length(colnames(exp)))
{
if(bb[i,j]==TRUE)
bb[i,j]=1
else
bb[i,j]=0
}
}
pathway_sample<-bb
后面代码一样
write.table(path_pair1,"path_pair1.txt",sep="\t")   #1233

实验结果

image.png
pathway_sample.txt
path_path.txt
qvalue.txt