任意基因的任意分組比較
0.輸入數(shù)據(jù)
rm(list=ls())
load("for_boxplot.Rdata")
這里面的數(shù)據(jù)：
exp是tumor-normal都有的表達(dá)矩陣，exprSet是只有tumor樣本的表達(dá)矩陣。meta是臨床信息表格，Group是tumor-normal分組信息。mut是突變信息，由maf文件讀取并取子集得到。
1.比較任意miRNA在tumor和normal樣本中的表達(dá)量
以hsa-mir-143為例畫(huà)圖，可替換為其他任意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在任意兩個(gè)分組中的表達(dá)量對(duì)比
只要是可以根據(jù)臨床信息查到或得到的分組，例如生死、人種、階段，都可以拿來(lái)做分組。
需要注意調(diào)整樣本順序，一一對(duì)應(yīng)。
#按照生死、人種、分期分組看看
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.根據(jù)某個(gè)基因是否突變分組比較某miRNA的表達(dá)量
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個(gè)樣本中,有331個(gè)有突變信息記錄,將這些樣本對(duì)應(yīng)的表達(dá)矩陣取出來(lái)。
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
#可以計(jì)算每個(gè)基因的p值，找找是不是有顯著的。
res.aov <- t.test(gene ~ as.factor(mut), data = dat)
res.aov$p.value
#> [1] 0.9086252
#自己試試能不能寫(xiě)出來(lái)這個(gè)循環(huán)？
