共突变通路识别(R语言实现)

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

一.实验内容

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


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个基因

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

四.实验代码

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
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 215,634评论 6 497
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,951评论 3 391
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 161,427评论 0 351
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,770评论 1 290
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,835评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,799评论 1 294
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,768评论 3 416
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,544评论 0 271
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,979评论 1 308
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,271评论 2 331
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,427评论 1 345
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,121评论 5 340
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,756评论 3 324
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,375评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,579评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,410评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,315评论 2 352

推荐阅读更多精彩内容