表达矩阵任意两个基因的相关性分析和绘图

花花写于2020-01-31 我家没青菜吃了。

1.输入数据和R包

1rm(list=ls())
2library(ggpubr)
3library(stringr)
4load(file = "for_boxplot.Rdata")

这个Rdata里面有三个数据：

expr和meta是miRNA的表达矩阵和临床信息，由GDC下载整理得到。

mut是突变信息，读取maf得到的数据框筛选了几列得到。

在生信星球公众号回复“box”即可获得。也可参照前面的笔记自己获得。（为什么是box呢，因为和箱线图用的同一个数据，我又懒得改了。）

 1expf = expr[,as.numeric(substr(colnames(expr),14,15)) < 10]
 2expf = expf[,!duplicated(str_sub(colnames(expf),1,12))]
 3x1=str_sub(colnames(expf),1,12)
 4x2=str_to_upper(meta$patient.bcr_patient_barcode)
 5meta = meta[x2 %in% x1 ,]
 6expm = expf[,x1 %in% unique(str_sub(mut$Tumor_Sample_Barcode,1,12))]
 7
 8VHL_mut=substr(as.character(
 9  as.data.frame( mut[mut$Hugo_Symbol=='VHL','Tumor_Sample_Barcode'])[,1] ),
10  1,12)

expf 是表达矩阵expr去掉正常样本后，仅保留肿瘤样本并去重。
expm 是有突变记录的肿瘤样本的表达矩阵。在上一节中有代码过程记录。

2.任意两个基因的相关性分析

2.1 简单绘图

使用ggpbur。

 1dat=data.frame(gene1=log2(expf['hsa-mir-10b',]+1),
 2               gene2=log2(expf['hsa-mir-143',]+1),
 3               stage=as.factor(meta$patient.stage_event.pathologic_stage))
 4
 5
 6sp1 <- ggscatter(dat, x = "gene1", y = "gene2",
 7                add = "reg.line",  # Add regressin line 
 8                add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
 9                conf.int = TRUE # Add confidence interval
10) + stat_cor(method = "pearson", label.x = 15, label.y = 20)
11sp1

2.2 按照stage分组

不仅stage，任意在meta信息中能找到或生产的分组都可以。

1sp2 <- ggscatter( dat, x = "gene1", y = "gene2",
2                color = "stage", palette = "jco",
3                add = "reg.line", conf.int = TRUE) +  stat_cor(aes(color = stage),label.x = 15 )
4sp2

2.3 按照是否突变来分组

理论上某个是否突变并不会改变某两个基因的相关性趋势，如果有这种特殊的突变，打乱了两个基因之间正常的相关关系，机制就有了。可以写循环试一下是否有这样的基因突变。

1dat=data.frame(gene1=log2(expm['hsa-mir-10b',]+1),
2               gene2=log2(expm['hsa-mir-143',]+1),
3               mut= substr(colnames(expm),1,12) %in% VHL_mut)
4sp3 <- ggscatter( dat, x = "gene1", y = "gene2",
5                color = "mut", palette = "jco",
6                add = "reg.line", conf.int = TRUE) +  stat_cor(aes(color = mut),label.x = 15 )
7sp3

3.强迫症闲的没事想给图分身

本以为会很麻烦，可以多浪费一点时间。打开了前年写的ggplot学习笔记，没想到异常简单。
https://www.jianshu.com/u/c93ac360691a（现在看看，写的什么玩意，不想承认是我写的了。）

1sp2  + facet_wrap(~stage, scales = "free_x")

1sp3  + facet_wrap(~mut, scales = "free_x")

4.闲的没事还喜欢玩拼图

自从有了patchwork，我拼图来就雄赳赳气昂昂了。

1library(patchwork)
2sp1 +sp2 +sp3