# MEGENA_pipeline

#### Won-Min Song

#### January 04, 2017

## MEGENA pipeline as of 10/25/2016

This is a routine MEGENA pipeline description encompassing from data correlation analysis to module plotting, and is based on version 1.3.6 https://CRAN.R-project.org/package=MEGENA.Please cite the paper below when MEGENA is applied as part of your analysis:

Song W.-M., Zhang B. (2015) Multiscale Embedded Gene Co-expression Network Analysis. PLoS Comput Biol 11(11): e1004574. doi: 10.1371/journal.pcbi.1004574.

For statistical mechanics aspects involved in MEGENA, please check:

Song W.-M., Di Matteo T.and Aste T., Building Complex Networks with Platonic Solids, Physical Reivew E, 2012 Apr;85(4 Pt 2):046115.

# calculate correlation

    rm(list = ls()) # rm R working space

    library(MEGENA)

    # input parameters
    n.cores <- 2; # number of cores/threads to call for PCP
    doPar <-TRUE; # do we want to parallelize?
    method = "pearson" # method for correlation. either pearson or spearman.
    FDR.cutoff = 0.05 # FDR threshold to define significant correlations upon shuffling samples.
    module.pval = 0.05 # module significance p-value. Recommended is 0.05.
    hub.pval = 0.05 # connectivity significance p-value based random tetrahedral networks
    cor.perm = 10; # number of permutations for calculating FDRs for all correlation pairs.
    hub.perm = 100; # number of permutations for calculating connectivity significance p-value.

    # annotation to be done on the downstream
    annot.table=NULL
    id.col = 1
    symbol.col= 2
    ###########

    data(Sample_Expression) # load toy example data

    ijw <- calculate.correlation(datExpr,doPerm = cor.perm,output.corTable = FALSE,output.permFDR = FALSE)

    ## i = 1
    ## i = 2
    ## i = 3
    ## i = 4
    ## i = 5
    ## i = 6
    ## i = 7
    ## i = 8
    ## i = 9
    ## i = 10

# calculate PFN

In this step, Planar Filtered Network (PFN) is calculated by taking significant correlation pairs, ijw. In the case of utilizing a different similarity measure, one can independently format the results into 3-column data frame with column names c(“row”,“col”,“weight”), and make sure the weight column ranges within 0 to 1. Using this as an input to calculate.PFN() will work just as fine.

    #### register multiple cores if needed: note that set.parallel.backend() is deprecated.
    if (doPar & getDoParWorkers() == 1)
    {
      cl <- parallel::makeCluster(n.cores)
      registerDoParallel(cl)
      # check how many workers are there
      cat(paste("number of cores to use:",getDoParWorkers(),"\n",sep = ""))
    }

    ## number of cores to use:2

    ##### calculate PFN
    el <- calculate.PFN(ijw[,1:3],doPar = doPar,num.cores = n.cores,keep.track = FALSE)

    ## ####### PFN Calculation commences ########
    ## [1] "343 out of 894 has been included."
    ## [1] "471 out of 894 has been included."
    ## [1] "556 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 249"
    ## [1] "663 out of 894 has been included."
    ## [1] "710 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 78"
    ## [1] "759 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 60"
    ## [1] "802 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 24"
    ## [1] "820 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 21"
    ## [1] "839 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 14"
    ## [1] "850 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 5"
    ## [1] "853 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 4"
    ## [1] "857 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 3"
    ## [1] "860 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 3"
    ## [1] "863 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 5"
    ## [1] "867 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 1"
    ## [1] "868 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 4"
    ## [1] "872 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 10"
    ## [1] "875 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 1"
    ## [1] "876 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 2"
    ## [1] "878 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 5"
    ## [1] "881 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 5"
    ## [1] "884 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 19"
    ## [1] "888 out of 894 has been included."
    ## [1] "Performing parallel quality checks on 2000"
    ## [1] "Qualified edges: 217"
    ## [1] "894 out of 894 has been included."
    ## [1] "PFG is complete."

    g <- graph.data.frame(el,directed = FALSE)

# perform clustering

MCA clustering is performed to identify multiscale clustering analysis. “MEGENA.output”" is the core output to be used in the down-stream analyses for summarization and plotting.

    ##### perform MCA clustering.
    MEGENA.output <- do.MEGENA(g,
     mod.pval = module.pval,hub.pval = hub.pval,remove.unsig = TRUE,
     min.size = 10,max.size = vcount(g)/2,
     doPar = doPar,num.cores = n.cores,n.perm = hub.perm,
     save.output = FALSE)

    ###### unregister cores as these are not needed anymore.
    if (getDoParWorkers() > 1)
    {
      env <- foreach:::.foreachGlobals
      rm(list=ls(name=env), pos=env)
    }

# summarize results

    summary.output <- MEGENA.ModuleSummary(MEGENA.output,
        mod.pvalue = module.pval,hub.pvalue = hub.pval,
        min.size = 10,max.size = vcount(g)/2,
        annot.table = annot.table,id.col = id.col,symbol.col = symbol.col,
        output.sig = TRUE)

    if (!is.null(annot.table))
    {
      # update annotation to map to gene symbols
      V(g)$name <- paste(annot.table[[symbol.col]][match(V(g)$name,annot.table[[id.col]])],V(g)$name,sep = "|")
      summary.output <- output[c("mapped.modules","module.table")]
      names(summary.output)[1] <- "modules"
    }

    print(head(summary.output$modules,2))

    ## $c1_2
    ##   [1] "C1QB"         "C1QA"         "SASH3"        "HLA-DPA1"
    ##   [5] "PLEK"         "NCKAP1L"      "EVI2B"        "HLA-DRA"
    ##   [9] "FCGR1A"       "CD53"         "PTPRC"        "DOCK2"
    ##  [13] "ITGB2"        "FCGR1C"       "ARHGAP9"      "CD37"
    ##  [17] "CD4"          "PTPN7"        "HLA-DMB"      "CD48"
    ##  [21] "FYB"          "CYBB"         "IKZF1"        "LILRB2"
    ##  [25] "HLA-DQA1"     "LILRB4"       "IL10RA"       "CCR5"
    ##  [29] "HLA-DPB1"     "LILRB1"       "CORO1A"       "SPN"
    ##  [33] "PSTPIP1"      "SNX20"        "MNDA"         "TFEC"
    ##  [37] "FGL2"         "CD74"         "C1QC"         "ITGAX"
    ##  [41] "FERMT3"       "CCR1"         "PRKCB"        "CD1E"
    ##  [45] "HLA-DQB1"     "IL12RB1"      "NCF1C"        "CLEC10A"
    ##  [49] "FPR3"         "WDFY4"        "SP140"        "FCGR3A"
    ##  [53] "CTSS"         "BIN2"         "SIGLEC7"      "TLR8"
    ##  [57] "HLA-DRB1"     "IFI30"        "GPR183"       "NCF2"
    ##  [61] "TMC8"         "MRC1"         "GIMAP5"       "CLEC17A"
    ##  [65] "CD86"         "PIK3AP1"      "CSF1R"        "CD163"
    ##  [69] "P2RY13"       "DOK2"         "CARD11"       "LILRB3"
    ##  [73] "APOC1"        "IRF8"         "CRTAM"        "PARP15"
    ##  [77] "AMICA1"       "SRGN"         "TNFRSF8"      "PLD4"
    ##  [81] "MPEG1"        "CD1A"         "TMEM176A"     "EMB"
    ##  [85] "RCSD1"        "FCGR1B"       "HCLS1"        "HLA-DOA"
    ##  [89] "CIITA"        "AOAH"         "SLA"          "KLHL6"
    ##  [93] "PIK3CG"       "CD84"         "SIGLEC10"     "CD1C"
    ##  [97] "NCF1"         "CD226"        "PTPN22"       "IGSF6"
    ## [101] "ZBP1"         "IL16"         "LYZ"          "HLA-DRB5"
    ## [105] "CCR4"         "CD209"        "CD200R1"      "TLR10"
    ## [109] "SIGLEC1"      "LOC100233209" "LILRA5"       "MYO1G"
    ## [113] "CD28"         "CYTIP"        "PTAFR"        "HCK"
    ## [117] "CSF1"         "TLR7"         "RAC2"         "ITGAL"
    ## [121] "STAP1"        "RASGRP2"      "PTPRO"        "CD180"
    ## [125] "P2RY12"       "EBI3"         "RASSF4"       "SIGLEC5"
    ## [129] "FCN1"         "OSCAR"        "OLR1"         "FCGR2B"
    ## [133] "NLRP2"
    ##
    ## $c1_3
    ##   [1] "CD3E"      "CD3D"      "CD2"       "CXCL11"    "ITK"
    ##   [6] "ACAP1"     "IL2RG"     "PTPRCAP"   "SLAMF6"    "SH2D1A"
    ##  [11] "SIRPG"     "CD5"       "NKG7"      "TIGIT"     "CD96"
    ##  [16] "CD247"     "CXCR3"     "ICOS"      "LY9"       "TBC1D10C"
    ##  [21] "UBASH3A"   "TAP1"      "LCK"       "GZMA"      "TBX21"
    ##  [26] "GPR174"    "CCL5"      "CD38"      "CXCL10"    "C5orf20"
    ##  [31] "IRF4"      "ZNF831"    "GBP5"      "LAX1"      "CD27"
    ##  [36] "SLA2"      "KLRK1"     "GZMB"      "BTLA"      "HLA-B"
    ##  [41] "SCML4"     "HLA-F"     "GZMK"      "HLA-H"     "SLAMF1"
    ##  [46] "LOC400759" "SLAMF7"    "ZBED2"     "PRF1"      "IDO1"
    ##  [51] "CD8B"      "IL4I1"     "ZAP70"     "GVIN1"     "FCRL3"
    ##  [56] "FASLG"     "APOL1"     "PSMB9"     "AQP9"      "TNFRSF9"
    ##  [61] "CXCL9"     "C8orf80"   "GBP1"      "HK3"       "ADAMDEC1"
    ##  [66] "CD80"      "CTLA4"     "SLAMF8"    "TIFAB"     "CLEC4D"
    ##  [71] "GPR18"     "SELL"      "APOL3"     "CLEC4E"    "PRKCQ"
    ##  [76] "CARD16"    "CCL4"      "KLRC3"     "IL7R"      "SIT1"
    ##  [81] "TRAT1"     "CD3G"      "GPR171"    "CD6"       "C16orf54"
    ##  [86] "CD8A"      "CD40LG"    "EOMES"     "GZMM"      "LTA"
    ##  [91] "MAP4K1"    "KIAA0748"  "SEPT1"     "CD52"      "PLA2G2D"
    ##  [96] "HCP5"      "MEI1"      "CST7"      "LTB"       "XCL2"
    ## [101] "AIM2"      "NLRC5"     "GNLY"      "IFNG"      "ETV7"
    ## [106] "CTSW"      "ZNF683"    "C11orf21"  "CD69"      "PNOC"
    ## [111] "IL2RB"     "TMEM150B"  "CCR8"      "LGALS2"    "LRMP"
    ## [116] "BCL11B"    "IKZF3"     "UBD"       "MCART6"    "CCL3L1"
    ## [121] "C15orf56"

    print(summary.output$module.table)

    ##       module.id module.size module.parent
    ## c1_2       c1_2         133          c1_1
    ## c1_3       c1_3         121          c1_1
    ## c1_4       c1_4          22          c1_1
    ## c1_5       c1_5          24          c1_1
    ## c1_6       c1_6         123          c1_2
    ## c1_8       c1_8          94          c1_3
    ## c1_10     c1_10          95          c1_6
    ##                                                                          module.hub
    ## c1_2  SASH3(34),CD53(33),PTPRC(26),CD86(20),NCKAP1L(19),FERMT3(19),PLEK(18),CD4(16)
    ## c1_3                                   CD2(33),CD3E(32),GBP5(20),ITK(19),SLAMF6(18)
    ## c1_4                                                                      MS4A1(12)
    ## c1_5                                                                      IFIT3(19)
    ## c1_6  SASH3(34),CD53(33),PTPRC(26),NCKAP1L(19),FERMT3(19),PLEK(18),CD86(17),CD4(16)
    ## c1_8                                            CD3E(32),CD2(32),ITK(19),SLAMF6(18)
    ## c1_10                   CD53(33),SASH3(29),FERMT3(19),PLEK(18),NCKAP1L(18),CD86(17)
    ##       module.scale module.pvalue
    ## c1_2            S2             0
    ## c1_3          <NA>             0
    ## c1_4            S2             0
    ## c1_5            S2             0
    ## c1_6          <NA>             0
    ## c1_8            S2             0
    ## c1_10           S1             0

# Plot some modules

You can generate refined module network plots:

    pnet.obj <- plot_module(output.summary = summary.output,PFN = g,subset.module = "c1_3",
        layout = "kamada.kawai",label.hubs.only = TRUE,
        gene.set = NULL,color.code =  "grey",
        output.plot = FALSE,out.dir = "modulePlot",col.names = c("magenta","green","cyan"),label.scaleFactor = 20,
        hubLabel.col = "black",hubLabel.sizeProp = 1,show.topn.hubs = Inf,show.legend = TRUE)

    ## Processing module:c1_3
    ##  - # of genes: 121
    ##  - # of hubs: 5
    ## - generating module subnetwork figure...

    #X11();
    print(pnet.obj[[1]])

# plot module hierarchy

    module.table <- summary.output$module.table
    colnames(module.table)[1] <- "id" # first column of module table must be labelled as "id".

    hierarchy.obj <- plot_module_hierarchy(module.table = module.table,label.scaleFactor = 0.15,
                                        arrow.size = 0.03,node.label.color = "blue")
    #X11();
    print(hierarchy.obj[[1]])