R语言跑生信脚本出错，求助各位大佬

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

outTab=data.frame()
> for(i in gene1){
+   for(j in gene2){
+     x=as.numeric(data[i,])
+     y=as.numeric(data[j,])
+     corT=cor.test(x, y, method = 'spearman')
+     cor=corT$estimate
+     pvalue=corT$p.value
+     if((abs(cor)>corFilter) & (pvalue<pvalueFilter)){
+       outTab=rbind(outTab, cbind(Gene1=i, Gene2=j, cor, pvalue))
+       #绘制相关性图形
+       df1=as.data.frame(cbind(x,y))
+       p1=ggplot(df1, aes(x, y)) +
+         xlab(i)+ylab(j)+
+         geom_point()+ geom_smooth(method="lm",formula = y ~ x) + theme_bw()+
+         stat_cor(method = 'spearman', aes(x =x, y =y))
+       p2=ggMarginal(p1, type = "density", xparams = list(fill = "orange"),yparams = list(fill = "blue"))
+       #输出相关性图形
+       pdf(file=paste0(i, "_", j, ".pdf"),width=5,height=4.8)
+       print(p2)
+       dev.off()
+     }
+   }
+ }
Error in data[j, ] : subscript out of bounds

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：R语言 estimate spearman numeric pearman

data的维数是多少，gene1和gene2分别是多少，看一看。一般是gene2的值超出了data的维数

00feng00 发表于 2022-8-19 10:41
data的维数是多少，gene1和gene2分别是多少，看一看。一般是gene2的值超出了data的维数

您好，这是原数据和原代码，实在看不出哪有问题，劳烦您看一下，谢谢您！

#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("limma")

#install.packages("ggplot2")
#install.packages("ggpubr")
#install.packages("ggExtra")


#引用包
library(limma)
library(ggplot2)
library(ggpubr)
library(ggExtra)

corFilter=0.4            #相关系数的过滤标准
pvalueFilter=0.001       #相关性检验pvalue的过滤标准
setwd("C:\\11.cor")     #设置工作目录

#读取表达数据文件,并对数据进行处理
rt=read.table("m6aGeneExp.txt", header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)

#删掉对照组样品
group=sapply(strsplit(colnames(data),"\\_"), "[", 2)
data=data[,group=="treat"]

#读取基因列表文件
geneRT=read.table("gene.txt", header=T, sep="\t", check.names=F, row.names=1)
gene1=row.names(geneRT)[geneRT[,1]=="erasers"]      #提取消码器相关的基因
gene2=row.names(geneRT)[geneRT[,1]=="writers"]      #提取编码器相关的基因

#相关性检验
outTab=data.frame()
for(i in gene1){
    for(j in gene2){
        x=as.numeric(data[i,])
        y=as.numeric(data[j,])
        corT=cor.test(x, y, method = 'spearman')
        cor=corT$estimate
        pvalue=corT$p.value
        if((abs(cor)>corFilter) & (pvalue<pvalueFilter)){
            outTab=rbind(outTab, cbind(Gene1=i, Gene2=j, cor, pvalue))
            #绘制相关性图形
            df1=as.data.frame(cbind(x,y))
            p1=ggplot(df1, aes(x, y)) +
                xlab(i)+ylab(j)+
                geom_point()+ geom_smooth(method="lm",formula = y ~ x) + theme_bw()+
                stat_cor(method = 'spearman', aes(x =x, y =y))
            p2=ggMarginal(p1, type = "density", xparams = list(fill = "orange"),yparams = list(fill = "blue"))
            #输出相关性图形
            pdf(file=paste0(i, "_", j, ".pdf"),width=5,height=4.8)
            print(p2)
            dev.off()
        }
    }
}

#输出相关性检验的结果
write.table(file="corResult.txt",outTab,sep="\t",quote=F,row.names=F)

数据看不了劳烦您加一下我的微信，我给您发（md798452941）