变异数据处理(3)

作者: 白米饭睡不醒 | 来源:发表于2022-02-06 16:01 被阅读0次

    任意基因的任意分组比较

    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
    #自己试试能不能写出来这个循环?
    

    相关文章

      网友评论

        本文标题:变异数据处理(3)

        本文链接:https://www.haomeiwen.com/subject/nzcnhrtx.html