ID转换以及LncRNA和mRNA相关性分析

基因注释

UCSC Xena的基因注释文件gencode.v22.annotation.gene.probeMap，下载方式，打开所研究癌种的count/FPKM页面，找到probeMap，下载即可

#读入下载于UCSC Xena的基因注释文件，怎么下载呢，上边说啦
probeMap <- read.table("TCGAdata/gencode.v22.annotation.gene.probeMap",sep = "\t" , header = T)  #加载TCGA官方基因注释文件probeMap
dim(probeMap)  #60483个基因探针对应文件     6列
probeMap[1:4,1:6]  #6列分别是 基因id，基因名，染色体位置，染色体起，止以及正反链
                                 id         gene chrom chromStart chromEnd strand
1 ENSG00000223972.5      DDX11L1  chr1      11869    14409      +
2 ENSG00000227232.5       WASH7P  chr1      14404    29570      -
3 ENSG00000278267.1    MIR6859-3  chr1      17369    17436      -
4 ENSG00000243485.3 RP11-34P13.3  chr1      29554    31109      +
nrDEG[1:4,1:4]  #之前我已用limma voom筛选到的LncRNA差异基因，行名为 Ensembl 基因 ID（带有版本号），第一次尝试奥，都是自己摸索的。
rownames(nrDEG) <- probeMap[rownames(nrDEG)%in%probeMap$id,2]  #用注释文件将nrDEG行名（Ensembl 基因 ID ）转换为gene symbol（probeMap第二列） 
head（rownames(nrDEG) ）#查看nrDEG行名前6个，判断转换是否成功
save(nrDEG，file = "anno_nrdeg.Rdata")   #保存为Rdata文件
load （"anno_nrdeg.Rdata"） #读入Rdata文件，即载入nrDEG，方便后续操作，也省得再次运行上述代码。
###读入坏死性凋亡相关基因名并命名为hs_hypotosis
hs_hypotosis <- data.table::fread("TCGAdata/gene.txt",header = T)
#为了进行后续的提取坏死性凋亡相关基因，我们对全部的表达矩阵进行ID转换一下
all_count_stad[1:4,1:4] #查看表达矩阵
                  TCGA-D7-5577-01A TCGA-D7-6818-01A TCGA-BR-4280-01A
ENSG00000242268.2                0                0                0
ENSG00000259041.1                0                0                0
ENSG00000270112.3                0                2                0
ENSG00000278814.1                0                0                0
                  TCGA-D7-8572-01A
ENSG00000242268.2                4
ENSG00000259041.1                0
ENSG00000270112.3                5
ENSG00000278814.1                0
#对表达矩阵进行ID转换
all_count_stad$symbol <- probeMap[probeMap$id%in%rownames(all_count_stad),2]
###利用limma包中的 avereps函数对重复基因取平均最大值并去重，转换为数据框赋值给TCGA_gset 
library(limma)
TCGA_gset = as.data.frame(avereps(all_count_stad[,-ncol(all_count_stad)],ID = all_count_stad$symbol) )   #经过数据检查，发现存在多个重复的基因名，这句代码是去重复的（重复原因：许多测序序列会比对到转录组的相同位置）
#提取坏死性凋亡相关基因
hs_relate_gene <- TCGA_gset[rownames(TCGA_gset)%in%hs_hypotosis$Genes,]

• 相关性分析，输入数据：LncRNAcount矩阵（行为样本，列为基因），mRNA count矩阵 （行为样本，列为基因），本例中分别为lncRNA_expr，hs_relate_mRNA
• 输入数据的准备

lncRNA_anno[1:4,1:4]  #注释后的DElncRNA
##lncRNA 表达矩阵准备  行 样本，列基因
lncRNA_expr <- TCGA_gset[rownames(TCGA_gset)%in%rownames(nrDEG),]
lncRNA_expr <- t(lncRNA_expr)
str(lncRNA_expr)
#坏死性凋亡相关基因矩阵准备  行样本，列基因
hs_relate_mRNA <- t(hs_relate_gene)
#查看数据
dim(hs_relate_mRNA)  # 547  10
dim(lncRNA_expr)  #547 1076
hs_relate_mRNA <- hs_relate_mRNA[,-4]  #检查数据后发现，第四列数据标准差为0，会导致无法计算相关性，所以这一步去掉它，原因呢，就推荐大家自行百度一下统计相关性计算部分

正式开始

cor_matrix<- Hmisc::rcorr(lncRNA_expr,hs_relate_mRNA)
##自定义函数将结果转换为四列，行，列，cor ，pvalue
flattenCorrMatrix <- function(cormat, pmat) {
  ut <- upper.tri(cormat)
  data.frame(
    lncRNA = rownames(cormat)[row(cormat)[ut]],
     gene= rownames(cormat)[col(cormat)[ut]],
    cor  =(cormat)[ut],
    p = pmat[ut]
  )
}

cor_lncRNA<- flattenCorrMatrix(cor_matrix$r,cor_matrix$P)
dim(cor_lncRNA)
cor_lncRNA$regulation <- ifelse(cor_lncRNA$cor>0,"positive","negative") #正相关为positive，负相关为negative
library(dplyr)
hs_relate_lncRNA <- cor_lncRNA[intersect(which(cor_lncRNA$p<0.001),which(cor_lncRNA$cor>0.4)),]#筛选条件 p<0.001 and cor>0.4
dim(hs_relate_lncRNA)  #查看数据维度
![image.png](https://img.haomeiwen.com/i27731909/b2f4f97db3a98b30.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

最后，结果还是有点小问题（😵），不过大概流程还可以哈! 共同学习进步，加油！