该博客介绍了如何通过斯皮尔曼和皮尔森相关性计算lncRNA与蛋白编码基因的共表达关系。首先，对原始数据进行预处理，然后分别计算上调和下调lncRNA的两种相关系数矩阵。通过对每行取最大值，筛选出共表达基因。结果显示，两种方法选出的相同蛋白编码基因有5个，共表达基因总数分别为92（斯皮尔曼）和103（皮尔森）。综合两种方法，共有174个不同的蛋白编码基因被选中。
题目回顾
计算与差异表达lncRNA共表达的蛋白编码基因（同时用斯皮尔曼相关及皮尔森相关），比较两种方法得到的显著共表达蛋白编码基因列表的差异。
在计算共表达蛋白编码基因之前，我们已经得到了差异表达的lncRNA和蛋白编码基因，并且已经区分了上调基因和下调基因。
相关操作解释及思路见识别差异表达蛋白编码基因
以下是上述处理用到的代码，如果已经完成，已经可直接略过
（碎碎念：我也想把这冗长的代码折叠起来，可是我不会）
##---读入文件---##
setwd("D:/QQ/1903786521/FileRecv")
pro.count <- read.table("pro_counts.txt",header = TRUE,row.names = 1,stringsAsFactors = F)
#> View(pro.count)
pro.fpkm <- read.table("protein_coding.fpkm.landscape.txt",header = TRUE,row.names = 1,stringsAsFactors = F)
lnc.fpkm <- read.table("lncRNA.fpkm.landscape.txt",header = TRUE,row.names = 1,stringsAsFactors = F)
##---pro预处理---##
pro.rest <- data.frame()
for(i in 1:nrow(pro.fpkm)){
if(all(pro.fpkm[i,]>0)){
pro.rest <- rbind(pro.rest,pro.fpkm[i,])
}
}
pro.log <- log2(pro.rest+0.05)
##---lnc预处理---##
lnc.rest <- data.frame()
for(i in 1:nrow(lnc.fpkm)){
if(all(lnc.fpkm[i,]>0)){
lnc.rest <- rbind(lnc.rest,lnc.fpkm[i,])
}
}
lnc.log <- log2(lnc.rest+0.05)
##---pro t检验---##
pro.pvalues <- vector()
for(i in 1:nrow(pro.log)){
pvalue <- t.test(pro.log[i,grep("ESC",colnames(pro.log))],
pro.log[i,grep("CM",colnames(pro.log))])$p.value
pro.pvalues <- c(pro.pvalues,pvalue)
}
pro.fdr <- p.adjust(pro.pvalues)
##---lnc t检验---##
lnc.pvalues <- vector()
for(i in 1:nrow(lnc.log)){
pvalue <- t.test(lnc.log[i,grep("ESC",colnames(lnc.log))],
lnc.log[i,grep("CM",colnames(pro.log))])$p.value
lnc.pvalues <- c(lnc.pvalues,pvalue)
}
lnc.fdr <- p.adjust(lnc.pvalues)
##---pro fold change---##
pro.mean_esc <- apply(pro.rest[,grep("ESC",colnames(pro.rest))], 1, mean)
pro.mean_cm <- apply(pro.rest[,grep("CM",colnames(pro.rest))], 1, mean)
fc <- pro.mean_cm / pro.mean_esc
pro.fc.log <- log2(fc)
##---lnc fold change---##
lnc.mean_esc <- apply(lnc.rest[,grep("ESC",colnames(lnc.rest))], 1, mean)
lnc.mean_cm <- apply(lnc.rest[,grep("CM",colnames(lnc.rest))], 1, mean)
fc <- lnc.mean_cm / lnc.mean_esc
lnc.fc.log <- log2(fc)
##---pro 合并---##
pro.result <- data.frame(esc=pro.mean_esc,cm=pro.mean_cm,log2fc=pro.fc.log,
pvalues=pro.pvalues,fdr=pro.fdr)
##---lnc 合并---##
lnc.result <- data.frame(esc=lnc.mean_esc,cm=lnc.mean_cm,log2fc=lnc.fc.log,
pvalues=lnc.pvalues,fdr=lnc.fdr)
##---pro tab---##
pro.result$lable <- ifelse(pro.result$pvalues > 0.01,"mid",
ifelse(pro.result$log2fc > 2,"up",
ifelse(pro.result$log2fc < -2,"down","mid")))
pro.up <- pro.result[grep("up",pro.result$lable),]
pro.down <- pro.result[grep("down",pro.result$lable),]
#ggplot(pro.result, aes(x=pro.fc.log, y=-log10(fdr)))+geom_point(aes(color=lable))+labs(title="Volcano")
##---lnc tab---##
lnc.result$lable <- ifelse(lnc.result$pvalues > 0.01,"mid",
ifelse(lnc.result$log2fc > 2,"up",
ifelse(lnc.result$log2fc < -2,"down","mid")))
lnc.up <- lnc.result[grep("up",lnc.result$lable),]
lnc.down <- lnc.result[grep("down",lnc.result$lable),]
lnc.up.fpkm <- lnc.fpkm[rownames(lnc.up),]
pro.up.fpkm <- pro.fpkm[rownames(pro.up),]
lnc.down.fpkm <- lnc.fpkm[rownames(lnc.down),]
pro.down.fpkm <- pro.fpkm[rownames(pro.down),]
spearman相关计算共表达蛋白编码基因
求相关系数矩阵
对于一个确定的lncRNA和一个确定的蛋白编码基因，能取出两行表达值，用这两行表达值可以求出两者的相关系数，所以每个lncRNA和每个蛋白编码基因之间都能得到一个相关系数，由此可以构建一个相关系数的矩阵，
比如：行名lncRNA,列名蛋白编码基因~
（本喵急着赶作业，就先不画图了，下次一定补上）
##---上调lnc spearman相关系数---##
up.spearman <- matrix(ncol=nrow(pro.up))
for(i in 1:nrow(lnc.up)){
a.lnc <- t(as.matrix(lnc.up.fpkm[i,]))
sprm <- vector()
for(j in 1:nrow(pro.up)){
a.pro <- t(as.matrix(pro.up.fpkm[j,]))
sprm <- cbind(sprm,cor(a.pro,a.lnc,method="spearman"))
}
up.spearman <- rbind(up.spearman,sprm)
}
up.spearman <- up.spearman[-1,]
rownames(up.spearman) <- rownames(lnc.up)
colnames(up.spearman) <- rownames(pro.up)
##---下调lnc spearman相关系数---##
down.spearman <- matrix(ncol=nrow(pro.down))
for(i in 1:nrow(lnc.down)){
a.lnc <- t(as.matrix(lnc.down.fpkm[i,]))
sprm <- vector()
for(j in 1:nrow(pro.down)){
a.pro <- t(as.matrix(pro.down.fpkm[j,]))
sprm <- cbind(sprm,cor(a.pro,a.lnc,method="spearman"))
}
down.spearman <- rbind(down.spearman,sprm)
}
down.spearman <- down.spearman[-1,]
rownames(down.spearman) <- rownames(lnc.down)
colnames(down.spearman) <- rownames(pro.down)
以上代码的up.spearman和down.spearman，分别是上调和下调的相关系数矩阵。
选共表达基因
然后选共表达基因，对于一个确定的lncRNA,它有一整行与蛋白编码基因之间的相关系数，那么应该选哪个呢？
这里采用的做法是每行取最大值，然后每个lncRNA就都找到了自己配对的蛋白编码基因
（当然这里有个bug：因为有可能某个蛋白编码基因既是lncRNA1这行的最大值，又是lncRNA2那行的最大值，那么这个蛋白编码基因就会被选中两次，我没觉得这有很大的问题，所以暂时不解决这个问题了）
##---up lnc spearman 筛选共表达基因---##
up.names <- colnames(up.spearman)
cov.up.spearman <- data.frame(lncName=rownames(up.spearman))
for(i in 1:nrow(up.spearman)){
maxx <- max(up.spearman[i,])
cov.up.spearman$spearman[i] <- maxx
proName <- up.names[grep(maxx,up.spearman[i,])]
cov.up.spearman$proName[i] <- proName
}
##---down lnc spearman 筛选共表基因---##
down.names <- colnames(down.spearman)
cov.down.spearman <- data.frame(lncName=rownames(down.spearman))
for(i in 1:nrow(down.spearman)){
maxx <- max(down.spearman[i,])
cov.down.spearman$spearman[i] <- maxx
proName <- down.names[grep(maxx,down.spearman[i,])]
cov.down.spearman$proName[i] <- proName
}
以上代码的cov.up.spearman和cov.down.spearman就是最后的筛选结果，三列分别是lncRNA、相关系数、蛋白编码基因，大概长这样：
pearson相关计算共表达蛋白编码基因
pearson和spearman的过程基本一样，改了个方法而已，就不做过多解释了。
##---上调lnc pearson相关系数---##
up.pearson <- matrix(ncol=nrow(pro.up))
for(i in 1:nrow(lnc.up)){
a.lnc <- t(as.matrix(lnc.up.fpkm[i,]))
prs <- vector()
for(j in 1:nrow(pro.up)){
a.pro <- t(as.matrix(pro.up.fpkm[j,]))
prs <- cbind(prs,cor(a.pro,a.lnc,method="pearson"))
}
up.pearson <- rbind(up.pearson,prs)
}
up.pearson <- up.pearson[-1,]
rownames(up.pearson) <- rownames(lnc.up)
colnames(up.pearson) <- rownames(pro.up)
##---下调lnc pearson相关系数---##
down.pearson <- matrix(ncol=nrow(pro.down))
for(i in 1:nrow(lnc.down)){
a.lnc <- t(as.matrix(lnc.down.fpkm[i,]))
prs <- vector()
for(j in 1:nrow(pro.down)){
a.pro <- t(as.matrix(pro.down.fpkm[j,]))
prs <- cbind(prs,cor(a.pro,a.lnc,method="pearson"))
}
down.pearson <- rbind(down.pearson,prs)
}
down.pearson <- down.pearson[-1,]
rownames(down.pearson) <- rownames(lnc.down)
colnames(down.pearson) <- rownames(pro.down)
##---up lnc pearson 筛选共表达基因---##
up.names <- colnames(up.pearson)
cov.up.pearson <- data.frame(lncName=rownames(up.pearson))
for(i in 1:nrow(up.pearson)){
maxx <- max(up.pearson[i,])
cov.up.pearson$pearson[i] <- maxx
proName <- up.names[grep(maxx,up.pearson[i,])]
cov.up.pearson$proName[i] <- proName
}
##---down lnc pearson 筛选共表基因---##
down.names <- colnames(down.pearson)
cov.down.pearson <- data.frame(lncName=rownames(down.pearson))
for(i in 1:nrow(down.pearson)){
maxx <- max(down.pearson[i,])
cov.down.pearson$pearson[i] <- maxx
proName <- down.names[grep(maxx,down.pearson[i,])]
cov.down.pearson$proName[i] <- proName
}
两种方式筛选出的蛋白编码基因比较
##---比较差异---##
cov.spearman <- rbind(cov.up.spearman,cov.down.spearman)
cov.pearson <- rbind(cov.up.pearson,cov.down.pearson)
cov.df <- data.frame(lncName=cov.spearman$lncName,spearman=cov.spearman$spearman,
proSpearman=cov.spearman$proName,pearson=cov.pearson$pearson,
proPearson=cov.pearson$proName)
以上代码只是将两个方法选出的蛋白编码基因的结果综合了一下，cov.df长这样：
下面稍微做一下比较
##---同一个lncRNA选出了相同的蛋白编码基因---##
length(cov.df[cov.df$proSpearman==cov.df$proPearson,]) #5
##---筛选出来总的蛋白编码基因数---##
nrow(cov.df) #103
length(unique(cov.df$spearman)) #92
length(unique(cov.df$pearson)) #103
length(unique(c(cov.df$proSpearman,cov.df$proPearson))) #174
发现同一个lncRNA选出了相同的蛋白编码基因的情况发生了五次，上调下调lncRNA共103个，用spearman选出的共表达蛋白编码基因92个，用pearson选出的共表达蛋白编码基因103个，把两种筛选结果综合一下剩下174个，说明其中存在92+103-174=21个两种方法都选中的蛋白编码基因。
1765
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
