Moonlight：一种识别生物标志物在不同肿瘤类型和分期中作为癌基因或肿瘤抑制因子的多种作用的方法 by 生信菜鸟团

Moonlight: an approach to identify multiple role of biomarkers as oncogene or tumorsuppressor in different tumor types and stages [1,2].

2020第一次发表（MoonlightR）Colaprico A, Olsen C, Bailey MH, et al. Interpreting pathways to discover cancer driver genes with Moonlight. Nat Commun. 2020;11(1):69. Published 2020 Jan 3. doi:10.1038/s41467-019-13803-0(2022-11-20 更新，包更新名为 Moonlight2R)
Astrid Saksager, Mona Nourbakhsh, Nikola Tom,et.al.An Automatized Workflow to Study Mechanistic Indicators for Driver Gene Prediction with Moonlight. bioRxiv 2022.11.18.517066; doi: https://doi.org/10.1101/2022.11.18.517066

1. 引入

Moonlight 包是2020发表在Nature communication(2022 更新版目前发表在BioRxiv，包名为Moonlight2R)。简单从包的名字上看，会让人摸不着头脑——好端端的生物分析的R包和 “”月光“有什么关系。实际上而言，就“Moonlight”单词检索Pubmed，会发现这个单词广泛应用于“生物学现象”的描述。稍稍检索发现：

Cite from 微信公众号“英语” 《老外常说的 "moonlight" 是什么意思？可不止是“月光”哦》

根据作者对Moonlight 一词的阐述：

The name refers to (i) the concept of protein moonlighting (or gene sharing) is a phenomenon by which a protein can perform more than one function, and (ii) casting genes in a new light can lead to improved treatment regimens and prognostic indicators.

就纯粹的望文生义的说，如果说月光是太阳光的反射，那么Moonlight 效应就可以类似于理解为旁路效应或者是说非主要效应。引申义，“兼职”的翻译或许也还能说的过去。

言归正传，标题已经大部分阐释了Moonlight R包的应用场景：回答生物分子是否为促癌基因或者抑癌基因的一个方法。那么，这个包解决问题的逻辑是怎样的呢？具体如下：

获取数据：函数getDataTCGA和getDataGEO，是作者封装的能够从TCGA中及筛选匹配的GEO数据集中下载数据，进行后续分析；
差异分析(DPA)：DPA函数是作者封装的能够进行差异分析的函数，TCGA数据依托于TCGAbiolink R 包的差异分析函数，GEO数据依托于edgR函数，分组包括正常和肿瘤、或正常和 I 期、正常和分子亚型等。
功能富集分析(FEA)：FEA 函数是作者封装用于富集分析的函数，作者从文献以及,Msigdb数据库整理获得一个“Diseaselist” 基因集，使用作者构建的方法进行富集分析
基因调控网络构建(GRN)：GRN 函数是通过Parmigene 包的K邻近聚类分析获取基因调控网络的工具，最终会obtaining a set of regulated genes for each DEG。
上游调节分析（URA）：这里可以说是这个R包的核心之一，其应用 z-score，计算参与特定功能的所有基因的所有预测效应的总和与所有基因数量的平方根之间的比率。
模式识别分析（PRA）：PRA函数完成的是根据预定义的生物学过程或者机器学习的方式识别候选 TSG（tumor suppressor genes,增值Down，凋亡UP）和 OCG（oncogenes , 增值Up，凋亡Down）。
驱动基因突变分析（DMA）:更新后Moonlight2R新增的功能，DMA函数是对PRA输出（OCG）进一步处理，主要内容是通过Cscape-somatic算法[3]分析MAF(Mutation Annotation Format)文件，对OCG基因是否为突变驱动基因或者为“passenger” 基因进行进一步评估。

Rogers MF, Gaunt TR, Campbell C. CScape-somatic: distinguishing driver and passenger point mutations in the cancer genome [published correction appears in Bioinformatics. 2021 Oct 25;:]. Bioinformatics. 2020;36(12):3637-3644. doi:10.1093/bioinformatics/btaa242

以上是整个分析的主要函数及流程，而R包中的moonlight 函数则是对上述函数的进一步封装。下图是其简要的流程图示，以及部分绘图函数的结果：

图1. Moonlight pipeline 图示

2. 安装及加载

MoonlightR 可以通过Bioconductor的方式进行安装，

if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install("MoonlightR")
library(edgeR)
library(EDASeq)
library(MoonlightR)

如果安装失败，可以尝试更换Rstudio默认的镜像，自测清华的镜像安装不太行。Biocondutor 的R包一般安装不会有什么太大的问题。这个R包需要”edgeR"包和"EDAseq"包，所以在运行前需要预先加载。
更新后的Moonlight2R函数提供的是github安装方式，如下

devtools::install_github(repo = "ELELAB/Moonlight2R")

如果想要看demo 文档，需要安装BiocStyle包

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("BiocStyle")
devtools::install_github(repo = "ELELAB/Moonlight2R", build_vignettes = TRUE)
vignette( "Moonlight2R", package="Moonlight2R")

3. 主要函数解析

3.1 数据下载（getDataTCGA and getDataGEO）

Moonlight 函数内嵌了 TCGAbiolink 包和GEOquery R包，所以在MoonlightR的应用环境中，可以直接使用函数获取TCGA的18种癌症数据库以及可以通过GEOquery 获取的GEO数据。在作者给定的脚本文件中给出了使用方法，具体如下：

3.1.1 TCGA数据下载

方式一：按癌症类型和数据类型搜索 [Gene expression]

dataFilt <- getDataTCGA(cancerType = "LUAD", 
                          dataType = "Gene expression",
                          directory = "data",
                          nSample = 4)
 class(dataFilt)
# [1] "matrix" "array" 
head(dataFilt)
##          TCGA-44-6146-01B-04R-A277-07 TCGA-55-A48X-01A-11R-A24H-07 
## 100134869                          251                           19                  
## 155060                             851                          512                   
## 8225                               216                         1111                      
## A1BG                               109                          116   
## A2LD1                              184                          133

方式二：按癌症类型和数据类型搜索 [Methylation]

dataFilt <- getDataTCGA(cancerType = "BRCA", 
                          dataType = "Methylation", 
                          directory = "data",nSample = 4)

深度解析getDataTCGA函数的代码，可以发现，其中主要内嵌的就是TCGAbiolink的下载流程和相关参数，而在这里作者大概是重点关注基因表达及甲基化相关的内容，所以只覆盖这两类数据的简单快速下载

getDataTCGA()
function (cancerType, dataType, directory, cor.cut = 0.6, qnt.cut = 0.25, 
          nSample, stage = "ALL", subtype = 0, samples = NULL) 
{
  DiseaseList <- get("DiseaseList")
  GDCprojects <- get("GDCprojects")
  geneInfo <- get("geneInfo")
  CancerProject <- paste0("TCGA-", cancerType)
  DataDirectory <- paste0(directory, "GDC_", gsub("-", "_", 
                                                  CancerProject))
  if (dataType == "Gene expression") {
    ……
                                     }
  else if (dataType == "Methylation") {
   ……
  }
}
<bytecode: 0x55b600b5cbc8>
  <environment: namespace:MoonlightR>

方式二：按癌症类型和数据类型搜索 [Mutation]
Moonlight2R 说明文档中新增MAF文件下载方式，简单来说就是依托TCGAbiolink R包下载TCGA突变数据。

library(TCGAbiolinks)
query <- GDCquery(
  project = "TCGA-LUAD", 
  data.category = "Simple Nucleotide Variation",
  data.type = "Masked Somatic Mutation",
  access = "open", 
  legacy = F,
  workflow.type = "Aliquot Ensemble Somatic Variant Merging and Masking")

GDCdownload(query, directory = "data")
dataMAF <- GDCprepare(query, 
                        directory = "data")

dataMAF <- dataMAF %>% sample_n(size = 3000, replace = FALSE)

TCGAbiolink R包提供的MAF的文件下载说明文档链接：https://bioconductor.org/packages/release/bioc/vignettes/TCGAbiolinks/inst/doc/mutation.html

3.1.2 GEO数据下载

GEO 函数命令如下，根据作者的描述，这里的GSE数据集仅限于作者GEO_TCGAtab表格中提供的GSE 数据集

dataFilt <- getDataGEO(GEOobject = "GSE20347",platform = "GPL571")
# Found 1 file(s)
# GSE20347_series_matrix.txt.gz
###  或者  ####
dataFilt <- getDataGEO(TCGAtumor = "ESCA")

getDataGEO 实际上就是个封装好GEOquery 在内的一个获取数据集的函数，而对于现在的大多数的GSE数据集来说，GEOquery 的可用性一般。

getDataGEO
function (GEOobject = "GSE39004", platform = "GPL6244", TCGAtumor = NULL) 
{
    GEO_TCGAtab <- get("GEO_TCGAtab")
    if (length(TCGAtumor) != 0) {
        GEOobject <- GEO_TCGAtab[GEO_TCGAtab$Cancer == TCGAtumor, 
            "Dataset"]
        platform <- GEO_TCGAtab[GEO_TCGAtab$Cancer == TCGAtumor, 
            "Platform"]
    }
    gset <- getGEO(GEOobject, GSEMatrix = TRUE, AnnotGPL = TRUE)
    if (length(gset) > 1) 
        idx <- grep(platform, attr(gset, "names"))
    else idx <- 1
    gset <- gset[[idx]]
    fvarLabels(gset) <- make.names(fvarLabels(gset))
    return(gset)
}
<bytecode: 0x557cab26d530>
<environment: namespace:MoonlightR>

GEO_TCGAtab 是作者整理好的18个TCGA数据集及其相对应的GSE序列号，以及一些差异基因信息

knitr::kable(GEO_TCGAtab, digits = 2, 
+             caption = "Table with GEO data set matched to one 
+             of the 18 given TCGA cancer types ",
+             row.names = TRUE)
# Table: Table with GEO data set matched to one 
#              of the 18 given TCGA cancer types 
#
# |   |Cancer |TP   |NT  |DEG. |Dataset  |TP.1 |NT.1 |Platform |DEG.. |Common |GEO_Normal      |GEO_Tumor               |
# |:--|:------|:----|:---|:----|:--------|:----|:----|:--------|:-----|:------|:---------------|:-----------------------|
# |1  |BLCA   |408  |19  |2937 |GSE13507 |165  |10   |GPL65000 |2099  |896    |control         |cancer                  |
# |2  |BRCA   |1097 |114 |3390 |GSE39004 |61   |47   |GPL6244  |2449  |1248   |normal          |Tumor                   |
# |3  |CHOL   |36   |9   |5015 |GSE26566 |104  |59   |GPL6104  |3983  |2587   |Surrounding     |Tumor

3.2 DPA 差异基因分析

##对于TCGA数据， 应用TCGAanalyze_DEA（基因表达数据） 或TCGAanalyze_DMC函数（甲基化数据）
##对于 GEO数据，则应用limma 的差异基因分析流程
dataFilt <- getDataTCGA(cancerType = "LUAD", 
                          dataType = "Gene expression",
                          directory = "data",
                          nSample = 4)
dataDEGs <- DPA(dataFilt = dataFilt,
                dataType = "Gene expression")  
#如 为“Methylation” 则相应修改；DPA函数内根据dataType参数设置了条件语句
# 函数来源于TCGAbiolink
head(dataDEGs)
#            logFC   logCPM       LR       PValue          FDR
# ABCA12  8.475687 2.588540 80.86120 2.421419e-19 1.803715e-15
# ABCA8  -2.756115 5.109457 14.23398 1.614294e-04 2.530602e-03
# ABCC2   6.763848 4.276506 17.73730 2.536046e-05 5.978166e-04
# ABCC3   3.258925 7.397783 19.43158 1.042686e-05 2.987296e-04
# ABP1    7.000924 6.369022 30.85150 2.785464e-08 2.901947e-06
# ACADL  -2.058819 4.093631 11.41242 7.295472e-04 8.178175e-03

DPA
function (dataType, dataFilt, dataConsortium = "TCGA", fdr.cut = 0.01, 
          logFC.cut = 1, diffmean.cut = 0.25, samplesType, colDescription, 
          gset, gsetFile = "gsetFile.RData") 
{
  if (dataConsortium == "GEO") {
    ……
    fit <- lmFit(gset, design)
    xContrast <- c("G1-G0")
    cont.matrix <- makeContrasts(xContrast, levels = design)
    fit2 <- contrasts.fit(fit, cont.matrix)
    fit2 <- eBayes(fit2, 0.01)
    tT <- topTable(fit2, adjust.method = "fdr", sort.by = "B", 
    ……
  }
  else if (dataConsortium == "TCGA") {
    if (dataType == "Gene expression") {
    ……
      dataDEGs <- TCGAanalyze_DEA(mat1 = dataFilt[, dataSmNT], 
                                  mat2 = dataFilt[, dataSmTP], Cond1type = "Normal", 
                                  Cond2type = "Tumor", fdr.cut = fdr.cut, logFC.cut = logFC.cut, 
                                  method = "glmLRT")
    ……
    }
    else if (dataType == "Methylation") {
     ……
      system.time(cancer.met <- TCGAanalyze_DMC(dataFilt, 
                                                groupCol = "shortLetterCode", group1 = "TP", 
                                                group2 = "NT", p.cut = fdr.cut, diffmean.cut = 0.25, 
                                                legend = "State", plot.filename = "test.png"))
   ……
    }
  }
}
<bytecode: 0x55b615cbdde8>
  <environment: namespace:MoonlightR>
差异分析的结果可以用TCGAbiolink 函数的绘制火山图
library(TCGAbiolinks)
TCGAVisualize_volcano(dataDEGs$logFC, dataDEGs$FDR,
                      filename = "DEGs_volcano.png",
                      x.cut = 1,
                      y.cut = 0.05,
                      names = rownames(DEGsmatrix),
                      color = c("black","red","dodgerblue3"),
                      names.size = 2,
                      show.names = "highlighted",
                      highlight = c("gene1","gene2"),
                      xlab = " Gene expression fold change (Log2)",
                      legend = "State",
                      title = "Volcano plot (Normal NT vs Tumor TP)",
                      width = 10)

图 2. TCGA volcano

3.3 FEA：功能富集分析

data(DEGsmatrix)  # DEGsmatrix 即 差异分析的结果
dataFEA <- FEA(DEGsmatrix = DEGsmatrix)

图3. dataFEA

这里如果按照作者给的绘图代码肯定是画不出来的，因为他写的源码有点小问题，但影响不大，稍修改一下图还是能画出来的。（更新后的Moonlight2R并没有解决这个问题）

#绘图
plotFEA(dataFEA = dataFEA, additionalFilename = "_exampleVignette", height = 20, width = 10)

plotFEA
function (dataFEA, topBP = 10, additionalFilename = NULL, height, 
          width, offsetValue = 5, angle = 90, xleg = 35, yleg = 5, 
          titleMain, minY = -5, maxY = 10, mycols = c("#8DD3C7", "#FFFFB3", 
                                                      "#BEBADA")) 
{
  titleMain <- "TCGA BRCA DEGs"  ## 这里写的也是有问题的，基本函数自带的titleMain几乎无用
  if (!is.null(additionalFilename)) {
    pdf(additionalFilename, width, height)  
    ### 这里pdf()不完整，因此图画不出来，需要修改成：
    ### pdf(file = paste0("plotFEA", additionalFilename, ".pdf"))
  }
  ……
  if (!is.null(additionalFilename)) {
    dev.off()
  }
}
# 或者代码修改成，这样才能输出pdf文件。
plotFEA(dataFEA = dataFEA, additionalFilename = "_exampleVignette.pdf", height = 20, width = 10)

图4. FEAplot

3.4 GRN：基因调节网络分析

该步骤的核心是 ‘parmigene’ 包的 knnmi.cross 函数，本质是K-邻近聚类的算法

dataGRN <- GRN(TFs = rownames(DEGsmatrix)[1:100], normCounts = dataFilt,
                   nGenesPerm = 10,kNearest = 3,nBoot = 10)
# TFs  基因名，normCounts 表达矩阵

图5. dataGRN

在作者的说明文档里并没有提出该怎么解读这个dataGRN的结果，这里按下不表。对于GRN数据结果，作者与Cosmic数据库相关联，绘制巢状网络图，knownDriverGenes即来自于Cosmic 数据库

data(knownDriverGenes)
head(knownDriverGenes)
# $TSG
# [1] "ABI1"     "ALDH2"    "AMER1"    "APC"      "ARHGAP26" "ARHGEF12" "ARID1A"   "ARID1B"   "ARID2"    "ASXL1"    "ATM"      "BAP1"     "BCOR"    
# [14] "BRCA1"    "BRCA2"    "CASP8"    "CBLB"     "CDC73"    "CDKN1B"   "CDKN2A"   "CDKN2C"   "CHEK2"    "CTCF"     "CYLD"     "DICER1"   "DNMT3A"  
# [27] "ERCC2"    "FAT1"     "FAT4"     "IKZF1"    "LZTR1"    "MEN1"     "NCOR2"    "PALB2"    "PBRM1"    "PTEN"     "PTPN13"   "RB1"      "RBM10"   
# [40] "SDHA"     "SETD2"    "SH2B3"    "SMAD2"    "SMAD3"    "SMAD4"    "SOCS1"    "TET2"     "TGFBR2"   "TP53"     "TSC1"     "TSC2"     "VHL"     
# [53] "XPA"      "XPC"      "ZFHX3"   

# $OCG
# [1] "ABL1"    "ABL2"    "ACKR3"   "ACVR1"   "AFF1"    "AFF3"    "AFF4"    "AKAP9"   "AKT1"    "AKT2"    "ALK"     "AR"      "ASPSCR1" "ATP1A1" 
# [15] "BCL2"    "BCL3"    "BCL5"    "BCL6"    "BCL7A"   "BCL9"    "BCR"     "BRAF"    "CALR"    "CCND1"   "CDK4"    "CDK6"    "CHD4"    "CTNNB1" 
# [29] "CXCR4"   "ERBB2"   "ERBB3"   "EZR"     "FLT3"    "FOXP1"   "FUS"     "GATA2"   "HIF1A"   "HLF"     "HMGA2"   "HOXA11"  "HOXA13"  "HOXA9"  
# [43] "HOXC11"  "HOXC13"  "HOXD11"  "HOXD13"  "HRAS"    "IDH1"    "IDH2"    "JAK2"    "KIT"     "KLF4"    "KMT2A"   "KRAS"    "MDM2"    "MDM4"   
# [57] "MECOM"   "MET"     "MLLT1"   "MLLT3"   "MPL"     "MTOR"    "MYC"     "MYCN"    "MYD88"   "NAB2"    "NFKB2"   "NRAS"    "PDGFB"   "PDGFRA" 
# [71] "PIK3CA"  "PML"     "PRKACA"  "PTPN11"  "RARA"    "REL"     "RET"     "ROS1"    "RUNX1"   "RUNX1T1" "STAT3"   "SYK"     "TRRAP"   "USP8"   

data(dataGRN)
plotNetworkHive(dataGRN, knownDriverGenes, 0.55)

图6. plotNetworkHive

3.5 URA:上游调节子分析

这个步骤需要同时用到GRN和FEA的分析结果，

data(dataGRN)
data(DEGsmatrix)

dataFEA <- FEA(DEGsmatrix = DEGsmatrix)

BPselected <- dataFEA$Diseases.or.Functions.Annotation[1:5]
dataURA <- URA(dataGRN = dataGRN,
               DEGsmatrix = DEGsmatrix,
               BPname = BPselected,  ##可以是指定的Biological process
               nCores=1)
head(dataURA)
#       abdnominal adenocarcinoma adenocarcinoma alzheimers disease amyloidosis apoptosis of tumore cell lines
# A1BG                   0.0000000      0.0000000                  0   0.0000000                      0.6666667
# A2M                    0.0000000      0.0000000                  0   0.0000000                     -0.2294157
# AASS                   0.0000000      0.0000000                  0   0.0000000                     -0.3779645
# ABAT                   0.0000000      0.0000000                  0   0.0000000                      0.0000000
# ABCA10                -0.2041241     -0.1740777                  0   0.0000000                     -0.9486833
# ABCA12                 0.0000000      0.0000000                  0  -0.5303301                     -0.8340577

3.6 PRA: 模式识别子分析

data(dataURA)
dataDual <- PRA(dataURA = dataURA,
                          BPname = c("apoptosis","proliferation of cells"),
                          thres.role = 0)
head(dataDual)
#$TSG
#    ACN9    ADCK5  AFAP1L1 
#1.123390 1.126988 1.112131 

#$OCG
#    ABCA3     ABCG1     ABCG2     ABHD6     ABTB1     ACADL      ACAN 
#3.3825395 1.1088252 3.3062860 1.6521305 2.7699272 2.8953050 1.6175689 
#     ACE2      ACP5      ACP6      ACRC     ACSS1     ACSS2     ACTG2 
#1.1016235 1.4461591 1.7583342 1.2523922 2.8238690 0.6931391 2.9399313 
#   ADAM19     ADAM6  ADAMTS12  ADAMTS16    ADHFE1     ADRB1    ADSSL1 
#2.0317063 2.9133937 1.5595271 1.0147188 3.7288269 2.3278862 0.7482893 
#   AGPAT4 
#3.2790185

PRA嵌套了画图函数，具体如下

plotURA(dataURA = dataURA[c(names(dataDual$TSG),
                            names(dataDual$OCG)),, drop = FALSE], 
        additionalFilename = "_exampleVignette")

图7. plotURA

这里绘图得到的数值和作者提供在示例文档中的略有区别，但整体的趋势是一致的。

3.7 DMA：Driver Mutation Analysis（驱动突变分析）

正如前述，DMA分析嵌套了CScape-somatic 算法，这个部分的输入数据包括三个数据， MAF矩阵，差异分析结果（DEGmatrix）,以及3.6步骤中PRA的输出结果（dataPRA，含OCG）,另外以及CScape的类似于索引文件输入（这个部分为什么不作为R包的内置嵌入，而要人工输入？）。这里提供的文件索引在R包中没有找到，后续在CScape-somatic的网页中找到了这个数据，链接如下http://cscape-somatic.biocompute.org.uk/#download；

data(dataPRA)
data(dataMAF)
dataDMA <- DMA(dataMAF = dataMAF,
               dataDEGs = DEGsmatrix, 
               dataPRA = dataPRA,
               results_folder = "DMAresults",
               coding_file = "/data/databases/CScape/CScape-20210624/css_coding.vcf.gz",
               noncoding_file = "/data/databases/CScape/CScape-20210624/css_noncoding.vcf.gz")

稍稍看了下，这两个文件确实有点大，实在是包装不了。

这个部分主要输出是DEG_Mutations_Annotations Oncogenic_mediators_mutation_summary，Moonligh2R包装了画图函数plotMoonlight，具体使用如下

data(DEG_Mutations_Annotations)
data(Oncogenic_mediators_mutation_summary)
plotMoonlight(DEG_Mutations_Annotations,
        Oncogenic_mediators_mutation_summary, 
        dataURA, gene_type = "drivers", n = 50)

图8. plotMoonlight

以及plotDMA

data(dataDMA)
data(DEG_Mutations_Annotations)
data(Oncogenic_mediators_mutation_summary)
plotDMA(DEG_Mutations_Annotations, 
        Oncogenic_mediators_mutation_summary, 
        type = 'complete', additionalFilename = "./DMAresults/")

以上是Moonlight R语言包的主要函数。

4. 使用示例

作者在说明文档里提供了几个案例研究，主要包括以下三类

Downstream analysis LUAD using RNA expression data
Expression pipeline Pan Cancer with five cancer types
Expression pipeline with stages I II III IV (BRCA)

4.1 案例一 :Downstream analysis LUAD using RNA expression data

dataFilt <- getDataTCGA(cancerType = "LUAD", 
                          dataType = "Gene expression",
                          directory = "data",
                          nSample = 4)
dataDEGs <- DPA(dataFilt = dataFilt,
                dataType = "Gene expression")

dataFEA <- FEA(DEGsmatrix = dataDEGs)

dataGRN <- GRN(TFs = rownames(dataDEGs)[1:100], 
               DEGsmatrix = dataDEGs,
               DiffGenes = TRUE,
               normCounts = dataFilt)

dataURA <- URA(dataGRN = dataGRN, 
              DEGsmatrix = dataDEGs, 
              BPname = c("apoptosis",
                         "proliferation of cells"))

dataDual <- PRA(dataURA = dataURA, 
               BPname = c("apoptosis",
                          "proliferation of cells"),
               thres.role = 0)

CancerGenes <- list("TSG"=names(dataDual$TSG), "OCG"=names(dataDual$OCG))
head(CancerGenes)
# $TSG
# [1] "ACAT1"    "ACP5"     "ADAMTS14" "ADCY8"    "ADH1A"    "ADH1B"    "AGTR1"    "ALAS2"    "ALDH1A2"  "ANGPT1"   "ANKRD58"  "ANXA8"    "AP1S2"   
# [14] "APLNR"    "APLN"     "APOC1"    "ARC"      "ARHGAP32" "ARHGAP42" "ASPA"     "ATP1B2"   "ATP8A2"   "AVPR2"   

# $OCG
# [1] "ABCC2"    "ADAM12"   "ADAM28"   "ADAMTSL5" "ADH4"     "AMN"      "ANKFN1"   "ANKRD22"  "ANKRD43"  "ANLN"     "ANO3"     "APC2"     "APOB"    
# [14] "AQP11"    "ATHL1"    "AURKA"

CancerGenes 即MoonlightR 所获得的 TSG（肿瘤抑制基因），OCG(促癌基因）。

4.2 案例二：Expression pipeline Pan Cancer 5 cancer types

moonlight 函数是对所有函数的进行一个统一包装，但是，这里没有对dataMAF进行统一处理，新版本的moonlight2R只是简单的把dataMAF简单的写在函数的input里面，也没有对这个函数做获取数据的嵌入，导致他的示例数据无法使用。MoonlightR 版本是可以运行示例数据的，但Moonlight2R这里运行不了，就是因为dataMAF没有很好的嵌入。按照更新版本的moonlight，理论上来说，pan cancer 的分析代码中应该嵌入5个函数MAF文件的获取（可以参考源code的for 循环），但这里显然是没有的，感兴趣的可以看一下原始函数。

# MoonlightR version 1
cancerList <- c("BLCA","COAD","ESCA","HNSC","STAD")
## moonlight 嵌套了上述MoonlightR包的主要函数
listMoonlight <- moonlight(cancerType = cancerList, 
                      dataType = "Gene expression",
                      directory = "data",
                      nSample = 10,
                      nTF = 100,
                      DiffGenes = TRUE,
                      BPname = c("apoptosis","proliferation of cells"))
save(listMoonlight, file = paste0("listMoonlight_ncancer4.Rdata"))
## 绘图
plotCircos(listMoonlight = listMoonlight, additionalFilename = "_ncancer5")
# Moonlight2R
moonlight
function (cancerType = "panCancer", dataType = "Gene expression", 
          directory = "GDCdata", BPname = NULL, cor.cut = 0.6, qnt.cut = 0.25, 
          Genelist = NULL, fdr.cut = 0.01, logFC.cut = 1, corThreshold = 0.6, 
          kNearest = 3, nGenesPerm = 10, DiffGenes = FALSE, nBoot = 100, 
          nTF = NULL, nSample = NULL, thres.role = 0, stage = NULL, 
          subtype = 0, samples = NULL, dataMAF, path_cscape_coding,  # 仅简单增加dataMAF及相关输入文件
          path_cscape_noncoding) 
{……
  driverGenes <- DMA(dataMAF = dataMAF, dataDEGs = dataDEGs, 
                       dataPRA = listCandidates, runCscape = TRUE, coding_file = path_cscape_coding, 
                       noncoding_file = path_cscape_noncoding, results_folder = "./DMAresults")
 ……
 }

图10. listMoonlight 是多个流程结果的List

plotCircos 是对多个结果的 TSGs 和OCGs 的综合展示。

图11. plotCircos

4.3 案例三：Downstream analysis BRCA with stages

在MoonlightR初始版本中，出现了以下问题

rm(list=ls())

listMoonlight <- NULL
for (i in 1:4){
  dataDual <- moonlight(cancerType = "BRCA", 
                        dataType = "Gene expression",
                        directory = "data",
                       nSample = 10,
                        nTF = 5,
                        DiffGenes = TRUE,
                        BPname = c("apoptosis","proliferation of cells"),
                        stage = i)
  listMoonlight <- c(listMoonlight, list(dataDual))
  save(dataDual, file = paste0("dataDual_stage",as.roman(i), ".Rdata"))
}
# [1] "-----------------------------------------"
# [1] "Get TCGA data"
# [1] "-----------------------------------------"
# [1] "cancer type: BRCA"
# --------------------------------------
# o GDCquery: Searching in GDC database
# --------------------------------------
# Genome of reference: hg19
# --------------------------------------------
# oo Accessing GDC. This might take a while...
# --------------------------------------------
# ooo Project: TCGA-BRCA
# Connected to your session in progress, last started 2022-Oct-25 10:13:22 UTC (14 minutes ago)
# --------------------
# oo Filtering results
# --------------------
# ooo By platform
# ooo By data.type
# ooo By file.type
----------------
# oo Checking data
# ----------------
# ooo Checking if there are duplicated cases
# Warning: There are more than one file for the same case. Please verify query results. You can use the command View(getResults(query)) in rstudio
# ooo Checking if there are results for the query
# -------------------
# o Preparing output
# -------------------
#Error in `$<-.data.frame`(`*tmp*`, "tumor_stage", value = character(0)) : 
#  replacement has 0 rows, data has 1098

按作者提供的代码，出现报错，复盘→ Debug ，


moonlight
function (cancerType = "panCancer", dataType = "Gene expression", 
    directory = "GDCdata", BPname = NULL, cor.cut = 0.6, 
    qnt.cut = 0.25, Genelist = NULL, fdr.cut = 0.01, logFC.cut = 1, 
    corThreshold = 0.6, kNearest = 3, nGenesPerm = 10, DiffGenes = FALSE, 
    nBoot = 100, nTF = NULL, nSample = NULL, thres.role = 0, 
    stage = NULL, subtype = 0, samples = NULL) 
{
    GDCprojects <- get("GDCprojects")
    if (length(cancerType) == 1 && cancerType == "panCancer") {
        cancerType <- sort(sapply(strsplit(grep("TCGA", 
            GDCprojects, value = TRUE), "TCGA-"), "[", 
            2))
    }
    res <- NULL
    for (cancer.i in cancerType) {
        print("-----------------------------------------")
        print("Get TCGA data")
        print("-----------------------------------------")
        print(paste("cancer type:", cancer.i))
        dataFilt <- getDataTCGA(cancerType = cancer.i, dataType = dataType, 
            directory = directory, cor.cut = cor.cut, qnt.cut = qnt.cut, 
            nSample = nSample, stage = stage, subtype = subtype, 
            samples = samples)
        print("-----------------------------------------")
        print("Differential phenotype analysis")
        print("-----------------------------------------")
        dataDEGs <- DPA(dataType = dataType, dataFilt = dataFilt, 
            fdr.cut = fdr.cut, logFC.cut = logFC.cut)
        print("-----------------------------------------")
        print("Gene regulatory network")
        print("-----------------------------------------")
        if (is.null(nTF)) {
            nTF <- nrow(dataDEGs)
        }
        if (is.null(Genelist)) {
            Genelist <- rownames(dataDEGs)[1:nTF]
        }
        dataGRN <- GRN(TFs = Genelist, normCounts = dataFilt, 
            DEGsmatrix = dataDEGs, DiffGenes = FALSE, nGenesPerm = nGenesPerm, 
            kNearest = kNearest, nBoot = nBoot)
        print("-----------------------------------------")
        print("Upstream regulator analysis")
        print("-----------------------------------------")
        dataURA <- URA(dataGRN = dataGRN, DEGsmatrix = dataDEGs, 
            BPname = BPname)
        print("-----------------------------------------")
        print("Get candidates")
        print("-----------------------------------------")
        listCandidates <- PRA(dataURA = dataURA, BPname = BPname, 
            thres.role = thres.role)
        res.i <- list(dataDEGs = dataDEGs, dataURA = dataURA, 
            listCandidates = listCandidates)
        res <- c(res, list(res.i))
    }
    names(res) <- cancerType
    return(res)
}
<bytecode: 0x000002312986ce50>
<environment: namespace:MoonlightR>

仔细浏览及逐一运行moonlight函数的代码，发现 # o Preparing output 是 getDataTCGA 函数内部GDCquery函数的输出提示文本，而下一个提示文本Differential phenotype analysis 还没有出现表明，报错出现在getDataTCGA。继续在getDataTCGA中debug

library(TCGAbiolinks)
getDataTCGA
function (cancerType, dataType, directory, cor.cut = 0.6, qnt.cut = 0.25, 
          nSample, stage = "ALL", subtype = 0, samples = NULL) 
{
  DiseaseList <- get("DiseaseList")
  GDCprojects <- get("GDCprojects")
  geneInfo <- get("geneInfo")
  CancerProject <- paste0("TCGA-", cancerType)
  DataDirectory <- paste0(directory, "GDC_", gsub("-", 
                                                  "_", CancerProject))
  if (dataType == "Gene expression") {
    FileNameData <- paste0(DataDirectory, "_stage_", 
                           stage, "_subtype_", subtype, "_", "IlluminaHiSeq", 
                           ".rda")
    if (!file.exists(FileNameData)) {
      query <- GDCquery(project = CancerProject, data.category = "Gene expression", 
                        data.type = "Gene expression quantification", 
                        platform = "Illumina HiSeq", file.type = "results", 
                        legacy = TRUE)
      samplesDown <- query$results[[1]]$cases
      dataSmTP <- TCGAquery_SampleTypes(barcode = samplesDown, 
                                        typesample = "TP")
      dataSmNT <- TCGAquery_SampleTypes(barcode = samplesDown, 
                                        typesample = "NT")

if (is.numeric(stage) == TRUE) {
+         dataClin <- GDCquery_clinic(project = CancerProject, 
+                                     type = "clinical_patient")
+         curStage <- paste0("Stage ", as.roman(stage))
+         dataClin$tumor_stage <- toupper(dataClin$tumor_stage)
+         dataClin$tumor_stage <- gsub("[ABCDEFGH]", 
+                                      "", dataClin$tumor_stage)
+         dataClin$tumor_stage <- gsub("ST", "Stage", 
+                                      dataClin$tumor_stage)
+         dataStg <- dataClin[dataClin$tumor_stage %in% 
+                               curStage, ]
+         message(paste(curStage, "with", nrow(dataStg), 
+                       "samples"))
+         dataStgC <- dataSmTP[substr(dataSmTP, 1, 12) %in% 
+                                dataStg$bcr_patient_barcode]
+         dataSmTP <- dataStgC
+       }
##  类似的报错在这里出现了 !!!
## Error in `$<-.data.frame`(`*tmp*`, tumor_stage, value = character(0)) :  。
## 替换数据里有0行，但数据有1098
#   继续逐一运行，发现报错主要来自这一行
dataClin$tumor_stage <- toupper(dataClin$tumor_stage)
###   Error in `$<-.data.frame`(`*tmp*`, tumor_stage, value = character(0)) : 
###   替换数据里有0行，但数据有1098
dataClin$tumor_stage
###  NULL

图12. dataClin

图13. dataClin

由此，Bug 主要来自于BRCA的数据 dataClin 没有 tumor_stage 这一列，而是ajcc_pathologic_stage 列。而后逐一查看“COAD”及‘BLCA’的临床数据集均出现问题。所以这部分的代码大概率是不太行了。
在Moonlight2R的更新版本中，moonlight函数并没有提供获取MAF矩阵的函数，因此这部分报错是预料之中的，但是让我意外的是，更新后的Moonlight 还是出现了和之前差不多的报错，这……作者上传之前都不检查一下自己写的代码能不能运行的嘛？？

> for (i in 1:4){
+   dataDual <- moonlight(cancerType = "BRCA", 
+                         dataType = "Gene expression",
+                         directory = "data",
+                         nSample = 10,
+                         nTF = 5,
+                         DiffGenes = TRUE,
+                         BPname = c("apoptosis","proliferation of cells"),
+                         stage = i)
+   listMoonlight <- c(listMoonlight, list(dataDual))
+   save(dataDual, file = paste0("dataDual_stage",as.roman(i), ".Rdata"))
+ }
[1] "-----------------------------------------"
[1] "Get TCGA data"
[1] "-----------------------------------------"
[1] "cancer type: BRCA"
--------------------------------------
o GDCquery: Searching in GDC database
--------------------------------------
Genome of reference: hg19
--------------------------------------------
oo Accessing GDC. This might take a while...
--------------------------------------------
ooo Project: TCGA-BRCA
--------------------
oo Filtering results
--------------------
ooo By platform
ooo By data.type
ooo By file.type
----------------
oo Checking data
----------------
ooo Checking if there are duplicated cases
Warning: There are more than one file for the same case. Please verify query results. You can use the command View(getResults(query)) in rstudio
ooo Checking if there are results for the query
-------------------
o Preparing output
-------------------
Error in `$<-.data.frame`(`*tmp*`, "tumor_stage", value = character(0)) : 
  replacement has 0 rows, data has 1098
In addition: Warning message:
In readBin(3L, raw(0), 32768L) :
  URL 'https://api.gdc.cancer.gov/legacy/files/?pretty=true&expand=cases,cases.samples.portions.analytes.aliquots,cases.project,center,analysis,cases.samples&size=15492&filters=%7B%22op%22:%22and%22,%22content%22:[%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22cases.project.project_id%22,%22value%22:[%22TCGA-BRCA%22]%7D%7D,%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22files.data_category%22,%22value%22:[%22Gene%20expression%22]%7D%7D,%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22files.data_type%22,%22value%22:[%22Gene%20expression%20quantification%22]%7D%7D,%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22files.platform%22,%22value%22:[%22Illumina%20HiSeq%22]%7D%7D,%7B%22op%22:%22in%22,%22content%22:%7B%22field%22:%22files.tags%22,%22value%22:[%22unnormalized%22]%7D%7D]%7D&format=JSON': Timeout of 60 seconds was reached

5. 总结

总体而言，MoonlightR 确实是提供了一个推断肿瘤抑制基因（TSG）或者促癌基因（OCG）的一个R包，后续新近更新版本的Moonlight2R本以为会修复一些Bug，但没想作者只是简单的将新函数加入，而没有对原先的Code 进行检查，并且新函数加入的也很粗糙，这就导致总体的体验是有点差的，或许作者也没有想到自己写的R包也会有人学习？实事求是的说，这篇示例文档写的并不太友好，更新版本的Moonlight2R尽管功能更加强大，用也是能用的，但用的就是很不爽，总体观感就不是太好。抛开代码及文档的因素不谈，这个代码流程确实很有意义，有值得学习的地方，也能够用于课题相关肿瘤驱动基因的挖掘中，但是对使用者的能力有更高的要求，能够成功Debug作者的不尽之处，给出正确的解决办法，也是对使用者的一个小考验。

- END -

ixxmu / mp_duty

Moonlight：一种识别生物标志物在不同肿瘤类型和分期中作为癌基因或肿瘤抑制因子的多种作用的方法 #2988