ID转换以及LncRNA和mRNA相关性分析
2022-05-09 本文已影响0人
monkey_study
基因注释
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) #查看数据维度
最后，结果还是有点小问题（😵），不过大概流程还可以哈! 共同学习进步，加油！
