ixxmu
commented
8 months ago mRNA表达矩阵和lncRNA表达矩阵批量相关性分析识别免疫相关lncRNA by 医学和生信笔记
💡专注R语言在🩺生物医学中的使用
设为“星标”，精彩不错过
之前已经给大家介绍了4种识别免疫相关lncRNA的方法：
今天再给大家介绍一种，也是最简单的一种（去年写的忘记发了🌚🌚）。
我们都知道，如果是想要获取免疫相关的mRNA，可以直接从immport这个网站下载，老牌网站，知名度很高。
直达网址：https://www.immport.org/shared/genelists
下载完成后，我们通过计算lncRNA和免疫相关mRNA的相关性，就可以得到免疫相关lncRNA。
下面进行实操，我们以TCGA-COAD的转录组数据为例。
下载数据
首先下载TCGA-COAD的转录组数据，直接用easyTCGA，1行代码搞定：
library(easyTCGA)
getmrnaexpr("TCGA-COAD")
这样我们就得到了TCGA-COAD的mRNA和lncRNA的表达矩阵，我们使用tpm数据进行演示。这个数据是没有经过log转换的。
load(file = "G:/easyTCGA_test/output_mRNA_lncRNA_expr/TCGA-COAD_lncrna_expr_tpm.rdata")
load(file = "G:/easyTCGA_test/output_mRNA_lncRNA_expr/TCGA-COAD_mrna_expr_tpm.rdata")
然后我们从网站直接下载免疫相关mRNA:
# 免疫相关mRNA是从网站直接下载的
imm.mrna <- read.table("E:/projects/yefd/TCGA_BLCA/files/GeneList.txt",sep="\t",header=T)
head(imm.mrna)
## Symbol ID Name Synonyms
## 1 AZGP1 563 alpha-2-glycoprotein 1, zinc-binding ZA2G|ZAG
## 2 B2M 567 beta-2-microglobulin IMD43
## 3 CALR 811 calreticulin CRT|HEL-S-99n|RO|SSA|cC1qR
## 4 CANX 821 calnexin CNX|IP90|P90
## 5 CD1A 909 CD1a molecule CD1|FCB6|HTA1|R4|T6
## 6 CD1B 910 CD1b molecule CD1|CD1A|R1
## Chromosome Category
## 1 7 Antigen_Processing_and_Presentation
## 2 15 Antigen_Processing_and_Presentation
## 3 19 Antigen_Processing_and_Presentation
## 4 5 Antigen_Processing_and_Presentation
## 5 1 Antigen_Processing_and_Presentation
## 6 1 Antigen_Processing_and_Presentation
我们只要其中的Symbol即可，所以去重一下：
imm.mrna <- unique(imm.mrna$Symbol)
length(imm.mrna)
## [1] 1509
head(imm.mrna)
## [1] "AZGP1" "B2M" "CALR" "CANX" "CD1A" "CD1B"
还剩1509个免疫相关的mRNA。
immune-related mRNA matrix
我们已经有了免疫相关mRNA：imm.mrna，接下来就是提取免疫相关mRNA的表达矩阵。然后使用这个免疫相关mRNA的表达矩阵以及lncRNA的表达矩阵就可以计算相关性了。
首先是提取免疫相关mrna的表达矩阵。
# 取个交集而已
index <- intersect(rownames(mrna_expr_tpm),imm.mrna)
imm_expr <- mrna_expr_tpm[index,]
dim(imm_expr)
## [1] 1098 524
lncRNA matrix
接下来这个lncRNA表达矩阵就有很多处理方式了。
如果你比较偷懒，或者不需要过滤，那你可以直接使用所有的lncRNA和imm_expr中的mRNA批量做相关性分析。
我们这个lncrna_expr_tpm一共有16882个lncRNA，这个数据计算量其实是非常大的！
所以我们也可以先对lncRNA做差异分析，然后再做相关性分析，也可以先做生存分析，也可以差异分析+生存分析一起做。
当然，方法很多，目的只是为了缩小一下目标范围而已，那我做WGCNA行不行？做免疫浸润行不行？做分子分型行不行？都行！
我们这里演示的是差异分析+生存分析，然后再做相关性分析。
先做差异分析，使用easyTCGA包，1行代码解决：
library(easyTCGA)
# 先过滤一下
lncrna_expr <- lncrna_expr_tpm[!apply(lncrna_expr_tpm,1,sd)<0.1,]
deg_res <- diff_analysis(exprset = lncrna_expr,
is_count = F,
logFC_cut = 1,
pvalue_cut = 0.05,
adjpvalue_cut = 0.05
)
## => log2(x+0.1) transform finished
## => Running limma
## => Running wilcoxon test
## => Analysis done.
# 用limma的结果
deg_limma <- deg_res$deg_limma
dim(deg_limma)
## [1] 1421 7
head(deg_limma)
## logFC AveExpr t P.Value adj.P.Val B
## AC092652.1 -3.836206 -2.479661 -35.07337 1.368028e-139 1.847659e-135 307.5981
## VSTM2A-OT1 -1.434307 -3.066762 -30.42429 7.111569e-118 4.802443e-114 258.0771
## LINC00682 -2.770612 -2.762584 -30.13310 1.762372e-116 7.934198e-113 254.8951
## LINC00974 -4.175620 -2.359253 -28.73307 9.966442e-110 3.365169e-106 239.4790
## AC005358.1 -2.800434 -2.632442 -28.56468 6.545520e-109 1.768076e-105 237.6123
## AC025470.2 -2.128380 -2.861899 -27.80674 3.221137e-105 7.250780e-102 229.1800
## genesymbol
## AC092652.1 AC092652.1
## VSTM2A-OT1 VSTM2A-OT1
## LINC00682 LINC00682
## LINC00974 LINC00974
## AC005358.1 AC005358.1
## AC025470.2 AC025470.2
差异分析后还剩1421个lncRNA，下面我们继续做生存分析，还是用easyTCGA,1行代码解决，而且默认是基于最佳截点，同时进行log-rank和单因素Cox两种批量生存分析：
# 加载生存信息，easyTCGA已经一起帮我们下载整理好了
load(file = "G:/easyTCGA_test/output_mRNA_lncRNA_expr/TCGA-COAD_clinical.rdata")
批量生存分析：
surv_res <- batch_survival(exprset = lncrna_expr,
clin = clin_info,
is_count = F, # 不是count数据
optimal_cut = T # 使用最佳截点
)
#save(surv_res, file = "../000files/surv_res.rdata")
查看结果：
load(file = "../000files/surv_res.rdata")
# 我们就用log-rank的结果，或者取logrank和cox的交集也可
log_rank_res <- surv_res$res.logrank
dim(log_rank_res)
## [1] 9578 3
head(log_rank_res)
## gene p.value cutpoint
## 1 MALAT1 0.02133476 5.2170446
## 2 NEAT1 0.26639139 2.7373944
## 3 SNHG25 0.04169925 0.4839386
## 4 AC103702.2 0.05712193 2.6707860
## 5 AL591895.1 0.01227174 6.8008582
## 6 AL354743.2 0.04968819 5.7748580
虽然我们使用了最佳截点，但是并不能保证P值一定是小于0.05的哈，所以我们还是要筛选一下，只要p<0.05的：
suppressMessages(library(dplyr))
log_rank_res <- log_rank_res %>%
filter(p.value < 0.05)
dim(log_rank_res)
## [1] 5157 3
差异分析和生存分析的结果取交集：
deg_logrank <- intersect(log_rank_res$gene, deg_limma$genesymbol)
length(deg_logrank)
## [1] 575
还剩下575个lncRNA，然后就可以提取表达矩阵，做相关性分析了：
lnc_expr <- lncrna_expr_tpm[deg_logrank,]
dim(lnc_expr)
## [1] 575 524
批量相关性分析
下面就是两个表达矩阵的批量相关性分析，这个运算量其实不小，如果你电脑不行，很容易卡死。
在此之前还是要做些小准备，首先是只保留肿瘤样本，然后对表达矩阵进行log2转换，我们在做差异分析和生存分析时并没有做，这是因为easyTCGA会自动帮我们做。
# 只保留肿瘤样本
sample_type <- ifelse(as.numeric(substr(colnames(imm_expr),14,15)) < 10,"tumor","normal")
lnc_expr <- lnc_expr[,sample_type == "tumor"]
imm_expr <- imm_expr[,sample_type == "tumor"]
# 进行log2转换
lnc_expr <- log2(lnc_expr+0.1)
imm_expr <- log2(imm_expr+0.1)
dim(lnc_expr)
## [1] 575 483
dim(imm_expr)
## [1] 1098 483
# 两个表达矩阵的样本顺序完全一样
identical(colnames(lnc_expr),colnames(imm_expr))
## [1] TRUE
1098个mRNA和575个lncRNA批量进行相关性分析，这个计算过程非常大，小电脑会直接卡死！使用linkET::correlate计算是我能找到的最快的方法，比普通的方法快了几十倍！
批量计算两个矩阵的相关性也是一个非常实用的技能哈，大家可以记录一下。
此方法我在之前的推文中也介绍过：免疫浸润常见的可视化方法
# 自定义一个函数
cormatrixes <- function(x,y){
tem <- linkET::correlate(t(x),t(y),engine = "WGCNA")
tem1 <- as.data.frame(tem[[1]])
tem1 <- cbind(rownames(tem1),tem1)
tem1_long <- reshape2::melt(tem1,value.name = "correlation")
tem2 <- as.data.frame(tem[[2]])
tem2 <- cbind(rownames(tem2),tem2)
tem2_long <- reshape2::melt(tem2,value.name = "pvalue")
result <- cbind(tem1_long,tem2_long$pvalue)
names(result) <- c("v1","v2","correlation","pvalue")
return(result)
}
这个函数接受两个表达矩阵，然后返回相关系数和P值：
cor_res <- cormatrixes(lnc_expr,imm_expr)
##
## Using rownames(tem1) as id variables
## Using rownames(tem2) as id variables
tail(cor_res)
## v1 v2 correlation pvalue
## 631345 LINC02868 TNFRSF6B NA NA
## 631346 AP001330.1 TNFRSF6B NA NA
## 631347 AL512358.2 TNFRSF6B NA NA
## 631348 AL391807.1 TNFRSF6B NA NA
## 631349 AC055874.1 TNFRSF6B NA NA
## 631350 AL772337.2 TNFRSF6B NA NA
v1是lncRNA，v2是mRNA，有些结果是NA，因为有的分子可能表达量在所有样本中都一样！我们把NA去掉即可。
最后筛选P值小于0.05和相关系数大于0.7的lncRNA（这个东西没有标准，只要你能解释得通就行！），就可以认为是免疫相关lncRNA了！
cor_res <- na.omit(cor_res)
imm_lnc <- cor_res %>%
filter(correlation > 0.7, pvalue < 0.05) %>%
distinct(v1) %>%
pull(v1)
length(imm_lnc)
## [1] 35
head(imm_lnc)
## [1] "LINC01614" "AC090559.1" "AC016747.1" "LINC02381" "MIR100HG"
## [6] "AC134312.5"
最终得到35个免疫相关的lncRNA，而且是差异表达的，并且对生存有意义的！
如果你嫌这个数字还是太多（或太少）怎么办？在演示过程中至少有3次可以设定阈值的地方，你随便试，直到找到你想要的数字即可！
简单的可视化
使用easyTCGA，我们做个简单的可视化。
# 随便取6个看看
set.seed(1)
markers <- sample(imm_lnc,6,replace = F)
markers
## [1] "LINC02381" "LINC01614" "AC015911.3" "AL365361.1" "AC103591.4"
## [6] "AC092376.2"
tmp <- log2(lncrna_expr_tpm+1)
首先是lncRNA在normal和tumor中的表达量箱线图：
plot_gene(exprset = tmp,marker = markers,group = sample_type,return_data = F)
配对的也可以：
plot_gene_paired(exprset = tmp,marker = markers,return_data = F)
在不同分期中的表达量箱线图，这时候就要把normal样本去掉了，而且由于我们并没有根据stage进行筛选，所以有些P值是大于0.05的：
tmp <- tmp[,sample_type=="tumor"]
# 这个最好再手动整理下
ajcc_satge <- clin_info[sample_type=="tumor","ajcc_pathologic_stage"]
plot_gene(exprset = tmp, marker = markers, group = ajcc_satge, return_data = F)
生存分析也很easy，不过不能像上面这样拼图了：
clin <- clin_info[sample_type=="tumor",c("days_to_last_follow_up","vital_status")]
names(clin) <- c("time","event")
clin$event <- ifelse(clin$event=="Dead",1,0)
# 默认也是基于最佳截点哦
plot_KM(exprset = tmp,marker = markers[1],
clin = clin, return_data = F)
真的是太棒了！
给大家介绍了这么多识别免疫相关lncRNA的方法，大家学会了吗？
联系我们，关注我们
免费QQ交流群1：613637742 免费QQ交流群2：608720452 公众号消息界面关于作者获取联系方式 知乎、CSDN、简书同名账号 哔哩哔哩：阿越就是我
https://mp.weixin.qq.com/s/Wj2dmk5hPhkbGIkdUylGOg