RESEARCHARTICLE Gene Network Rewiring to Study Melanoma Stage Progression and Elements Essential for Driving Melanoma AbhinavKaushik1,YashumaBhatia1,ShakirAli2,DineshGupta1* 1 BioinformaticsLaboratory,StructuralandComputationalBiologyGroup,InternationalCentreforGenetic EngineeringandBiotechnology,NewDelhi,110067,India,2 DepartmentofBiochemistry,Facultyof Science,JamiaHamdard,NewDelhi,110062,India *[email protected] Abstract Metastaticmelanomapatientshaveapoorprognosis,mainlyattributabletotheunderlying heterogeneityinmelanomadrivergenesandalteredgeneexpressionprofiles.Thesechar- acteristicsofmelanomaalsomakethedevelopmentofdrugsandidentificationofnovel drugtargetsformetastaticmelanomaadauntingtask.Systemsbiologyoffersanalternative OPENACCESS approachtore-explorethegenesorgenesetsthatdisplaydysregulatedbehaviourwithout Citation:KaushikA,BhatiaY,AliS,GuptaD(2015) beingdifferentiallyexpressed.Inthisstudy,wehaveperformedsystemsbiologystudiesto GeneNetworkRewiringtoStudyMelanomaStage ProgressionandElementsEssentialforDriving enhanceourknowledgeabouttheconservedpropertyofdiseasegenesorgenesets Melanoma.PLoSONE10(11):e0142443. amongmutuallyexclusivedatasetsrepresentingmelanomaprogression.Wemeta-ana- doi:10.1371/journal.pone.0142443 lysed642microarraysamplestogeneratemelanomareconstructednetworksrepresenting Editor:RogerChammas,UniversidadedeSão fourdifferentstagesofmelanomaprogressiontoextractgeneswithalteredmolecularcir- Paulo,BRAZIL cuitrywiringascomparedtoanormalcellularstate.Intriguingly,amajorityofthemelanoma Received:July20,2015 network-rewiredgenesarenotdifferentiallyexpressedandthediseasegenesinvolvedin Accepted:October21,2015 melanomaprogressionconsistentlymodulateitsactivitybyrewiringnetworkconnections. Wefoundthattheshortlisteddiseasegenesinthestudyshowstrongandabnormalnetwork Published:November11,2015 connectivity,whichenhanceswiththediseaseprogression.Moreover,thedeviatednetwork Copyright:©2015Kaushiketal.Thisisanopen propertiesofthediseasegenesetsallowranking/prioritizationofdifferentenriched,dysre- accessarticledistributedunderthetermsofthe CreativeCommonsAttributionLicense,whichpermits gulatedandconservedpathwaytermsinmetastaticmelanoma,inagreementwithprevious unrestricteduse,distribution,andreproductioninany findings.Ouranalysisalsorevealspresenceofdistinctnetworkhubsindifferentstagesof medium,providedtheoriginalauthorandsourceare metastasizingtumorforthesamesetofpathwaysinthestatisticallyconservedgenesets. credited. Thestudyresultsarealsopresentedasafreelyavailabledatabaseathttp://bioinfo.icgeb. DataAvailabilityStatement:ProcessedDataare res.in/m3db/.Theweb-baseddatabaseresourceconsistsofresultsfromtheanalysispre- availableathttp://dx.doi.org/10.6084/m9.figshare. sentedhere,integratedwithcytoscapewebanduser-friendlytoolsforvisualization,retrieval 1577546. andfurtheranalysis. Funding:DGconceivedandreceivedagrantfrom theDepartmentofBiotechnology(DBT)-Indiagrants (Grantno.BT/BI/25/066/2012).AKreceivedresearch fellowshipfromCouncilofScientificandIndustrial Research(CSIR)-India. CompetingInterests:Theauthorshavedeclared thatnocompetinginterestsexist. PLOSONE|DOI:10.1371/journal.pone.0142443 November11,2015 1/30 GeneNetworkRewiringtoStudyMelanomaStageProgression Introduction Advancemalignantmelanomarepresentsadeadlycancerstateduetoitshighdissemination potentialandincreasingtherapyresistance[1].Ingeneral,melanomainitiateswithmarkeddis- ruptionofcellularhomeostaticmechanismsleadingtoamalignanttransformationofskin melanocytes(cutaneousnon-metastatic;CnM).Melanomainitiationisfollowedbyitsprolifer- ationtodifferentlayersofskinviamultiplegrowthphases(cutaneousmetastasis;CM)and finallytolymphnodes(LN),whichmetastasizetumourtodifferentorgans[2]. Aftertheadventofthe‘omics’era,despiteseveralattemptstoidentifygenesignaturesand geneticchangesinmelanoma,inconsistenciesintheresultsstillexists[3–5].Thelikelyreason fortheinconsistencyisthediversehistologicaltypesandmultivariatenatureofmelanoma, makingconsensusviewaboutmelanomageneexpressionbehaviour,achallenge,whichalso hindersdevelopmentofnovelanti-melanomatherapies. Historically,high-throughput(HT)melanomastudieshavebeenmainlyconfinedtoidenti- ficationandvalidationofmutatedgenesandmultitudeofDifferentiallyExpressed(DE)genes forvariousclassifiedmelanomasormelanomasub-types.AlthoughtheinformationaboutDE genesenhanceourunderstandingofdiseases,recentstudiesindicatethatalterationingene activityisnotlimitedtoitsfoldchangeinexpressionaloneandtheinteractomesundergore- wiringasaresultofcellularoradaptiveresponse[6].Systemsbiologyoffersarangeofalterna- tivemethodsthatallowsidentificationofsuchdiseaselinkedgenescharacterisedbyabnormal linkedbehaviour.Thenewsystemsbiologymethodsbasedongeneco-expressionvaluesallows predictionofgeneessentialitybeyondmodulesandhubs;whichevenhelpsgeneratenovel hypothesesfordiseasemechanisms[7]. Ithasbeenfoundthatinanormalcellularhomeostasis,thegenesdisplayarobustsetofcor- relatedoscillations,hereinreferredasnetworkconnections,whicharedisruptedindiseases[7, 8].Thus,dysregulationlinkedwithdiseasestatescausedifferentialcorrelationofvariousgenes duetotheiralteredgeneexpressionpatternleadingtoreformationofnetworkwiringcircuitry [9].Itisexpectedthatthegenesthatamplifyitsconnectivityorexpressionindiseaseascom- paredtothatincontrolstagesaremorelikelytobeinvolvedindiseaseprogression.Insuch stressconditions,severalgenesloseconnectionsintheperturbednetworkwhilesomemay evengainalternatesetsofconnections[10].Fig1illustratesthedifferentialconnectivityinnet- work-1ascomparedtonetwork-2,whichshowsthatsomeconnectionsarelostwhilenewcon- nections(evennewnodes)maybeestablishedinthealterednetwork-2.Additional experimentalevidenceforthisconceptisareportbyHudsonetal.[11],whoconducteddiffer- entialwiringanalysisofexpressiondatatocorrectlypredictmyostatinasthegenecontaining thecausalmutation,characterisedbydifferentialwiringratherthandifferentialexpression. Anglanietal.describedthelossofgeneconnectionsasameasureofitsmoleculardysfunction infivedifferentcancers,whichoccurswhengenesloseitsroleasaneffector[12].Thestudy highlightedthatoncogenicalterationsincodingregionsandpost-translationalmodifications canmodifygeneconnectionsandhencefunctions,withoutanychangesinexpressionlevels.In anotherstudy,Kimetal.correlatedthegeneessentialitywithrewiringinyeastandmouse interactomeandsuggestedthatcircuitryrewiringallowgenestointegrateintoalternativepath- ways[13].Asimilarstudyutilizesthenumberofconnectionsgainedandlostbyagenetomea- sureitsoverallrewiringscore[8].Theanalysisrevealedthatthediseasegenestendtorewireits networkcircuitandhencegene-rewiringscoresaresuitableforgeneprioritization.Moreover, previouslyitwasfoundthatERK-MAPKactivityinmelanomaleadstopathwayrewiring whichresultsinitsconnectionofERKtoJNKpathways-thisobservationleadtoidentification ofadditionalmelanomatherapeutictargets[14].Thus,theavailableliteraturesuggeststhatnot onlygenesbuttheassociatedpathwaystoogetdysregulatedafterconnectionrewiring, PLOSONE|DOI:10.1371/journal.pone.0142443 November11,2015 2/30 GeneNetworkRewiringtoStudyMelanomaStageProgression Fig1.Geneconnectionsinbiologicalnetworksaredynamicandmayshowalteredco-expression undervariableconditions,evenifgenesshowsnosignificantdifferentialco-expression.Gene connectionsaremarkedaslines,newgenesmarkedashollowcircles(markedwitharrows)anddottedlines representnewconnections. doi:10.1371/journal.pone.0142443.g001 understandingofwhichisacriticalstepindesigningnovelandmoreeffectivetherapeutics [15].Anadditionaladvantageofexploitingthenetworkbasedsystembiologyapproachforthe predictionofgenesinvolvedindiseaseisthatthesensitivityofreconstructednetworkincreases withheterogeneityinexpressiondataset,whichisacommonfeatureofmelanomasamples[5, 16]. Thereisnoreportofanystudy,whichattemptstoidentifyandmeasurethenetworkdepen- dentdynamicityandreproducibilityofsuchpotentialandnoveldiseasegenesorgenesets(i.e. pathways)involvedinmelanomaprogression,basedontheanalysisofheterogeneousmela- nomadatasets.Inthisstudyweexploitedthepubliclyavailablecross-platformmelanomatran- scriptomicdata(samplesize>600,groupedintofourdifferentstagespersemelanoma progression),todecipherdifferentiallywiredgenesineachstage. Inthisstudywehaveidentifiedgenesetsthatarerewiredandconnectedsimilarlyinmutu- allyexclusivedatasets.Moreover,wehaveinvestigatedthedifferentialnetworkpropertiesof functionalclustersdevelopedfromtheshortlisteddiseasegenes.Thedeviatednetworkprofile analysisassistedinidentificationofimportantstagespecificpathwaytermsinvolvedinstage transition.Wedemonstratethatthereisasuddenincreaseinquantitativecomplexityofthe predictedpathwayswithtumourprogressionespeciallyaftertheonsetofmetastasis.Suchpri- oritized/rankedgeneclusterswerefoundtopossessstrongdiseaseassociationonthebasisof diseaselinkageanalysis.Wealsoattemptedtounderstandthedynamicityandflowofranked pathwaytermsenrichedduringmelanomaprogression.Thecompleteresultsarepresentedas aSQLdatabaseusingDHTMLandcytoscape-web[17]foruserfriendlyvisualization,retrieval andfurtheranalysis.Thedatabaseresourceisfreelyavailableathttp://bioinfo.icgeb.res.in/ m3db/. PLOSONE|DOI:10.1371/journal.pone.0142443 November11,2015 3/30 GeneNetworkRewiringtoStudyMelanomaStageProgression Table1. Datasetdescription. Melanomastage Datasetsize Series Platforms† Normal 97 13 07[a,b,c,d,e,f,g] Cutaneousnon-Metastatic 183 09 05[d,e,f,h,g] CutaneousMetastatic 246 10 05[d,e,f,h,g] LymphNodeMetastatic 116 07 05[a,c,d,e,h,i] †Microarrayplatformsusedbysubmitter. a:GPL10558 b:GPL4133 c:GPL5175 d:GPL570 e:GPL571 f:GPL6883 g:GPL96 h:GPL1708 i:GPL6884 doi:10.1371/journal.pone.0142443.t001 Results Melanomageneexpressiondatacollectedfrompubliclyavailabledatarepositorieswasmanu- allystratifiedintofourdifferentstages,intheorderofmelanomaprogressioninvolvingnormal skinmelanocytes(N)andthethreetransitionevents{N!CnM!CM!LN}.Table1shows thenumberofcross-platformgeneexpressiondatasetsizesandseries,representingeachstage anddifferentsources(fordetailsofeachofthestudydataset,seeS1File).Weperformedesti- mationandremovalofbatcheffectresultingfromdifferentexperimentalconditions,priorto thedataanalysis.TheresultsaresummarizedinS1Fig,whichsuggeststhatalargeproportion ofgeneswereaffectedwithbatcheffect,removedbyComBatbatchadjustment(seematerials andmethods). Meta-analysisrevealedanovelDEmelanomagenesignature Weinitiatedtheanalysisbycheckingthedatasetqualityandvalidationbypredictionofknown melanomamarkersamongststatisticallysignificantDEgenes(logFC>1andadj.p<0.05) betweencontrol(N)anddifferentstagesofmelanoma(CnM,CMandLN).Comparativeanal- ysishighlightsseveralstagespecificchanges(Fig2).Weidentified324,622and1398DEgenes inCnM,CMandLN,respectively.Weobservedthatdespitethelowfoldchangethreshold, only55genesareexpressedacrossallthestagesofmelanomaandindicatethatthegenesare essentiallyrequiredformelanomaprogression. Welookedformelanomabiomarkersandknowncancerrelatedgenesintheconserved geneset(n=55)andfoundseveraluniquesetsofgenespreviouslyreportedtobemelanomaor cancerlinked.Thegenesetisenrichedwithseveralmelanoma-associatedgenesincluding MAGEA1,MAGEA12,IL8,FOXD1,IL1B,POSTN,PRAME,MMP9,SERPINEandCTGF. Genesconservedbetweenmetastaticstages(n=67)tooholdsseveralkeyregulatorsofmela- nomalikeVEGFA,GDF15,SPP1,UPP1,SPRY4,FGF2,CENPN,SERPINA3,BUB1.Thecom- pletelistofbothgenesetsalongwithitsdescriptionisgivenintheS2File.Despite heterogeneityintheconditionsusedtogeneratethedatasets,themeta-analysisshowsthat thesefewbiomarkersandknowncancerlinkedgenesshowconsistentupregulation,whichnot onlystatisticallyrevalidatestheirsignificanceandessentialsimilaritiesinthedatasetbutalso revealsnovelmelanomaassociatedgenes. PLOSONE|DOI:10.1371/journal.pone.0142443 November11,2015 4/30 GeneNetworkRewiringtoStudyMelanomaStageProgression Fig2.A.Scatterplotoflog2ratio(foldchange;FC)versusadjustedp-valuetomeasuretheassociatedsignificantdifferentialexpressionofeach gene.Geneswhicharesignificantlydifferentiallyexpressed(DE;log(FC)>1or<1;adj.pvalue<0.05)areshownindarkbluecolorB.Significantly DEgenecount(up-anddown-regulated)ineachstage.C.Numberofup-regulatedgenesconservedacrossdifferentmelanomastages. doi:10.1371/journal.pone.0142443.g002 Tofurthervalidatethereliabilityofthestudydatasetsandanalysisresults,wecomparedthe foldchangeanalysisofDEgenesacrossdifferentstages(Fig2A).Theaveragefoldchanges observedforCnM,CMandLNare~0.4,~0.5and~0.7,respectively.Thetotalnumberofup- regulatedgenesinLNisalsomuchhigherthananyoftheotherstages,asevidentfromthe majordriftinfoldchangevolcanoplot,whereas,CnMcontainstheleastnumberofDEgenes. Foldchangeanalysissuggestsasignificantdown-regulationofgeneexpressionincutaneous melanomastages,especiallyduringdiseasemetastasis(S2Fig).CloserinvestigationoftheDE genelistforeachstagerevealedseveralmelanomabiomarkersandassociatedgenes.Forexam- ple,theMAGEfamilygenes(fewofwhichareconsideredasmelanomabiomarkers)andIL8(a well-knowncancerrelatedgene)areamongthetop-rankedgenesenrichedinoursignificantly up-regulatedgenelist[18–21].Table2highlightsthetop20DEup-regulatedgenesindecreas- ingorderoffoldchangeineachstage.Itisnotasurprisethataknownmelanomabiomarker genePRAME[22]isconsistentlyfoundtobetopmostup-regulatedgeneinallthemelanoma stages.AmongtheMAGEfamilygenes,MAGEA12isconsistentlyhighlyup-regulatedinall themelanomastages,whereasMAGEA1isup-regulatedonlyinthemetastaticmelanoma stages.CnMisalsocharacterisedbytheexpressionofseveralkeratinfamilygenesinvolvedin buildingstructuralframeworkofepithelialcellsalongwithS100familygenes,involvedin PLOSONE|DOI:10.1371/journal.pone.0142443 November11,2015 5/30 GeneNetworkRewiringtoStudyMelanomaStageProgression Table2. DifferentiallyExpressedGenes(inthedecreasingorderoffoldchange). CnM CM LN PRAME PRAME PRAME S100A7 SPP1 SPP1 S100A8 IL8 IL8 S100A9 MAGEA12 MAGEA12 SPRR1B FOXD1 POSTN KRT14 POSTN IGFBP3 KRT16 IGFBP3 TNC KRT6B MAGEA1 GJC1 POSTN FAM198B CTGF FAM198B SULF1 LILRB1 SERPINB3 IL13RA2 HLA-DRA CXCL10 CTGF IL13RA2 SERPINB4 GDF15 FAM198B FOXD1 HAS2 SERPINA3 LPPR4 CDH2 EIF2S3 CXCL9 PLAT ID3 CSTA SNX10 CD74 IL8 CCL20 MAGEA1 KRT5 SDC3 SULF1 MAGEA12 IL1B CCL5 doi:10.1371/journal.pone.0142443.t002 calciumbindingandcellularimmuneresponses.Apartfromthegenesmentionedabove,sev- eralothergenesthatareknowntobeexpressedinmelanomaarealsorepresentedinthetop-20 genelist,whichincludesGDF-15,CDH2,CTGF,SERPINA3andPOSTN.Interestingly,POSTN (Periostin)isreportedtoplayaroleinacceleratingmelanomametastasisviatheintegrin/mito- gen-activatedproteinkinase(MAPK)signallingpathway[23]. WealsoinvestigatedthelikelyglobalchangesduetoDEgenes,byperforminggeneontology (GO)enrichmentanalysis(usingathresholdp<0.01).Afterthefirsttransition,themost over-representedGOtermcorrespondstoskindevelopmentandimmuneresponse(S3Fig). However,asmelanomaprogresses,therepresentationofGOtermsassociatedwithbloodvessel developmentandadhesionincreases,whereastherepresentationofskindevelopmentrelated GOtermsdecreases.InLN,themostover-representedGOtermsforup-regulatedgenesare immuneresponserelatedprocessesandcellcycle. Inanutshell,theDEanalysisshowsexpectedoutcomesandfromthelargesetofDEgenes onlyasmallnumberofgenesareconserved,howevertheconservedgenesetcontainsseveral knownmelanomaandcancerassociatedgeneswhichvalidatestheauthenticityofourdatapro- cessingprocedureandqualityoftheprocesseddata. Co-expressionnetworks Forfurtheranalysisandpredictionofnoveldysregulatedgenes,fourgeneinteractionmaps wereconstructedusinggenepair-wisePearsonCorrelationCoefficient(PCC;seematerialsand methods).ThereasonforchoosingPCCisthatitisthemostwidelyusedmethodtoidentify thelinearrelationshipbetweentworandomvariablesandcapturedsignificantrelationships evenwhenthegeneshavenon-normaldistribution[24].Weidentifiedthecorrelationcoeffi- cient(r)distributionandalsoestimatedthep-valueandFalseDiscoveryRate(FDR)associated witheachcorrelationcoefficientwhichissummarizedinS4Fig.Moreover,theanalysisoffree PLOSONE|DOI:10.1371/journal.pone.0142443 November11,2015 6/30 GeneNetworkRewiringtoStudyMelanomaStageProgression Table3. Fundamentalnetworkproperties. Normal CnM CM LN Samplesize 97 183 246 116 Edges 4,269,800 368,136 1,605,642 1,835,739 Unconnectednodes 488 2406 2546 1228 Clusteringcoefficient 0.57 0.39 0.53 0.53 Density 0.04 0.004 0.15 0.02 Centralization(degree) 0.166 0.06 0.12 0.12 doi:10.1371/journal.pone.0142443.t003 nodecountatdifferentnetworkthresholdshelpedusindecidingthereferencenetworkthresh- old,i.e.r>=0.5andFDR<=0.05(seeS5Fig). Next,weevaluatedthenetworkqualitiesbycomparingitstopologicalpropertieswithran- domnetworksofthesameaveragedegreeandgenes.Randomnetworksweregeneratedusing Erdos-Renyimodel[25],whichcreateseverypossibleedgewithequalprobability'p',indepen- dentofotheredges.ThedegreeofavertexinsucharandomgraphfollowsaPoissondistribu- tionwithmostnodeshavingapproximatelysamenumberoflinks.Kolmogorov-Smirnovtest [26],whichevaluatesthesignificanceofdeviation,suggestsignificantdeviationoftheinferred genenetworkfromrandomnetworks(P<2.2e-16,FigureAinS6Fig). Generallyinbiologicalsystems,afewgeneshaveasignificantlyhighernumberofconnec- tionsascomparedtotheothergenesinthenetworki.e.theyactashubgenes.Suchscalefree topologysignifiesthenetworkrobustnessanditsresistancetorandomgeneknockout.Any scalefreenetworkisgenerallycharacterizedby“heavytailed”natureofnodedegreedistribu- tion(FigureBinS6Fig)[27].Weobserveasimilardistributionofdecreasingdegreedistribu- tionwithincreasinglinks,markingthepresenceofgeneswithfewconnectionstothe'hub' genes.Wethenperformedastatisticalmeasuretocheckscalefreetopologybytestingpower lawdistributionfittoeachofthenetworkmodelusingbootstrappinghypothesistest[28].We foundthattheNormalandCMnetworksdonotsignificantlyfittoapowerlawmodel (P<0.1),andhencethenetworksarenotscale-free.However,“heavytailed”degreedistribution oftheentirenetworkandtheirdeviationfromtherandomnetworksprovidesthecharacteristic ofatruecomplexbiologicalnetwork.Otherfundamentalnetworktopologicalparametersare givenintheTable3.Itisobservedthatdatasizehasnosignificantimpactonfundamental topologicalnetworkpropertieswithnormalbeingthemostconnectednetworkmodeland CnM,theleastconnected.Althoughasignificantvariabilityamongthefournetworkmodels exists,yetmetastaticsamplesarecloselyrelatedtoeachotherowningtosimilartopological propertieslikeedges,clusteringcoefficientanddegreecentralization.Overall,thesetopological propertiessuggestthedifferencesintheexpressionorconnectivityprofilesofthesamegenes undervariableconditions. Differentialnetworkanalysis Mostofthedifferentialnetworkanalysis,notspecifictomelanoma,reliesonthegenedegree difference(DifferentialConnectivity;DC)acrosstwoormorenetworkstoelucidatealtered geneset.InordertostudyDCinmelanomagenes,wemeasuredeachgeneDCacrossdifferent stagesofmelanomawithrespecttothenormalnetworkanditsstatisticalsignificanceinterms ofp-valueandfalsediscoveryrate(adj.porFDRadjustedp<0.01).Asexpected,asignificantly highdegreelossisobservedindifferentmelanomastages(Fig3A).Foreachmelanomastage, morethanatwo-foldincreaseinnumberofgeneswithconnectivitylossthangainexists.Ana- lysingthedifferentialconnectivityasafunctionofp-valuesuggeststhatamajorityofthegenes PLOSONE|DOI:10.1371/journal.pone.0142443 November11,2015 7/30 GeneNetworkRewiringtoStudyMelanomaStageProgression Fig3.A.Numberofgenessignificantlydifferentiallyconnected(DC;eitherlossorgain)ineachtumorstagewithrespecttothenormalstage network.B.ForeachgeneDCisplottedwithitsstatisticalsignificance(FDRadj.p-value)andwithdiseaseprogression,numberofgeneswith increasednetworkconnectivityincreasesunderthesignificantp-value(p<0.01)asrevealedbyhorizontallineonthepositivex-axis.C.Foreach gene,thecalculatedgainedco-expression(Rw)isplottedwithitsoverallDCscore.Thegraphexplainsthatwithdiseaseprogression,genestends PLOSONE|DOI:10.1371/journal.pone.0142443 November11,2015 8/30 GeneNetworkRewiringtoStudyMelanomaStageProgression tointensifythealterationsintheconnectiontypesbygainingaswellasloosingconnectionswhichmayormaynotalteritsDC.Fourquadrants suggestthatgenesmaygainnoveledgeswithincreasedDC(quadrantI)orgainnovelconnectiondespitelosingoverallDC(quadrantII)orlose edgeswithdecreasedDC(quadrantIII)orshowincreasedDCdespitelosingconnections(quadrantIV).Inthegraphsinglegenewasplottedtwo timestocorrelateRw(+ve)versusDCandRw(-ve)versusDCasbothofthetermswerefoundtobeindependent. doi:10.1371/journal.pone.0142443.g003 withhighdifferentialconnectivityhavesignificantlylowp-value(Fig3B).Thehigherlossof connectionswithlowp-valuestatisticallyconfirmstheedgelossfromnormaltocancertissue. Also,theasymmetricalnatureofthegraphssuggeststhepresenceofgeneswithahighernum- beroflostconnectionsthangained.Boththemetastaticstagesarecharacterizedbyanearly symmetricalgraph,whereasthenon-metastaticstageisasymmetric. Differentialconnectivityvsedgesrewiring Wedecipheredthealteredconnectiontypesineachofthemelanomastagesascomparedto thatinthenormal,byoverlappingtherespectivenetworks.Inordertoscrutinizetheassump- tionwhetheredgerewiringproducesdifferentresultsthanDC,wecomparedtheresultspre- sentedbythetwoapproaches. Foreachofthemelanomastages,wefoundahighnumberofgeneswithnovelconnections in-spiteoftheoveralldegreelossascomparedtothatinthecontrol(Fig3C;IIquadrant).Each comparisonischaracterisedbyadiagonalrelationshipbetweenconnectivitylossandnodes lost(EL);connectivitygainandnovelnodesgain(EG).Thedensenatureofthegraphreflects highproportionsofedgerewiringduringmetastasis,LNthemostgene-rewiredstate.Thus, melanocytetomelanomatransformationbeginswithsignificantlossofnetworkrobustness, however,progressiontometastatictransformationoccursonlyaftermulti-edgegains.We believethatthemulti-edgegain,displayedbyagenesetcouldbeanimportantphenomenon fortransformedbiologicalactivityandcanhelpidentificationofgenesinthemostaberrantly rewirednetworks.Therefore,intheupcomingpartofthismanuscript,wewillprimarilyfocus ontheshortlistedDCgeneswithsignificantlygainednovelconnections. Multi-edgerewiringinmelanomaprogression Forelucidationofaberrantlynetworkedgenes,networkconnectivityrewiringanalysiswasper- formedforastage-wiseassessmentoftheedgerewiring.Theanalysisrevealsthatmelanocyte tomelanomaconversioninitiateswithasignificantedgesloss(highestEL=2059)ascompared tothenumberofedgesgained(highestEG=422).Weobserveasuddenincreaseingenesdis- playingEG,whilemaintainingitsconnectivitylossafterthesecondtransition(Fig4A).Box- plotanalysisclearlyillustratestheeffectofmelanomaongeneconnections,thenumberof edgeslostisfoundtobeashighas~2000inallthestagesofmelanoma(Fig4B).Inorderto confirmifthegeneswithlostconnectionsinCnMretainitsconnectivityprofileinlaterstages, wecomparedthegenesthatlostmorethan20edges(arbitrarilyselectedvalue)acrossdifferent stagesofmelanoma.Morethan90%oftheCnMgenesthatdisplayedELcontinuetoloose edgesinCMandLN(Fig4C).Thus,globalanalysissuggeststhatalargesetofgenesactively rewiresitsconnectiontypes,manyofwhichgainnovelconnectionsand,ifnotsignificantly down-regulated,mayactivelybeinvolvedinaberrantnetworksinmelanomaprogression, especiallyinthecellsprogrammedformetastasis. Acloseranalysisoftheresultsrevealthatthegenethatgainedmaximumedgesduringthe firsttransitionisKIDINS220(EG=422,EL=211)whilethegenewiththehighestlossofcon- nectionsisC1orf216(EG=4,EL=2059).Duringthenon-metastatictometastatictransition, thegeneLRP5gainedthehighestnumberofconnections(EG=1283,EL=113),while C1orf216continuedtobethegenewiththehighestedgeloss(EG=0,EL=2065).Afterthe PLOSONE|DOI:10.1371/journal.pone.0142443 November11,2015 9/30 GeneNetworkRewiringtoStudyMelanomaStageProgression Fig4.A.Linegraphtoobservethedifferencebetweenedgeslost(color:green)andgain(color:red)byasinglegeneindifferentmelanoma stages.Theobservedtrendsuggestsconsistentandsignificantincreaseinnoveledgesgain,aftertheonsetofmetastasis.B.Boxplotshowing thedifferencesinthenumberofnoveledgesgained(color:red)andlost(color:green)bygenesineachstage.C.Comparisonforthegenesthat displayedconnectivityloss(edgeloss>20)acrossdifferentmelanomastages.Venndiagramsuggestsremarkablesimilarityofgenebehavior acrossmutuallyexclusivedatasetsandlossofconnectionsisaconservedphenomenon. doi:10.1371/journal.pone.0142443.g004 thirdtransition,geneASNSD1gainedthehighestnumberofedges(EG=1254,EL=32) whereasgeneARMC9(EG=19,EL=1924)lostmaximumnumberofedges.However,the analysisoftopfiverewiredgenesdonotlistanyknownrecognizedmelanomamarker,thus suchgenesmaybeexploredaspotentialtargets. Diseasegenesaresignificantlyrewired Bycomparingthedatasetofeachstagewiththenormaldataset,wegeneratedagenelistofdis- easegenesthatareeithersignificantlyupregulatedorabnormallyconnected(seematerialsand methods).Ouraimwastoderiveanovelgenesetthatcouldpossiblyaccountforstage-wise networkaberrationsduringmelanomaprogression.Interestingly,themetastaticstagesinclude averyhighnumberofshortlistedgenes(i.e.diseasegenes).Thehighestnumberofputativedis- easegenesisfoundinLNascomparedtothatinotherstages(Fig5A–5E).Itisevidentthatthe highnumberofdiseasegenesmayhelpustounderstandthemetastaticbehaviourofthegenes inthedatasetandidentifynovelassociatedbiomarkers.Itisfoundthat60%oftheCnMgenes retainitsactivityprofileinoneoftheothertwometastaticstages.Wefoundthat59%ofthe uniquelystagespecificdiseasegenesinLNthatmaximallyshares1317geneswithCM. PLOSONE|DOI:10.1371/journal.pone.0142443 November11,2015 10/30