200块的代码我的学徒免费送给你,GSVA和生存分析

(现在学习量和弹幕都非常少,大家的机会来了哦!)

https://www.bilibili.com/video/av81874183

前奏

最近做的生存分析都是奇奇怪怪的,从来没有重复出作者的图。哈哈哈,但是我相信自己的代码是没有错的,只是参数设置跟作者不同而已。

最近跟着Jimmy老师的视频,学习用GSVA做生存分析,反反复复做了几次都无法复现作者的结果。

第一次是我找的基因集跟作者不一致,第二次可能参数不同,也只有一个显著性结果。

看了这篇推文,帮助你节约200块钱,可以多搓一顿火锅了。

看上图👆,不是我骗你,真的有人代码卖200块。好啦,正文开始啦!

学习前必须要知道的生物学知识

01

生存分析(Survival analysis)是指根据试验或调查得到的数据对生物或人的生存时间进行分析和推断,研究生存时间和结局与众多影响因素间关系及其程度大小的方法,也称生存率分析存活率分析。常见分析的用寿命表法、Kaplan-Meier,Cox回归等等. 构建在于病人分组后进行比较!

如果是多个基因,可以通过GSVA来进行分组再去做生存分析。

02

GSVA (The Gene Set Variation Analysis package for microarray and RNA-seq data),算法太复杂,详情搜对应的R包。

我是大标题-1.读文献

第一步需要对文献进行解析,文献来源为Cell. 2018 May  :Chemoresistance Evolution in Triple-Negative Breast Cancer Delineated by Single-Cell Sequencing

暂时跳过了这一步,因为时间有限,就先看了Jimmy老师的推文使用单细胞多组学探索TNBC病人的新辅助化疗疗效

大标题-2.数据库查询指定基因集的基因列表

  • msigdb数据库(找出基因集的基因)

  • 根据一些feature分类;可能是unique,可能是multi

  • 下载msigd数据库的all.gmt数据后,用Linux处理

2-1.下载gmt数据

打开对应的官网(http://software.broadinstitute.org/gsea/msigdb/index.jsp)

👆就是对应的gmt文件,我们下载基因名称的gmt文件

然后打开终端,cd到文件夹

ls -lh

结果如👇

这里做了两次,我分两次讲解,第二次会总结第一次的错误

下面是第一次的workflow

找到你想要的基因集的基因

###寻找基因集
grep -i epithelial_mesenchymal *
#h.all.v6.2.entrez.gmt:HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION    http://www.broadinstitute.org/gsea/msigdb/cards/HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION  1281    1290    12892200    1277    2335    1293    5054    1278    1282    1284    1462    3486    7045    66784060    3915    4015    3918    1490    6876    1294    4017    1292    3685    7058    13074837    7168    1000    4313    1301    7057    2191    633 871 11167   10631   70706696    3371    7980    22795   3693    4314#   7431    4016    10516   1303    2006    10091311    6695    649 9235    3909    7076    5768    7078    7412    3491    10085   800 57692   2192    6443    1893    3908    10272   7169    3624    1601    2014    10409   36784256    7422    2919    7474    6382    5352    5118    26585   3688    50509   388 56542247    6641    1647    4232    131578  4982    966 59  30008   4147    26577   52702817    25890   2517    6586    284217  56937   1296    2201    3485    5217    6385    960 4616    3576    11010   290 64175   7424    4323    6444    5351    4148    10398   6535813 5396    51330   2331    3398    2669    5329    4638    7040    6422    8985    3569333 2199    4487    5806    8325    3725    10979   22943   6591    667 7171    16342697    5376    3487    1314    4035    3673    2316    8076    5744    7049    6424    39565999    1004    6303    4907    1809    5479    7052    6445    3690    8572    115908  18429244    374 3600    4176    2619    5645    23705   5021    7857    6372    4312    7128822 10486   25878   2303    50863   2026    355 627 8038    5817    6387    51599353    4853    79709   2882    7456

这里我下载成entrezID.gmt,后面重新下载。正确基因集的结果如👇

grep -i epithelial_mesenchymal *
#HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION    http://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITIONCOL3A1   COL5A2  COL5A1  FBN1    COL1A1  FN1 COL6A3  SERPINE1    COL1A2  COL4A1  COL4A2  VCAN    IGFBP3  TGFBI   SPARC   LUM LAMC1   LOX LAMC2   CCN2    TAGLN   COL7A1  LOXL2   COL6A2  ITGAV   THBS2   COL16A1 NNMT    TPM1    CDH2    MMP2    COL11A1 THBS1   FAP BGN SERPINH1    FSTL1   POSTN   THY1    SPP1    TNC TFPI2   NID2    ITGB5   MMP3    VIM LOXL1   FBLN5   COL12A1 ELN CDH11   COMP    SPOCK1  BMP1    IL32    LAMA3   TIMP1   QSOX1   TIMP3   VCAM1   CCN1    EDIL3   CALD1   MAGEE1  FBLN1   SGCB    ECM1    LAMA2   FSTL3   TPM2    INHBA   DAB2    EMP3    BASP1   ITGA5   MGP VEGFA   CXCL1   WNT5A   SDC1    PLOD2   PCOLCE  GREM1   ITGB1   COL5A3  RHOB    HTRA1   FGF2    SNTB1   GADD45A MEST    LRRC15  TNFRSF11B   CD59    ACTA2   EFEMP2  MATN2   PCOLCE2 SERPINE2    GPC1    ABI3BP  FUCA1   SLIT3   LAMA1   PMEPA1  COL8A2  FBN2    IGFBP2  PFN2    SDC4    CD44    GADD45B CXCL8   GLIPR1  ANPEP   P3H1    VEGFC   MMP14   SGCD    PLOD1   MATN3   MYL9    SLC6A8  CALU    PRRX1   TNFRSF12A   FMOD    ID2 GEM PLAUR   MYLK    TGFB1   SFRP1   PLOD3   IL6 APLP1   FBLN2   MSX1    PTX3    FZD8    JUN FERMT2  DKK1    SNAI2   DST TPM4    DCN GJA1    PMP22   IGFBP4  COPA    LRP1    ITGA2   FLNA    MFAP5   PTHLH   TGFBR3  SFRP4   LGALS1  RGS4    CDH6    SAT1    NT5E    DPYSL3  PPIB    TGM2    SGCG    ITGB3   PDLIM4  CTHRC1  ECM2    CRLF1   AREG    IL15    MCM7    GAS1    PRSS2   CADM1   OXTR    SCG2    CXCL6   MMP1    TNFAIP3 CAPG    CAP2    MXRA5   FOXC2   NTM ENO2    FAS BDNF    ADAM12  PVR CXCL12  PDGFRB  SLIT2   NOTCH2  COLGALT1    GPX7    WIPF1
grep -i HALLMARK_HYPOXIA *
#msigdb.v7.0.symbols.gmt:HALLMARK_HYPOXIA    http://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_HYPOXIAPGK1   PDK1    GBE1    PFKL    ALDOA   ENO2    PGM1    NDRG1   HK2 ALDOC   GPI MXI1    SLC2A1  P4HA1ADM    P4HA2   ENO1    PFKP    AK4 FAM162A PFKFB3  VEGFA   BNIP3L  TPI1    ERO1A   KDM3A   CCNG2   LDHA    GYS1    GAPDH   BHLHE40 ANGPTL4 JUN SERPINE1    LOX GCK PPFIA4  MAFF    DDIT4   SLC2A3  IGFBPNFIL3  FOS RBPJ    HK1 CITED2  ISG20   GALK1   WSB1    PYGM    STC1    ZNF292  BTG1    PLIN2   CSRP2VLDLR  JMJD6   EXT1    F3  PDK3    ANKZF1  UGP2    ALDOB   STC2    ERRFI1  ENO3    PNRC1   HMOX1   PGF GAPDHS  CHST2   TMEM45A BCAN    ATF3    CAV1    AMPD3   GPC3    NDST1   IRS2    SAP30   GAA SDC4    STBD1IER3   PKLR    IGFBP1  PLAUR   CAVIN3  CCN5    LARGE1  NOCT    S100A4  RRAGD   ZFP36   EGFR    EDN2    IDS CDKN1A  RORA    DUSP1   MIF PPP1R3C DPYSL4  KDELR3  DTNA    ADORA2B HS3ST1  CAVIN1  NR3C1   KLF6    GPC4    CCN1    TNFAIP3 CA12    HEXA    BGN PPP1R15A    PGM2    PIM1    PRDX5   NAGK    CDKN1B  BRS3    TKTL1MT1E   ATP7A   MT2A    SDC3    TIPARP  PKP1    ANXA2   PGAM2   DDIT3   PRKCA   SLC37A4 CXCR4   EFNA3   CP  KLF7    CCN2    CHST3   TPD52   LXN B4GALNT2    PPARGC1A    BCL2    GCNT2   HAS1    KLHL24  SCARBSLC25A1    SDC2    CASP6   VHL FOXO3   PDGFB   B3GALT6 SLC2A5  SRPX    EFNA1   GLRX    ACKR3   PAM TGFBIDCN    SIAH2   PLAC8   FBP1    TPST2   PHKG1   MYH9    CDKN1C  GRHPR   PCK1    INHA    HSPA5   NDST2   NEDD4TPBG   XPNPEP1 IL6 SLC6A6  MAP3K1  LDHC    AKAP12  TES KIF5A   LALBA   COL5A1  GPC1    HDLBP   ILVBLNCAN   TGM2    ETS1    HOXB9   SELENBP1    FOSL2   SULT2B1 TGFB3
grep -i HALLMARK_ANGIOGENESIS *
#msigdb.v7.0.symbols.gmt:HALLMARK_ANGIOGENESIS    http://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_ANGIOGENESIS  VCAN    POSTN   FSTL1   LRPAP1  STC1    LPL VEGFA   PF4 THBD    FGFR1   TNFRSF21    CCND2COL5A2 ITGAV   SERPINA5    KCNJ8   APP JAG1    COL3A1  SPP1    NRP1    OLR1    PDGFA   PTK2    SLCO2A1 PGLYRP1 VAV2    S100A4  MSX1    VTN TIMP1   APOH    PRG2    JAG2    LUM CXCL6
grep -i HALLMARK_PI3K_AKT_MTOR_SIGNALING *
#msigdb.v7.0.symbols.gmt:HALLMARK_PI3K_AKT_MTOR_SIGNALING    http://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_PI3K_AKT_MTOR_SIGNALING   MAPK8   PIK3R3  GRB2    NFKBIB  MAP2K6  MAPK9   AKT1    MAPK1   PLCG1   TRIB3GSK3B  MAP2K3  CDKN1A  RAC1    RIPK1   AKT1S1  ACTR2   PRKAR2A YWHAB   HRAS    PDK1    PIKFYVE TBK1    ACTR3E2F1   MYD88   ITPR2   SQSTM1  RPS6KA1 PTPN11  MAPKAP1 PLCB1   RAF1    CAMK4   RPTOR   CFL1    CDK4    TRAF2GNGT1  UBE2N   ADCY2   CDKN1B  VAV3    FGF6    ECSIT   RALB    ARF1    MKNK1   CDK1    PTEN    ARHGDIA GRK2    FGF17   DDIT3   IRAK4   TIAM1   CDK2    SFN PRKCB   GNA14   EIF4E   CLTC    TSC2    FGF22   PPP1CA  DUSP3HSP90B1    IL4 STAT2   SLA EGFR    PLA2G12A    MAPK10  CALR    THEM4   RIT1    MKNK2   PPP2R1B CAB39ARPC3  PITX2   NCK1    IL2RG   PFN1    FASLG   NOD1    DAPP1   UBE2D3  CAB39   AP2M1   MAP3K7  PRKAG1  CSNK2PRKAA2 ATF1    SLC2A1  PIN1    TNFRSF1A    LCK RPS6KA3 NGF CXCR4   ACACA   SMAD2   PAK4
grep -i BIOCARTA_ECM_PATHWAY *
#msigdb.v7.0.symbols.gmt:BIOCARTA_ECM_PATHWAY    http://www.gsea-msigdb.org/gsea/msigdb/cards/BIOCARTA_ECM_PATHWAY   MAPK3   PIK3CG  HRAS    MYL2    PIK3R1  MAP2K1  RAF1    RHOA    ROCK1   FYN MAPK1   SHC1    TLN1    PFN1    GSN DIAPH1  PIK3CA  ARHGAP5 ITGB1

这里有一个基因集找不到,只有查看来源文献Loss of E-cadherin promotes metastasis via multiple downstream transcriptional pathways

结果是文献套文献,套娃么。。放弃掉,这5个应该也够了。

下面是用哪个metabric下载表达矩阵

2-2.metabric下载表达矩阵

官网为[http://www.cbioportal.org/datasets],这里下载需要的BRCA的数据

下载完成解压后内容如👇

主要用到临床信息、表达矩阵与体细胞变异文件

后面就可以放在R语言了完成🌶

2-3.制作GSVA的输入

主要是制作表达矩阵与临床信息

###读入表达矩阵与临床信息
rm(list = ls())
wkdir=getwd() #为工作目录赋值
options(stringsAsFactors = F)
##临床信息
clin=read.table(file.path(wkdir,'brca_metabric','data_clinical_patient.txt'),
                  header = T,sep='\t')
##表达矩阵
expr=read.table(file.path(wkdir,'brca_metabric','data_expression_median.txt'),
                  header = T,sep='\t')
expr[1:4,1:4]
#  Hugo_Symbol Entrez_Gene_Id  MB.0362  MB.0346
1        RERE            473 8.676978 9.653589
2      RNF165         494470 6.075331 6.687887
3    CD049690             NA 5.453928 5.454185
4    BC033982             NA 4.994525 5.346010
#需要将第一列作为行名并去掉第一、二列
rownames(expr)=expr$Hugo_Symbol
expr=expr[,-c(1,2)]
###这里Jimmy老师为了节省空间,只取了取前两位数
  expr <- apply(expr, 2, function(x){
    as.numeric(format(x,digits = 2)) #digit代表数字位数
  })
  rownames(expr)=gs
  expr=na.omit(expr) #忽略NA值

接下来还需要做一些变化,因为表达矩阵的样本名是以"."分割

而临床信息是"-"分割

clin$PATIENT_ID = gsub('-','.',clin$PATIENT_ID) ##将“-”变为“.”
expr=clin[match(colnames(expr),clin$PATIENT_ID),] #防止错位进行配对矫正
###取TNBC亚型
clin1=clin[clin$CLAUDIN_SUBTYPE%in%c("Basal","claudin-low"),]
clin1=expr[,colnames(expr)%in%clin1$PATIENT_ID]
save(clin1,clin1,file = "GSVA_input.Rdata")

可以看到有427个病人属于这个亚型。如果不看Jimmy老师的视屏,我不会知道哪些是TNBC的。😃,生物学背景真的很重要,见推文没有生物学背景的数据分析很危险

2-4.gmt文件制作

用Linux将前面得到的基因集复制到一个文档,后缀为gmt

cat > GSVA_input.gmt
msigdb.v7.0.symbols.gmt:HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION    http://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITIONCOL3A1   COL5A2  COL5A1  FBN1    COL1A1  FN1 COL6A3  SERPINE1    COL1A2  COL4A1  COL4A2  VCAN    IGFBP3  TGFBI   SPARC   LUM LAMC1   LOX LAMC2   CCN2    TAGLN   COL7A1  LOXL2   COL6A2  ITGAV   THBS2   COL16A1 NNMT    TPM1    CDH2    MMP2    COL11A1 THBS1   FAP BGN SERPINH1    FSTL1   POSTN   THY1    SPP1    TNC TFPI2   NID2    ITGB5   MMP3    VIM LOXL1   FBLN5   COL12A1 ELN CDH11   COMP    SPOCK1  BMP1    IL32    LAMA3   TIMP1   QSOX1   TIMP3   VCAM1   CCN1    EDIL3   CALD1   MAGEE1  FBLN1   SGCB    ECM1    LAMA2   FSTL3   TPM2    INHBA   DAB2    EMP3    BASP1   ITGA5   MGP VEGFA   CXCL1   WNT5A   SDC1    PLOD2   PCOLCE  GREM1   ITGB1   COL5A3  RHOB    HTRA1   FGF2    SNTB1   GADD45A MEST    LRRC15  TNFRSF11B   CD59    ACTA2   EFEMP2  MATN2   PCOLCE2 SERPINE2    GPC1    ABI3BP  FUCA1   SLIT3   LAMA1   PMEPA1  COL8A2  FBN2    IGFBP2  PFN2    SDC4    CD44    GADD45B CXCL8   GLIPR1  ANPEP   P3H1    VEGFC   MMP14   SGCD    PLOD1   MATN3   MYL9    SLC6A8  CALU    PRRX1   TNFRSF12A   FMOD    ID2 GEM PLAUR   MYLK    TGFB1   SFRP1   PLOD3   IL6 APLP1   FBLN2   MSX1    PTX3    FZD8    JUN FERMT2  DKK1    SNAI2   DST TPM4    DCN GJA1    PMP22   IGFBP4  COPA    LRP1    ITGA2   FLNA    MFAP5   PTHLH   TGFBR3  SFRP4   LGALS1  RGS4    CDH6    SAT1    NT5E    DPYSL3  PPIB    TGM2    SGCG    ITGB3   PDLIM4  CTHRC1  ECM2    CRLF1   AREG    IL15    MCM7    GAS1    PRSS2   CADM1   OXTR    SCG2    CXCL6   MMP1    TNFAIP3 CAPG    CAP2    MXRA5   FOXC2   NTM ENO2    FAS BDNF    ADAM12  PVR CXCL12  PDGFRB  SLIT2   NOTCH2  COLGALT1    GPX7    WIPF1
msigdb.v7.0.symbols.gmt:HALLMARK_HYPOXIA    http://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_HYPOXIAPGK1   PDK1    GBE1    PFKL    ALDOA   ENO2    PGM1    NDRG1   HK2 ALDOC   GPI MXI1    SLC2A1  P4HA1ADM    P4HA2   ENO1    PFKP    AK4 FAM162A PFKFB3  VEGFA   BNIP3L  TPI1    ERO1A   KDM3A   CCNG2   LDHA    GYS1    GAPDH   BHLHE40 ANGPTL4 JUN SERPINE1    LOX GCK PPFIA4  MAFF    DDIT4   SLC2A3  IGFBPNFIL3  FOS RBPJ    HK1 CITED2  ISG20   GALK1   WSB1    PYGM    STC1    ZNF292  BTG1    PLIN2   CSRP2VLDLR  JMJD6   EXT1    F3  PDK3    ANKZF1  UGP2    ALDOB   STC2    ERRFI1  ENO3    PNRC1   HMOX1   PGF GAPDHS  CHST2   TMEM45A BCAN    ATF3    CAV1    AMPD3   GPC3    NDST1   IRS2    SAP30   GAA SDC4    STBD1IER3   PKLR    IGFBP1  PLAUR   CAVIN3  CCN5    LARGE1  NOCT    S100A4  RRAGD   ZFP36   EGFR    EDN2    IDS CDKN1A  RORA    DUSP1   MIF PPP1R3C DPYSL4  KDELR3  DTNA    ADORA2B HS3ST1  CAVIN1  NR3C1   KLF6    GPC4    CCN1    TNFAIP3 CA12    HEXA    BGN PPP1R15A    PGM2    PIM1    PRDX5   NAGK    CDKN1B  BRS3    TKTL1MT1E   ATP7A   MT2A    SDC3    TIPARP  PKP1    ANXA2   PGAM2   DDIT3   PRKCA   SLC37A4 CXCR4   EFNA3   CP  KLF7    CCN2    CHST3   TPD52   LXN B4GALNT2    PPARGC1A    BCL2    GCNT2   HAS1    KLHL24  SCARBSLC25A1    SDC2    CASP6   VHL FOXO3   PDGFB   B3GALT6 SLC2A5  SRPX    EFNA1   GLRX    ACKR3   PAM TGFBIDCN    SIAH2   PLAC8   FBP1    TPST2   PHKG1   MYH9    CDKN1C  GRHPR   PCK1    INHA    HSPA5   NDST2   NEDD4TPBG   XPNPEP1 IL6 SLC6A6  MAP3K1  LDHC    AKAP12  TES KIF5A   LALBA   COL5A1  GPC1    HDLBP   ILVBLNCAN   TGM2    ETS1    HOXB9   SELENBP1    FOSL2   SULT2B1 TGFB3
msigdb.v7.0.symbols.gmt:HALLMARK_ANGIOGENESIS    http://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_ANGIOGENESIS  VCAN    POSTN   FSTL1   LRPAP1  STC1    LPL VEGFA   PF4 THBD    FGFR1   TNFRSF21    CCND2COL5A2 ITGAV   SERPINA5    KCNJ8   APP JAG1    COL3A1  SPP1    NRP1    OLR1    PDGFA   PTK2    SLCO2A1 PGLYRP1 VAV2    S100A4  MSX1    VTN TIMP1   APOH    PRG2    JAG2    LUM CXCL6
msigdb.v7.0.symbols.gmt:HALLMARK_PI3K_AKT_MTOR_SIGNALING    http://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_PI3K_AKT_MTOR_SIGNALING   MAPK8   PIK3R3  GRB2    NFKBIB  MAP2K6  MAPK9   AKT1    MAPK1   PLCG1   TRIB3GSK3B  MAP2K3  CDKN1A  RAC1    RIPK1   AKT1S1  ACTR2   PRKAR2A YWHAB   HRAS    PDK1    PIKFYVE TBK1    ACTR3E2F1   MYD88   ITPR2   SQSTM1  RPS6KA1 PTPN11  MAPKAP1 PLCB1   RAF1    CAMK4   RPTOR   CFL1    CDK4    TRAF2GNGT1  UBE2N   ADCY2   CDKN1B  VAV3    FGF6    ECSIT   RALB    ARF1    MKNK1   CDK1    PTEN    ARHGDIA GRK2    FGF17   DDIT3   IRAK4   TIAM1   CDK2    SFN PRKCB   GNA14   EIF4E   CLTC    TSC2    FGF22   PPP1CA  DUSP3HSP90B1    IL4 STAT2   SLA EGFR    PLA2G12A    MAPK10  CALR    THEM4   RIT1    MKNK2   PPP2R1B CAB39ARPC3  PITX2   NCK1    IL2RG   PFN1    FASLG   NOD1    DAPP1   UBE2D3  CAB39   AP2M1   MAP3K7  PRKAG1  CSNK2PRKAA2 ATF1    SLC2A1  PIN1    TNFRSF1A    LCK RPS6KA3 NGF CXCR4   ACACA   SMAD2   PAK4
msigdb.v7.0.symbols.gmt:BIOCARTA_ECM_PATHWAY    http://www.gsea-msigdb.org/gsea/msigdb/cards/BIOCARTA_ECM_PATHWAY   MAPK3   PIK3CG  HRAS    MYL2    PIK3R1  MAP2K1  RAF1    RHOA    ROCK1   FYN MAPK1   SHC1    TLN1    PFN1    GSN DIAPH1  PIK3CA  ARHGAP5 ITGB1

结果如👇

大标题-3.GSVA

有了上面的文件就可以利用GSVA计算了

library(GSVA)
library(stringr)
gs=readLines("GSVA_input.gmt")#读入基因集
gs1 <- lapply(gs, function(x)  
  {x=strsplit(x,"\t")[[1]]   #通过\t分割
   y=x[3:length(x)]          #不取前两行,因为不是基因
   return(y)})               #返回y值
X=expr
es.max <- gsva(X, gs1,mx.diff=FALSE, verbose=FALSE, parallel.sz=1)
rownames(es.max)=unlist(lapply(gs, function(x) 
{x=strsplit(x,"\t")[[1]][1]
x=strsplit(x,":")[[1]][2]}))
es.dif <- gsva(X, gs1, mx.diff=TRUE, verbose=FALSE, parallel.sz=1)
rownames(es.dif)=rownames(es.max)       #行名命名
save(es.max,es.dif,gs1,file = "es_output.Rdata")
save(es.max,es.dif,gs1,file = "es_output.Rdata")

👆是两种算法的差异结果,具体我也没有看太懂,根据文章设置参数就好。

大标题-4.生存分析

两部分:es.max/es.dif

es.max

library(survival)
library(survminer)
phe$event=ifelse(phe$OS_STATUS=='DECEASED',1,0) #分为死亡与生存2种情况
phe$time=as.numeric(phe$OS_MONTHS) #根据月份划分
rownames(phe)=gsub("-",".",phe$PATIENT_ID) #同上面一样需要变一下分割符
ids=intersect(rownames(phe),colnames(es.max)) #取交集
phe=phe[ids,] #取子集
es.max=es.max[,ids] #取子集

筛选完成后只剩398个样本

根据文献大于0.1为“high,小于”-0.1“为low,方法就是用的GSVA,所以这里我们对其进行分组

for (i in 1:5){
phe1=phe[abs(es.max[i,]) > 0.1,] 
es.max1=es.max[,abs(es.max[i,]) > 0.1]
phe1$group_list=ifelse(es.max1[i,] > 0.1,'high','low')
# 利用ggsurvplot快速绘制漂亮的生存曲线图
sfit <- survfit(Surv(time, event)~group_list, data=phe1)
p=ggsurvplot(sfit, conf.int=F, pval=TRUE,legend.title=rownames(es.max)[i])
print(p$plot)
ggsave(filename = paste0("0.1",rownames(es.max)[i],".pdf"))}
dev.off()

结果如👇,没有一个显著

但当我把参数调为0.2的时候居然第一个基因集ECM显著,但是并不是文章中显著的2个基因集(AKT、HYPOXIA)

下面我用另一种方法试一下

es.dif

方法跟👆差不多

phe=clin1[clin1$OS_STATUS %in% c('DECEASED','LIVING'),]
phe$event=ifelse(phe$OS_STATUS=='DECEASED',1,0)
phe$time=as.numeric(phe$OS_MONTHS)
colnames(phe)
rownames(phe)=gsub("-",".",phe$PATIENT_ID)
ids1=intersect(rownames(phe),colnames(es.dif))
phe=phe[ids1,]
es.dif=es.dif[,ids1]
for (i in 1:5){
  phe2=phe[abs(es.dif[i,]) > 0.1,] 
  es.dif1=es.dif[,abs(es.dif[i,]) > 0.1]
  phe2$group_list=ifelse(es.dif1[i,] > 0.1,'high','low')
  # 利用ggsurvplot快速绘制漂亮的生存曲线图
  sfit <- survfit(Surv(time, event)~group_list, data=phe2)
  p=ggsurvplot(sfit, conf.int=F, pval=TRUE,legend.title=rownames(es.max)[i])
  print(p$plot)
  ggsave(filename = paste0("0.1_i",rownames(es.dif)[i],".pdf"))} 

结果如👇

参数为0.1,只有第一个基因集ECM显著

参数为0.2时,也是只有第一个基因集ECM显著


大标题-插曲

前面在下载gmt数据那里做了两次,第一次在上面已经讲解了,下面是第二次的探索过程,第二次会总结第一次的错误。

那么为什么有如此大的差别呢?方法与文中一致,那么只可能是数据出现的问题。检查文章,果然发现筛选的基因集不一致。前面我们主要筛选的是HALLMARK的基因集,是由多个已知的基因集构成的超基因集

而文章里用的是以前文章发表的基因集

那应该怎么搜索基因集呢?我找到了规律,只需要搜索作者名字+基因集名称就可以了

grep HARRIS_HYPOXIA *

结果如下面

这就是我们所需要的基因集了,其他基因集方法类似

grep HARRIS_HYPOXIA *
#msigdb.v7.0.symbols.gmt:HARRIS_HYPOXIA    http://www.gsea-msigdb.org/gsea/msigdb/cards/HARRIS_HYPOXIA PLAUR   PGK1    JUN VEGFA   XRCC5   RP1 HIF1A   CA9 PGF BNIP3   PKM ALDOA   BNIP3L  BHLHE40 HGF TFF3    ENO1    HK2 F3  CP  SPP1    FOS HMOX1   DDIT3   HDAC9   AK3 EDN1    P4HA1   TH  EPAS1   IGFBP3  CCL2    STC1    LDHA    COL5A1  FLT1    CA12    PDGFB   TEK IGF2    PTGS2   CD99    TXN CDKN1A  PFKL    XRCC6   SAT1    TGFB1   FTL TFRC    LRP8    CXCL8   NFKB1   TF  TAGLN   CCNG2   TGM2    EDN2    IL6 BIK PFKP    ANGPT2  TGFA    EPO SLC2A1  ADM VIM SLC2A3  IGFBP1  L1CAM   GAPDH   CDKN1B  MIF TGFB3   IGFBP2  APEX1   FGF3    PRPS1   HK1 MMP13   ENPEP
grep CREIGHTON_AKT1_SIGNALING *
#msigdb.v7.0.symbols.gmt:CREIGHTON_AKT1_SIGNALING_VIA_MTOR_UP    http://www.gsea-msigdb.org/gsea/msigdb/cards/CREIGHTON_AKT1_SIGNALING_VIA_MTOR_UP   NEDD8   BRMS1   CDK16   SPINT1  DDR1    LASP1   CLSTN1  MVK PRKCD   UBE2M   AC004156.1  NEU1    CORO1B  CYB561  KCTD5   SYNJ2BP CDC34   TJP3    HNRNPAB SLC37A1 NECTIN2 ADIPOR1 BSG BIK AKT1    DUSP10  MMP15   POR GPX4    ARHGEF16    CLDN3   TMED10  TOLLIP  PMPCA
#msigdb.v7.0.symbols.gmt:CREIGHTON_AKT1_SIGNALING_VIA_MTOR_DN    http://www.gsea-msigdb.org/gsea/msigdb/cards/CREIGHTON_AKT1_SIGNALING_VIA_MTOR_DN   TNFRSF12A   RGL2    PFKL    PPP4C   CIB1    ATP6V0B PPP2R1A KRT8    DHCR7   MRPS7   TUBB4B  YWHAB   ALDOA   PAFAH1B3    GOT1    TOM1    ATP6AP1 CTSA    ATP6V0C MIF GPI ATP6V1F TSPAN1
grep ONDER_CDH1 *
#msigdb.v7.0.symbols.gmt:ONDER_CDH1_TARGETS_3_UP    http://www.gsea-msigdb.org/gsea/msigdb/cards/ONDER_CDH1_TARGETS_3_UP    TSC22D3 KDELR3  ZBTB16  MAGED1  STC2    CYP1B1  MAGED2  LEPR    ALDH6A1 TCEA2   WIPI1   KCNMA1  RNASE4  S100A8  CRIP2   FTL PRKCA   KLF9
#msigdb.v7.0.symbols.gmt:ONDER_CDH1_TARGETS_3_DN    http://www.gsea-msigdb.org/gsea/msigdb/cards/ONDER_CDH1_TARGETS_3_DN    CXCL1   EPN3    PTGS2   DSC2    NAV3    EDN1    SOX9    KRT13   CWH43   LCN2    TGFA    SCEL    HS3ST2  C1orf116    FST HBEGF   CEACAM6 SPRR2D  IL1R2   MFAP5   TMPRSS11E   CRCT1   IVL KLK10   IL36G   EHF ADGRE2  SERPINB2    P2RY2   SERPINB13   ATP12A  IL1B    CXCL8   KRT8    MYO5C   KLK11   KLK7    SPAG1   ARL4C   MARC1   CYB5R2  CXCL2   ST6GALNAC5  ANO1    KRT15   PI3 GPRC5A  SERPINB1    TP63    PTHLH   CSF3    DEPP1   DUSP6   ARHGAP25    CXCL3   CSF2    VGLL1   MMP10   S100A7

其他的3个基因集文章没有并给文献,只能默认用HALLMARK的基因集

下面是第二次的workflow

其实跟第一次类似,不同的地方做注释,相似的地方不在赘述。这里以hypoxia基因集作为🌰,其他的基因集类似

rm(list = ls())
options(stringsAsFactors = F)
load("data/metabric_clinical.Rdata")
load("data/metabric_expression.Rdata")
clin1=clin[clin$CLAUDIN_SUBTYPE%in%c("Basal","claudin-low"),]
class(expr)
class(expr1)
library(GSVA)
library(stringr)
gs=readLines("harris_hypoxia.gmt") #记得这里要替换成不同基因集的gmt
gs1 <- lapply(gs, function(x) 
{x=strsplit(x,"\t")[[1]]
y=x[3:length(x)]
return(y)})
X=expr
hy.max <- gsva(X, gs1,mx.diff=FALSE, verbose=FALSE, parallel.sz=1)#max法
rownames(hy.max)=unlist(lapply(gs, function(x) 
{x=strsplit(x,"\t")[[1]][1]
x=strsplit(x,":")[[1]][2]}))
hy.dif <- gsva(X, gs1, mx.diff=TRUE, verbose=FALSE, parallel.sz=1)#diff法
rownames(hy.dif)=rownames(hy.max)
save(hy.max,hy.dif,gs1,file = "hy_hypoxia.Rdata")#存储一下
load("hy_hypoxia.Rdata")
library(survival)
library(survminer)
library(stringr)
dim(clin)
clin[1:4,1:4]
# 终点事件(outcome event)又称失效事件(failure event) 或死亡事件(death event)  
# 这种分组资料的生存分析常采用寿命表法(life-table method)
# 生存分析也经常采用Kaplan-Meier曲线及log-rank检验
table(clin$VITAL_STATUS)
table(clin$OS_STATUS)
# 理论上这个时候需要区分 Died of Disease和Other Caushy
phe=clin1[clin1$OS_STATUS %in% c('DECEASED','LIVING'),]
phe$event=ifelse(phe$OS_STATUS=='DECEASED',1,0)
phe$time=as.numeric(phe$OS_MONTHS)
colnames(phe)
rownames(phe)=gsub("-",".",phe$PATIENT_ID)
ids=intersect(rownames(phe),colnames(hy.max))
phe=phe[ids,]
hy.max=as.matrix(hy.max)#因为只有一列,需要转换成矩阵
hy.max=hy.max[,ids]
hy.max=as.matrix(hy.max)#因为只有一列,需要转换成矩阵
phe[1:4,1:4]
phe1=phe[abs(hy.max[1,]) > 0.1,] 
hy.max1=hy.max[,abs(hy.max[1,]) > 0.1]
hy.max1=t(as.matrix(hy.max1)) #因为只有一列,需要转换成矩阵并转置,不然会报错
phe1$group_list=ifelse(hy.max1[1,] > 0.1,'high','low')
# 利用ggsurvplot快速绘制漂亮的生存曲线图
sfit1 <- survfit(Surv(time, event)~group_list, data=phe1)
p1=ggsurvplot(sfit1, conf.int=F, pval=TRUE,legend.title="hypoxia")
p1
dev.off()

###diff
phe=clin1[clin1$OS_STATUS %in% c('DECEASED','LIVING'),]
phe$event=ifelse(phe$OS_STATUS=='DECEASED',1,0)
phe$time=as.numeric(phe$OS_MONTHS)
colnames(phe)
rownames(phe)=gsub("-",".",phe$PATIENT_ID)
ids1=intersect(rownames(phe),colnames(hy.dif))
phe=phe[ids1,]
hy.dif=hy.dif[,ids1]
hy.dif=as.matrix(hy.dif)
phe2=phe[abs(hy.dif[1,]) > 0.1,] 
hy.dif1=hy.dif[,abs(hy.dif[1,]) > 0.1]
hy.dif1=t(as.matrix(hy.dif1))
phe2$group_list=ifelse(hy.dif1[1,] > 0,'high','low')
# 利用ggsurvplot快速绘制漂亮的生存曲线图
sfit2 <- survfit(Surv(time, event)~group_list, data=phe2)
p2=ggsurvplot(sfit2, conf.int=F, pval=TRUE,legend.title="hypoxia")
p2
dev.off()

最终结果如👇

0.1的阈值,MAX法结果只有文章中显著的1个基因集(AKT)显著,另一个HYPOXIA并不显著

0.1的阈值,diff法结果文章中显著的1个基因集(AKT)显著,另一个HYPOXIA并不显著,另外ECM基因集也显著

总结

综上所述,生存分析就是任人打扮的小姑娘,不同参数,不同数据集,不同基因集,不同方法结果就会不一样,只要细心调整参数,总有你想要的结果,虽然听起来有点不寒而栗!

(0)

相关推荐