2022-09-28 看看基因之间的相关性

2022-09-28  本文已影响0人  学习生信的小兔子
load(file = "data/TCGA_BRCA_exprSet_plot.Rdata")
data=exprSet[,-c(1,2,3)]

#选取基因 PDCD1
gene="PDCD1"
genedata=as.numeric(data[,"PDCD1"])

#提取基因列表
genelist=colnames(data)

#定义一个数据框
correlation <- data.frame()
for (i in 1:length(genelist)){
  print(i)
  dd=cor.test(genedata,as.numeric(data[,i]),method="spearman")
  correlation[i,1]=gene
  correlation[i,2]=genelist[i]
  correlation[i,3]=dd$estimate
  correlation[i,3]=dd$p.value
}
colnames(correlation) <- c("gene1","gene2","cor","p.value")
p值矫正
correlation$padjust = p.adjust(correlation$p.value,method = "BH")
rm(list = ls())
load(file = "data/TCGA_ggplot.Rdata")

my_comparisons <- list(
  c("LumA", "Normal"),
  c("LumB", "Normal"),
  c("Her2", "Normal"),
  c("Basal", "Normal")
)
library(ggpubr)
ggboxplot(
  exprSet, x = "subgroup", y = "ESR1",
  color = "subgroup", palette = c("#00AFBB", "#E7B800", "#FC4E07","#A687D8", "#89E4A4"),
  add = "jitter"
)+
  stat_compare_means(comparisons = my_comparisons, method = "t.test")


## 改写一下,方便调用
subgroup_plot <- function(gene){
  ggboxplot(
    exprSet, x = "subgroup", y = gene,
    color = "subgroup", palette = c("#00AFBB", "#E7B800", "#FC4E07","#A687D8", "#89E4A4"),
    add = "jitter"
  )+
    stat_compare_means(comparisons = my_comparisons, method = "t.test")
  
  
}
## "OR4F5"    "SAMD11"   "BRCA1"    "ESR1" 

subgroup_plot("SAMD11")

subgroup_plot("BRCA1")

上一篇下一篇

猜你喜欢

热点阅读