任意基因的任意分组比较
0.输入数据
rm(list=ls())
load("for_boxplot.Rdata")
这里面的数据：
exp是tumor-normal都有的表达矩阵，exprSet是只有tumor样本的表达矩阵。meta是临床信息表格，Group是tumor-normal分组信息。mut是突变信息，由maf文件读取并取子集得到。
1.比较任意miRNA在tumor和normal样本中的表达量
以hsa-mir-143为例画图，可替换为其他任意miRNA。
table(Group)
#> Group
#> normal tumor
#> 71 522
library(ggstatsplot)
dat = data.frame(gene = exp["hsa-mir-143",],
group = Group)
ggbetweenstats(data = dat, x = group, y = gene,title = "hsa-mir-143")
image.png
2.任意miRNA在任意两个分组中的表达量对比
只要是可以根据临床信息查到或得到的分组，例如生死、人种、阶段，都可以拿来做分组。
需要注意调整样本顺序，一一对应。
#按照生死、人种、分期分组看看
table(meta$patient.vital_status)
#>
#> alive dead
#> 358 158
table(meta$patient.stage_event.pathologic_stage)
#>
#> i ii iii iv
#> 254 55 124 83
table(meta$patient.race)
#>
#> asian black or african american white
#> 8 56 445
dat = data.frame(gene = exprSet["hsa-mir-143",],
vital_status = meta$patient.vital_status,
stage = meta$patient.stage_event.pathologic_stage,
race = meta$patient.race)
p1 = ggbetweenstats(data = dat, x = vital_status, y = gene,title = "hsa-mir-143")
p2 = ggbetweenstats(data = dat, x = stage, y = gene,title = "hsa-mir-143")
p3 = ggbetweenstats(data = dat, x = race, y = gene,title = "hsa-mir-143")
library(patchwork)
p1/p2/p3
image.png
3.根据某个基因是否突变分组比较某miRNA的表达量
dim(exprSet)
#> [1] 552 516
head(mut)
#> Hugo_Symbol Chromosome Start_Position Tumor_Sample_Barcode t_vaf
#> 1: HNRNPCL2 chr1 13115853 TCGA-G6-A8L7-01A-11D-A36X-10 0.2148148
#> 2: ERMAP chr1 42842993 TCGA-G6-A8L7-01A-11D-A36X-10 0.1650165
#> 3: FAAH chr1 46394349 TCGA-G6-A8L7-01A-11D-A36X-10 0.3114754
#> 4: EPS15 chr1 51448116 TCGA-G6-A8L7-01A-11D-A36X-10 0.1677852
#> 5: HMGCS2 chr1 119764248 TCGA-G6-A8L7-01A-11D-A36X-10 0.2539683
#> 6: NOS1AP chr1 162367063 TCGA-G6-A8L7-01A-11D-A36X-10 0.2098765
#> pos
#> 1: chr1:13115853
#> 2: chr1:42842993
#> 3: chr1:46394349
#> 4: chr1:51448116
#> 5: chr1:119764248
#> 6: chr1:162367063
library(stringr)
length(unique(str_sub(mut$Tumor_Sample_Barcode,1,12)))
#> [1] 336
k = str_sub(colnames(exprSet),1,12) %in% unique(str_sub(mut$Tumor_Sample_Barcode,1,12));table(k)
#> k
#> FALSE TRUE
#> 185 331
#516个样本中,有331个有突变信息记录,将这些样本对应的表达矩阵取出来。
expm = exprSet[,k]；
#挑选VHL突变了的病人ID
VHL_mut=str_sub(as.character(
as.data.frame( mut[mut$Hugo_Symbol=='VHL','Tumor_Sample_Barcode'])[,1] ),
1,12)
library(dplyr)
VHL_mut = mut %>%
filter(Hugo_Symbol=='VHL') %>%
as.data.frame() %>%
pull(Tumor_Sample_Barcode) %>%
as.character() %>%
str_sub(1,12)
#false 是未突变样本，true是突变样本
tail(rownames(expm))
#> [1] "hsa-mir-944" "hsa-mir-95" "hsa-mir-96" "hsa-mir-98" "hsa-mir-99a"
#> [6] "hsa-mir-99b"
dat=data.frame(gene=expm['hsa-mir-98',],
mut= str_sub(colnames(expm),1,12) %in% VHL_mut)
ggbetweenstats(data = dat, x = mut, y = gene)
image.png
#可以计算每个基因的p值，找找是不是有显著的。
res.aov <- t.test(gene ~ as.factor(mut), data = dat)
res.aov$p.value
#> [1] 0.9086252
#自己试试能不能写出来这个循环？
