R语言批量计算一个基因的表达相关基因—科研工具箱

请关注公众号【叨客学习资料】 在使用网站的过程中有疑问,请来公众号进行反馈哦

代码不一定是最优化的,根据实际情况修改使用

rm(list = ls())
library(dplyr)
##读取基因ID,Ensembl ID注释文件
idmap <- read.delim(\"G:\\\\GTeX TCGA/gencode.v23.annotation.gene.probemap\",as.is=T)
#读取表达矩阵
rt <- data.table::fread(\"G:/GTeX TCGA/tcga_RSEM_gene_tpm/split/LUAD.csv\",data.table = F)
rt <- as.matrix(rt) #变为matrix
rownames(rt) <- rt[,1]
exp <-as.data.frame( rt[,2:ncol(rt)]) 
exp <- apply(exp, 2, as.numeric)
rownames(exp) <- rownames(rt) 
exp <- exp[,substr(colnames(exp),14,15)<10]
#将表达为0的样本赋值为NA
exp[exp== -9.9658]<-NA
# 求相关性
target <- \"FBXO22\"
#ID转换
gene <- idmap$id[which(idmap$gene==target)]
outTab=data.frame()
for (i in rownames(exp)){
  x=exp[gene,]
  y=exp[i,]
  cordata<- data.frame(x=as.numeric(x),y=as.numeric(y)) %>% na.omit()
  if (sd(cordata$y)==0|nrow(cordata)<10){
    next
  }
  corT=cor.test(cordata$x,cordata$y)
  cor=corT$estimate
  cor=round(cor,3)
  pvalue=as.numeric(corT$p.value)
  if(pvalue<0.001){
    pvalue=format(pvalue,scientific = TRUE)
  }else{
    pvalue=round(pvalue,3)
  }

    outTab=rbind(outTab,cbind(lncRNA=gene,mRNA=i,cor,pvalue))
}
##合并 需要,换一下列名
colnames(outTab)[2] <- \"Ensembl\"
outTab$Ensembl <- substr(outTab$Ensembl,1,15)
##另一个整理好的注释文件,包含RNAleixing
biomart_table <- readRDS(paste0(\"E:/apps/9. IDconvert/app_data/biomart_table_\",\"homo\",\".Rds\"))
outTab <- as.data.frame(outTab)
#ID转化
converted<-biomart_table %>%
  filter_all(., any_vars(. %in% outTab[,2])) %>%
  distinct(Ensembl, .keep_all = TRUE) %>%
  arrange(Ensembl)
conv.data<-merge(outTab,biomart_table,by=\"Ensembl\")[-1]
##保存数据
write.csv(conv.data,paste0(choose.files(),\"/\",target,\"_cor.csv\"))

将上述代码包装成函数:

corr_tcga <- function(tcga=\"LUAD\",target=\"TP53\"){
  rt <- data.table::fread(paste0(\"G:/GTeX TCGA/tcga_RSEM_gene_tpm/split/\",tcga,\".csv\"),data.table = F)
  rt <- as.matrix(rt) #变为matrix
  rownames(rt) <- rt[,1]
  exp <-as.data.frame( rt[,2:ncol(rt)]) 
  exp <- apply(exp, 2, as.numeric)
  rownames(exp) <- rownames(rt) 
  exp <- exp[,substr(colnames(exp),14,15)<10]
  exp[exp== -9.9658]<-NA
  # 求相关性
  gene <- idmap$id[which(idmap$gene==target)]
  outTab=data.frame()
  for (i in rownames(exp)){
    x=exp[gene,]
    y=exp[i,]
    cordata<- data.frame(x=as.numeric(x),y=as.numeric(y)) %>% na.omit()
    if (sd(cordata$y)==0|nrow(cordata)<10){
      next
    }
    corT=cor.test(cordata$x,cordata$y)
    cor=corT$estimate
    cor=round(cor,3)
    pvalue=as.numeric(corT$p.value)
    if(pvalue<0.001){
      pvalue=format(pvalue,scientific = TRUE)
    }else{
      pvalue=round(pvalue,3)
    }
    
    outTab=rbind(outTab,cbind(Target=target,Ensembl=i,cor,pvalue))
  }
  outTab$Ensembl <- substr(outTab$Ensembl,1,15)
  biomart_table <- readRDS(paste0(\"E:/apps/9. IDconvert/app_data/biomart_table_\",\"homo\",\".Rds\"))
  outTab <- as.data.frame(outTab)
  converted<-biomart_table %>%
    filter_all(., any_vars(. %in% outTab[,2])) %>%
    distinct(Ensembl, .keep_all = TRUE) %>%
    arrange(Ensembl)
  conv.data<-merge(outTab,biomart_table,by=\"Ensembl\")[-1]
}
conv.data <- corr_tcga(\"LUSC\",\"FBXO22\")
write.csv(conv.data,choose.files())
图片[1]-R语言批量计算一个基因的表达相关基因—科研工具箱-叨客学习资料网
© 版权声明
THE END
喜欢就支持一下吧
点赞0 分享
评论 抢沙发
头像
请输入有效评论哦,肆意灌水或者乱打评论是不会通过的,会影响您评论后获得资源哦~~
提交
头像

昵称

取消
昵称表情

    暂无评论内容