单基因批量相关性分析 TCGA基因相关性分析 单基因批量相关性分析的妙用

首先，做这个相关性分析，在这里需要安装几个R包。
https://blog.csdn.net/leianuo123/article/details/102613945

https://mp.weixin.qq.com/s?__biz=MzIyMzA2MTcwMg==&mid=2650733008&idx=1&sn=b66e3fd527f99ddf19dcf6c2501e5be3&chksm=f029aa79c75e236f8951b87e17a51dc6a7dfeb555b983d8dc6d5a8c756b3eacc21829f19dc86&cur_album_id=1336402123646631937&scene=190#rd

使用场景
1.已经确定研究的基因，但是想探索他潜在的功能，可以通过跟这个基因表达最相关的基因来反推他的功能，这种方法在英语中称为guilt of association，协同犯罪。

2.我们的注释方法依赖于TCGA大样本，既然他可以注释基因，那么任何跟肿瘤相关的基因都可以被注释，包括长链非编码RNA

BiocManager::install(c("tidyr","dplyr","ggstatsplot"))
##加载的数据，就是我们常见的TCGA的表达矩阵，只是需要先将其进行转置。
##我们在此用得是自己的一个Test的数据。其格式如下:
————————————————

注意，在此处，基因名一定要作为列名，因为后面需用ggstatsplot进行散点图绘制的时候，其需要用到列名中的基因名字。

该数据，列名就是基因名，行名就是样本名。为了缩短数据处理时间，我们在此只取前500列作为后续的实验分析。

自己构建数据 artificial data

#单基因相关性分析






#随机函数
??runif
as.matrix()

runif(5*500,min = 0.1,max = 10.5)

test3=matrix(runif(5*500,min = 0.1,max = 10.5),byrow = TRUE,
                nrow=5,ncol = 500,
             dimnames = list(c(paste0("sample",seq(1,5,1))),
                             c(paste0("gene",seq(1,500,1)))
                             )
            )
head(test3)[1:3,1:4]


y <- as.numeric(test3[,"gene1"])        ## 本次操作中，选取A1BG基因作为需要分析的单基因，批量求取它与该表达矩阵中其他基因的相关性
head(y)
colnames <- colnames(test3)
cor_data_df <- data.frame(colnames)    ##构建这些基因的数据框

##利用 for循环来求A1BG与其他基因之间的相关性
library（tidyr)
library(dplyr)

for (i in 1:length(colnames))
{
  test <- cor.test(as.numeric(test3[,i]),y,type="spearman")
  cor_data_df[i,2] <- test$estimate
  cor_data_df[i,3] <- test$p.value
}

names(cor_data_df) <- c("Symbol","correlation","pvalue")
head(cor_data_df)

##提取p值小于0.05的数据
cor_data_sig <- cor_data_df %>% 
  filter(pvalue < 0.05) %>% 
  arrange(desc(abs(correlation)))%>% 
  dplyr::slice(1:500)


head(cor_data_sig)


#BiocManager::install('ggside')
library(ggstatsplot)
##利用该包进行绘制相关性的散点图和条形图
ggscatterstats(data =as.data.frame(test3), 
               y = gene1, 
               x = gene2,
               centrality.para = "mean",                              
               margins = "both",                                         
               xfill = "#CC79A7", 
               yfill = "#009E73", 
               marginal.type = "histogram",
               title = "Relationship between A1BG and SH2D3C")

library(clusterProfiler)
#获得基因列表，说明由于我的测试数据质量不是特别好，因此取到的p<0.05的基因只有24条，在此，难以用于GO分析，所以我拿我的全部的基因来做GO分析
library(stringr)
gene <- str_trim(cor_data_sig$symbol,'both')##这是取出的相关性结果好的情况下
 gene <- str_trim(Symbol,'both') ##这里Symbol是所有的基因名
#基因名称转换，返回的是数据框
gene = bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
go <- enrichGO(gene = gene$ENTREZID, OrgDb = "org.Hs.eg.db", ont="all")

##将BP,CC,MF 三个GO类别按照ONTLOGY进行分类，将三个组分绘制在一张图上。
barplot(go, split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free")##条形图
dotplot(go, split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free")##气泡图
View(go)