Underexpression of Specific Interferon Genes Is Associated with Poor Prognosis of Melanoma

15 Pages·2017·1.83 MB·English

RESEARCHARTICLE
Underexpression of Specific Interferon Genes
Is Associated with Poor Prognosis of
Melanoma
AamirZainulabadeen1,2☯,PhilipYao1,3☯,HabilZare1*
1DepartmentofComputerScience,TexasStateUniversity,SanMarcos,Texas,UnitedStatesofAmerica,
2DepartmentofComputerScience,PrincetonUniversity,Princeton,NewJersey,UnitedStatesofAmerica,
3ElectricalEngineeringandComputerScience,UniversityofMichigan,AnnArbor,Michigan,UnitedStatesof
America
a1111111111
a1111111111 ☯Theseauthorscontributedequallytothiswork.
a1111111111 *[email protected]
a1111111111
a1111111111
Abstract
Becausetheprognosisofmelanomaischallengingandinaccuratewhenusingcurrentclini-
calapproaches,cliniciansareseekingmoreaccuratemolecularmarkerstoimproverisk
OPENACCESS
models.Accordingly,weperformedasurvivalanalysison404samplesfromTheCancer
Citation:ZainulabadeenA,YaoP,ZareH(2017)
GenomeAtlas(TCGA)cohortofskincutaneousmelanoma.Usingourrecentlydeveloped
UnderexpressionofSpecificInterferonGenesIs
AssociatedwithPoorPrognosisofMelanoma. genenetworkmodel,weidentifiedbiologicalsignaturesthatconfidentlypredicttheprogno-
PLoSONE12(1):e0170025.doi:10.1371/journal. sisofmelanoma(p-value<10−5).Ourmodelpredicted38casesaslow–riskand54cases
pone.0170025
ashigh–risk.Theprobabilityofsurvivingatleast5yearswas64%forlow–riskand14%for
Editor:RogerChammas,UniversidadedeSao high–riskcases.Inparticular,wefoundthattheoverexpressionofspecificgenesinthe
Paulo,BRAZIL
mitoticcellcyclepathwayandtheunderexpressionofspecificgenesintheinterferonpath-
Received:September26,2016 wayarebothassociatedwithpoorprognosis.Weshowthatourpredictivemodelassesses
Accepted:December26,2016 theriskmoreaccuratelythanthetraditionalClarkstagingmethod.Therefore,ourmodelcan
helpcliniciansdesigntreatmentstrategiesmoreeffectively.Furthermore,ourfindingsshed
Published:January23,2017
lightonthebiologyofmelanomaanditsprognosis.Thisisthefirstinvivostudythatdemon-
Copyright:©2017Zainulabadeenetal.Thisisan
stratestheassociationbetweentheinterferonpathwayandtheprognosisofmelanoma.
openaccessarticledistributedunderthetermsof
theCreativeCommonsAttributionLicense,which
permitsunrestricteduse,distribution,and
reproductioninanymedium,providedtheoriginal
authorandsourcearecredited.
DataAvailabilityStatement:Allrelevantdataare
Introduction
withinthepaperanditsSupportingInformation
files.Alternatively,datacanalsobeaccessedatthe Cutaneousmelanomaisamalignancyofmelanocytes.Itisthemostcommontypeofskincan-
followinglink:https://tcga-data.nci.nih.gov/docs/ cer.TheAmericanCancerSocietyestimatesthatover73,000newcaseswerediagnosedin
publications/skcm_2015/Alternatively,datacanbe
2015intheUnitedStatesandabout10,000deathsarecausedbymelanomaeachyear[1].The
accessfromTCGADataPortalorfromthe
prognosisofmelanomaishighlyvariable[2].Forinstance,the5–yearoverallsurvivalratecan
followingpage(doi:10.1016/j.cell.2015.05.044):
beashighas97%forstageIandaslowas3%forstageIV[3,4].Almostallcommontreatment
https://tcga-data.nci.nih.gov/docs/publications/
skcm_2015/. optionsformelanoma,includingsurgery,chemotherapy,andradiationtherapy,haveharmful
andseveresideeffects.Therefore,itiscriticaltoidentifypatientswhoarenotatasignificant
Funding:ThisstudywassupportedbytheNational
riskofmetastasisanddeathduetothedisease.Thepredictivepowerofclinicalfactorsislim-
ScienceFoundation(NSF)throughtheResearch
ExperiencesforUndergraduates(REU)program ited[3,5,6](e.g.,stagingbasedonthetumorsizeandthenumberofmetastaticsentinel
PLOSONE|DOI:10.1371/journal.pone.0170025 January23,2017 1/15
IFNPathwayandMelanomaPrognosis
(CNS1358939),andalsothroughtheinfrastructure lymphnodes[7]),thereforecliniciansareseekingmoreaccuratemolecularmarkersto
grant(CRI1305302). improveriskmodelsandtoavoidunnecessarytreatmentoflow-riskpatients[8–10].
CompetingInterests:Theauthorshavedeclared Geneexpressionprofilesignatureshaveusefulinformationonthemolecularstatusofcells
thatnocompetinginterestsexist. andtheycanpredicttheprognosisofmanycancers[11–14],includingmelanoma[10,15–18].
Forexample,Onkenetal.discoveredageneexpressionprognosticsignaturethatsignificantly
improvedtheclassificationofuvealmelanomacomparedtotraditionalstaging[15].Thatis,
theyshowedthatthederegulationoftheLEDA,FZD6,andENPP2genespredictsmetastatic
death(p-value<10−4).Inthefollow–upstudies,theyextendedtheirtesttoinclude15genes
[19]andtheirextendedtestcorrectlyclassified446(97%)ofthe459studiedcasesintolow–
risk(i.e.,atleast95%chanceof5–yearmetastasis–freesurvival)andhigh–risk(i.e.,notmore
than20%chanceof5–yearmetastasis–freesurvival)groups[20].
Recently,Geramietal.performedameta-analysisonseveralpublishedgenomicanalysesof
cutaneousmelanomatumors[15,21–27].Basedonthegeneontology[28]ofthefrequently
reportedgenes,theyidentified28discriminantgenesincludingBAP1,MGP, SPP1,CXCL14,
CLCA2,S100A8,BTG1,SAP130,ARG1,KRT6B,GJA1,ID2,EIF1B,S100A9,CRABP2,KRT14,
ROBO1,RBM23,TACSTD2,DSC1,SPRR1B,TRIM29,AQP3,TYRP1,PPL,LTA4H, andCST6
[4].Theyusedtheexpressionofthesegenestotrainageneralizedlinearmixedmodelwitha
radialbasisfunctionkernel.Theresultingsignatureclassified268primarycutaneousmela-
nomatumorsintolow-riskandhigh-riskgroups,with5–yeardisease–freesurvivalratesof
97%and31%,respectively[4].Thistestiscommerciallyavailableasadiagnostictoolcalled
DecisionDx-Melanoma[29,30],butitisnotyetrecommendedbytheNationalComprehen-
siveCancerNetwork[10].Itsbenefitstothepatientsmustbeconfirmedusingprospective
clinicaltrialsthatincludesignificantlylargercohortswithmorerepresentativemetastaticchar-
acteristics[18,31].
Wehypothesizedthatthereisroomforsignificantimprovementinthemelanomaprognos-
ticteststhroughrigid,unbiased,andcomprehensiveanalysisofgeneexpressionprofiles[10,
32].Accordingly,weappliedarobustlarge-scalenetworkanalysistothegeneexpressiondata
of404samplesfromTheCancerGenomeAtlas(TCGA)cohortofskincutaneousmelanoma
[33].Theaimofourstudywastoidentifythemolecularsignaturesthatpredicttheprognosis
ofmelanomaandtosegregatepatientsintolow–,medium–,andhigh–riskgroups.Our
approachisbasedonco-expressionnetworkanalysis,andweuseeigengenesasinformative
prognosticssignatures[34].
MaterialsandMethods
TheTCGADataset
SeveralmRNAexpressionprofilingdatasetshavebeenproducedtostudymelanomaprognosis
[10,23,27,33,35–40].WeusedtheTCGA2STATpackagetodownloadgeneexpressiondata
fromTheCancerGenomeAtlas(TCGA)repository[41].Specifically,wedownloadedRNA-
Seqdatafromtheskincutaneousmelanomacohortof473patients[33]andusedRPKMval-
ues(i.e.,readsperkilobaseoftranscriptpermillionmappedreads)[42]asameasureofgene
expression.Wemanuallydownloadedthecorrespondingclinicaldata,including:a)informa-
tionregardingthelaststatusofeachcase(i.e.whethertheeventofdiseaserecurrenceorpro-
gressionoccurred),b)thelengthofthedisease–freetimeperiod(i.e.,thetimefromtheinitial
melanomadiagnosisuntilthiseventoruntilthelastfollow–updateiftheeventdidnotoccur),
andc)theClarkscalestageofthemelanoma,whichwasdeterminedusingclinicopathological
featuressuchasthesize,number,andlocationofmetastases[43].
WecomputedtheSpearman’srankcorrelationbetweenthedisease–freetimeandgene
expression[44].Spearman’srankcorrelationismorerobustthanPearsoncorrelationanditis
PLOSONE|DOI:10.1371/journal.pone.0170025 January23,2017 2/15
IFNPathwayandMelanomaPrognosis
therecommendedapproachforskeweddistributions[45].Forinstance,Mukakashowedthat
afterremovaloftheoutliers,thechangeintheSpearman’scorrelationcanbenegligibleunlike
thePearsoncorrelation[46].Consistentwiththeapproachtakenbyotherscholars[47],we
usedonlythetopthirdofgenes(6,834)thatweremostcorrelatedwiththedisease–freetimein
ournetworkanalysis.Weconsideredanyprogressedorrecurredtumorashighrisk(n=263).
Intherestofthetumors,whichwerenotreportedtorecurduringthefollow–upperiod,any
casethathadatleastfiveyearsoffollow–updatawasconsideredlow-risk(n=33)[48–51].All
resultspresentedherecanbeconvenientlyreproducedusingoursupplementarycode(S4
File).
Validationdatasets
ToconfirmthefindingsthatweobtainedusingtheTCGAdataset,wevalidatedthemusing
twoindependentdatasets.Specifically,weusedtheLeedsmelanomageneexpressionset1
[39],whichispubliclyavailablethroughEuropeanMolecularBiologyLaboratory–European
BioinformaticsInstitute(EMBL–EBI,accessionnumber:E-MTAB-4725).Forbrevity,werefer
tothisdatasetasLEEDinthispaper.Thiscohortcompriseswhole–genomemRNAexpression
of204primarymelanomatumors,whicharemeasuredusingIlluminaDASLHT12.4.Kolesni-
kovetal.normalizedthegeneexpressionvalueswithquantilemethodafterbackgroundcor-
rection.WemanuallydownloadedthegeneexpressionandclinicaldatafromtheEMBI–EBI
ArrayExpressdatabase(http://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-4725/)[52].
Weusedthe“lastfollowup”(time)and“viability”(deathevent)columnsfromtheclinical
data.
WealsousedasimilarcohortthatwasproducedattheLundUniversityin2015[40].This
datasetispubliclyavailablethroughEMBL–EBI(accessionnumber:E-GEOD-65904)andalso
throughGeneExpressionOmnibus(GEO)(accessionnumber:GSE65904)[53].Forbrevity,
werefertothisdatasetasLUNDinthispaper.Cirenajwisetal.extractedtotalRNAfrom214
fresh–frozenmelanomatumorsandperformedgenome-wideexpressionprofilingusingIllu-
minaHumanHT-12V4.0BeadChiparrays.Wedownloadedthegeneexpressiondatausing
GEOquerypackage(Version2.40.0)[54].Wedownloadedthecorrespondingclinicaldata
fromtheEMBI–EBIdatabase(http://www.ebi.ac.uk/arrayexpress/experiments/E-GEOD-
65904/)andusedthe“diseasespecificsurvival”(time)and“diseasespecificdeath”(event)
columns.
Genenetworkanalysis
Weappliedtheweightedgeneco-expressionnetworkanalysis(WGCNA)(Version1.51)toall
473availablesamplestobuildagenenetworkandtoclusterthegenesintogenemodules(clus-
ters)[55].Specifically,weusedthecalculate.betafunctionwiththedefaultparameters
toinferthatthesoft-thresholdingpowerfornetworkconstructionwas7.Weidentified13
genemodulesusingthewgcna.one.stepfunctionfromthePigengenepackage[34],
whichisawrapperfortheblockwiseModulesfunction,withpower=7andleftthe
remainingargumentsasdefaults.WGCNAcouldnotconfidentlyassign1,404genestoanyof
themodules,becausethesegeneshadlittlecorrelationwiththeothergenes.Wecallthesetof
theseoutliergenesModule0.
Computingeigengenes
Aneigengeneofamoduleisaweightedaverageoftheexpressionofallthegenesinthatmod-
ule.Theseweightsareadjustedsothatthelossofbiologicalinformationisminimized[56,57].
Weusedprincipalcomponentanalysis(PCA)tocomputeeigengenes.First,webalancedthe
PLOSONE|DOI:10.1371/journal.pone.0170025 January23,2017 3/15
IFNPathwayandMelanomaPrognosis
numberofhigh–riskandlow–riskcasesusingoversamplingsothatbothgroupshadcompara-
blerepresentativesintheanalysis.Specifically,werepeatedthedataofeachhigh–riskand
low–riskcase6timesand45times,respectively.Thisapproachprovideduswith1,485and
1,578samplesfromeachgroup,respectively.Oversamplingwasnecessaryforcomputing
theeigengeneofamodule.Becauseaneigengeneisthefirstprincipalcomponentofthe
module,itwouldbebiasedtowardsthehigh–riskgroup,whichhasaroundeighttimesmore
samplesthanthelow–riskgroup.Oversamplingresolvesthisissue.Then,weappliedthe
moduleEigengenes()functionfromtheWGCNApackagetotheoversampleddata.This
functioncomputedthefirstprincipalcomponentofeachmodule,whichmaximizedthe
explainedvariance,thusensuringaminimumlossofbiologicalinformation.Weusedthe
project.eigenfunctionfromthePigengenepackage(Version0.99.23)toinferthevalues
ofeigengenesforallofthe473samplesintheTCGA,the204samplesintheLEEDS,andthe
214samplesintheLUNDSdatasets(S1File)[34].
Survivalanalysis
Weusedthe14inferredeigengenesascovariates(prognosticfeatures),andweincludedonly
the404samplesforwhichthefinalstatusandthesurvivaltimewereavailable.Weusedthe
glmnet()functionfromtheglmnetpackage(Version2.0-5)[58]toperformapenalized
Coxregressionanalysis[59,60].Wesetα=1tousetheleastabsoluteshrinkageandselection
operator(Lasso)[61].TheLasso,alsoknownasL regularization,enforcesmostofthecoeffi-
1
cientsofthecovariates(eigengenes)intheCoxproportionalhazardsmodeltobezero.
Thereby,itidentifiesthemodulesthatarethemostassociatedwithsurvival.
Toevaluatethesignificanceoftheselectedmodulesinpredictingthesurvivaltime,we
fittedanacceleratedfailuretime(AFT)modeltotheselectedeigengenes[62].Weusedthe
survregfunctionfromthesurvivalpackage(Version2.39-4)[63],settheWeibulldistribu-
tionwithscale=1asthebaselinehazardfunction,andusedthedefaultvaluesfortherestof
theparameters.Weusedthefittedacceleratedfailuretimemodeltopredictthesurvivaltime
ofeachsample.Wechosetwothresholdsforthepredictedvaluesthatmaximizedtheprecision
oflow–andhigh–riskpredictions.Thesamplesthathadapredictedsurvivaltimebetweenthe
twothresholdswereconsideredmedium–risk.Weusedthesurvfitfunctiontoobtaina
Kaplan-Meiersurvivalcurveforeachoftheriskgroups[64].Weusedthesurvdifffunction
totestwhetherthesurvivalcurvesthatcorrespondtohigh–riskandlow–riskgroupsdiffersig-
nificantly.Thisfunctioncomputedthelog-rankp-valueofthecorrespondingMantel-Haens-
zeltest[65].
Cross–validation
Weperformed5–foldcross–validationtoconfirmthattheselectionofthemodulesbythe
penalizedCoxregressionisrobustwithrespecttochoosingthesamples.Weusedallofthe14
eigengenescorrespondingtothemodulesthatwereidentifiedbyourgenenetworkanalysis.
Wedidnotrecomputethemodulesoreigengenes,instead,werepeatedthepenalizedCox
regressionmodelasfollows.Therewere404samplesforwhichthefinalstatusandthesurvival
timewereavailable.Werandomlydividedthese404samplesintofivedivisions.Thesedivi-
sionshadanalmostequalsizeandequalnumberofhigh–risksamples.Wesetasideonedivi-
sionandperformedapenalizedCoxregressionanalysisontherestofthesamplesusingall14
eigengenes.Werecordedtheselectedmodulesandrepeatedthisprocedurefivetimes.Weran
thisexperiment10timeswithdifferentseedsandcountedthefrequencyoftheselectedmod-
ulesineachrun(S2File).Onaverage,thethreemostfrequentmoduleswereselected4.9,4.8,
PLOSONE|DOI:10.1371/journal.pone.0170025 January23,2017 4/15
IFNPathwayandMelanomaPrognosis
and4.2times,respectively.Incontrast,allothermoduleswereselected0.6timesorless,on
average.
Results
Usingcoexpressionnetworkanalysis,weidentified13modulesofhighlycoexpressedgenes.
Thesizeoftheesmodulesrangesfrom48to2,247geneswithameanof418,amedianof88,
andastandarddeviationof629(S1Fig).Wecomputedaneigengeneforeachmodule,which
summarizesthebiologicalinformationofthemoduleintoonevaluepersample.Weused
theseeigengenesasbiologicalsignatures(features)toperformasurvivalanalysis.
ThepenalizedCoxregressionconsistentlyselectedthreemodulesasthemostassociated
moduleswithdisease–freesurvival(Methods).Thesethreegenemodulesincludetheoutlier
module(with1,404genes,ME0),theninthlargestmodule(with58genes,ME9),andthe
twelfthlargestmodule(with52genes,ME12)(S3File).Theoutliermoduleconsistsofgenes
thathavetoosmallacorrelationwithothergenestobeincludedinanyofthemodules.Hyper-
geometrictestsrevealedthattheothertwoselectedmodulesareassociatedwiththemitoticcell
cycleandtheinterferon(IFN)pathway,respectively[66](S2Fig).Thatis,15(29%)ofthe58
genesinthelargermodulearemembersoftheReactomemitoticcellcyclepathway,whichhas
454genes(adjustedp-value<10−9)[67].ThesegenesincludeAURKB,CCNE1,CDCA8,
CDK4,CENPO,GINS2,H2AFZ,LIG1,PKMYT1,PLK1,PTTG1,SKA1,TUBA1B,TUBA1C,
andTYMS.AccordingtotheCoxmodel,overexpressionofthesegenesisassociatedwithpoor
prognosis,whichisexpected[68,69].
Thesmallermodule(ME12)has52genes.Interestingly,16(31%)ofthesegenesaremem-
bersoftheReactomeinterferonsignalingpathway,whichhas158genes(adjustedp-
value<10−21).ThegenesintheoverlapincludeDDX58,EIF2AK2,GBP3,HERC5,IFIT1,
IFIT2,IFIT3,OAS1,OAS2,OAS3,PSMB8,SP100,STAT1, UBE2L6,USP18,andXAF1.Our
modelindicatesthattherelativelyhigherexpressionofthesegenesisassociatedwithagood
prognosisinmelanoma.OurinsilicooverrepresentationanalysisshowedthattypeIinterferon
signalingpathway(GeneOntologyaccessionnumber:60337[28]),whichincludes63genes,is
thebiologicalprocessesthathasasignificantoverlapwiththisgeneset.Specifically,theoverlap
consistsof12genes,whichis75timesmorethanexpected(p-value=10−15)[70].Thismodule
isverystablewithrespecttoselectionofsamples.Toconfirmthis,wereconstructedthecoex-
pressionnetwork10timesusingonly426(90%)randomlyselectedsamples.Alloftheresult-
ingnetworkshadasimilarmodule,i.e.,meanandmedianoftheJaccard[71](Tanimoto[72])
similaritywere0.92and0.93,respectively.
Tofurthervalidatetheassociationofthesethreeselectedmoduleswithmelanomaprogno-
sis,wefittedanacceleratedfailuretime(AFT)modeltothecorrespondingeigengenes,and
classifiedthepatientsintolow–,medium–,andhigh–riskgroups(Methods)[62].Wecom-
paredtheKaplan-Meier(KM)curvesofthesegroups[64](Fig1).OurAFTmodelpredicted
that38caseswerelow–risk.Thesecaseshadasignificantlyhighersurvivalratethanthe54
casesthatwerepredictedtobehigh–risk(log-rankp-value<8×10−6[73]).Forinstance,the
probabilityofsurvivingforatleastfiveyearswas0.64forlow–riskand0.14forhigh–riskcases.
Excludingtheinterferonmodulehasanegativeimpactonthestatisticalpowerandthep-value
(Fig1).Specifically,whilethenumberofcasesinthepredictedlow–riskgroupdoesnot
increasedramatically(i.e.,twocases,only5%improvement),thenumberofpredictedhigh–
riskcasesdecreasesconsiderablyto35(i.e.,35%decline).Thatis,usingtheinterferonpathway
module,themodelcanidentifymorehigh–riskcaseswithoutsacrificingtheaccuracy.
Asexpected,thecasesthatwerepredictedtobehigh–riskhadgenerallymoreadvanceddis-
easeaccordingtothetraditionalClarkscalestageofmelanoma(Table1).Inparticular,the
PLOSONE|DOI:10.1371/journal.pone.0170025 January23,2017 5/15
IFNPathwayandMelanomaPrognosis
Fig1.Kaplan–Meiersurvivalcurves.Thep-valuesindicatethatthedifferencebetweenthelow–riskgroup(green)andthehigh–riskgroup(red)is
statisticallysignificant.Usingallthethreemodules,whichareassociatedwiththeinterferonpathway,mitoticcellcycle,andoutliers;resultsinabetterp-
value(a)comparedtoamodelwithouttheinterferonpathway(b).Theorangehorizontallinesindicatethatbothmodelshavesimilaraccuracies.However,
includingtheinterferonpathwayimprovesthep-value,becausemoresamplesareclassifiedintotal(i.e.,38low–riskplus54high–riskcasesin(a),
comparedto40low–riskplus35high–riskcasesin(b)).
doi:10.1371/journal.pone.0170025.g001
majorityofthehigh–riskcaseswereinstageIVorV(36cases,66%).Incontrast,only13
(35%)ofthelow–riskcaseswereinstageIVorV.Interestingly,6(46%)ofthese13caseswere
disease–freeformorethaneightyears,whichindicatesthat,forthesecases,theClarkscale
stagingislessaccuratethanourpredictions.Also,7caseshadstageI,II,orIIImelanoma,but
theywereclassifiedashigh–riskbyourmodel.Only2(29%)ofthesecasessurvivedmorethan
2years.
Tofurthervalidatetheassociationbetweentheinterferonmoduleandmelanomaprogno-
sis,weusedtheLEEDSandLUNDgeneexpressiondatasets,whichinclude204and214
melanomasamples,respectively.WeinferredtheeigengenesandfittedanAFTmodeltothe
eigengenestoclassifythesamplesintolow–,medium–,andhigh–riskgroupsineachdataset
Table1.ThedistributionofmelanomaClarkstagesineachriskgroup.Also,thepercentageofeachstageclassineachriskgroupisshown.Thelow–
riskgroupisenrichedinpatientsatstagesIIIandIV.Thehigh–riskgroupisenrichedinpatientsatstagesIVandV.
Stage I II III IV V Unknown Total
(%) (%) (%) (%) (%) (%)
Low–risk 0 1 11 12 1 13 38
(0) (3) (29) (32) (3) (34)
Medium–risk 4 16 54 105 37 96 312
(1) (5) (17) (34) (12) (31)
High–risk 1 0 6 24 12 11 54
(2) (0) (11) (44) (22) (20)
Total 5 17 71 141 50 120 404
doi:10.1371/journal.pone.0170025.t001
PLOSONE|DOI:10.1371/journal.pone.0170025 January23,2017 6/15
IFNPathwayandMelanomaPrognosis
(Methods).SimilartotheanalysisonTCGAdataset,weusedthethreeeigengenescorrespond-
ingtotheoutlier,mitoticcellcycleandinterferonpathwaymodules.
OurAFTmodelpredicted96casesintheLUNDdatasettobelow–riskand45casestobe
high–risk.Thepredictedlow–riskgrouphadasignificantlyhighersurvivalratethanthehigh–
riskgroup(log-rankp-value2×10−3,S3aFig).Excludingtheinterferonpathwaymodulefrom
theAFTmodelresultedinreducingthenumberofpredictedlow–riskto49casesandpre-
dictedhigh–riskto25cases,andalso,alesssignificantp-value(3×10−2).Thatis,withoutthe
interferonpathway,thenumberofpredictedlow–andhigh–risksampleswouldbereducedby
almosthalfleaving67(33%)moresamplesinthemedium–riskgroup.Thisshowstheneces-
sityoftheinterferonpathwaymoduleinpredictingtheprognosis.TheresultsintheLEEDS
datasetfollowedasimilarpattern.Specifically,excludingtheinterferonpathwaymodulefrom
theAFTmodelresultedinreducingthenumberofpredictedlow–risksamplesfrom34cases
to31caseswhileslightlydecreasingthenumberofpredictedhigh–risksamplesfrom29to28.
Nevertheless,thelog-rankp-valueincreasedbyalmosttwoordersofmagnitude,from9×10−6
to7×10−4,indicatingthatthepredictionofsurvivalislessaccuratewithouttheinterferon
pathwaymodule(S3candS3dFig).
Discussion
Predictingtheprognosisofmelanomaisclinicallyusefulandimportant[10].Todate,mostof
thestudiesthataimatpredictingmelanomasurvivalbasedongeneexpressionhavebeenlim-
itedintheirnumberofgenes,numberofsamples,ortheirfollow–uptime[4,15–27,30,74].
Weperformedgenenetworkanalysison470melanomacasestoextendthepreviousstudies
andtoidentifynovelprognosticsignatures.
Thisisthefirststudytoshowthattheunderexpressionofspecificgenesfromtheinterferon
pathwayinmelanomatissuesisasignofpoorprognosis.Theroleoftheinterferonpathwayin
othercancerswerestudiedbyothers[75–80].Ingeneral,defectsininterferonsignalingresults
indysfunctionoftheimmunesystem[81].However,itsassociationwithmelanomawasprevi-
ouslyshownonlyinvitro[81–84].
Interestingly,ourinterferonmodulehas17genesincommonwiththe274genesthatHoek
etal.reportedtobedownregulatedinmelanomacelllines[82].Thisisasignificantoverlap(p-
valueofthehypergeometrictest<10−19).Similarly,thelistofthetop25genesthatCritchley
etal.reportedtobedifferentiallyexpressedinperipheralbloodmononuclearcell(PBMC)
samplesofmelanomapatientshasasignificantoverlapwithourinterferonmodule(13genes,
p-value<10−27).Comparedtoinvitroexperiments,ouranalysisprovidesmuchstrongerevi-
dencefortheroleoftheinterferonpathwayinmelanoma,becauseourstudyisbasedonthe
survivalanalysisofarelativelylargecohortofpatientswithanextendedfollow–uptime.The
totalnumberofpatientsclassifiedaslow–riskorhigh–riskincreasesfrom75to92(a23%
improvement)whenweincludetheinterferonmoduleinourpredictivemodel.Thisimprove-
ment,aswellasthedecreaseinp-value,indicatethattheinformationintheinterferonpath-
wayisessentialforpredictingtheprognosis.
However,functionalstudieswillbeneededtodeterminethemechanismandimpactofthe
interferonpathwayonmelanomaprognosis.Onechallengeindesigningsuchastudyispossi-
blytherelativelylownumberofsamplesthatcouldbeassociatedwiththegenesintheinter-
feronmodule.Forexample,inourstudyonthe404TCGAcases,includingthesegenesinthe
predictivemodelledtoclassificationofonly17(4%)morecases.Therefore,afollow-upfunc-
tionalstudymostlikelyneedstoinvestigateatleasthundredsofsamplesinordertoincludea
fewsamplesthatareassociatedwiththeinterferonpathway.
PLOSONE|DOI:10.1371/journal.pone.0170025 January23,2017 7/15
IFNPathwayandMelanomaPrognosis
Theprobabilityofsurvivingforatleastfiveyearsis0.64and0.14forourpredictedlow–risk
andhigh–riskgroups,respectively.Theaccuracyofourpredictivemodeliscomparablewith
previousstudiesonskincutaneousmelanomathatwerebasedongeneexpression[10].For
instance,Sivendranetal.developedalog-rankMantel–Coxtestbasedona21-geneexpression
signatureandulceration[74].Theirtestislessspecificthanoursinpredictinglong-term
poorprognosis.Thatis,theytested48patientsandreportedthattheprobabilityofsurviving
foratleastfiveyearswas60%and35%fortheirpredictedlow–riskandhigh–riskgroups,
respectively.
Geramietal.reportedthattheDecisionDx-Melanomatest[29,30]classified268primary
cutaneousmelanomatumorsintolow-riskandhigh-riskgroupswithatleast5–yeardisease–
freesurvivalratesof97%and31%,respectively[4].Theirlow–riskspecificityismorethan
ours,andtheirhigh–riskspecificityislessthanours.However,thereportedassessmentofthe
sensitivityandspecificityofthistestiscontroversial,becausetheircohortwasnotrepresenta-
tiveofthegeneralprimarymelanomapatientpopulation[10,18].Ourresultsarenotcompa-
rablewiththeresultsreportedbyOnkenetal.[20],becausetheystudieduvealmelanoma,
whichisdifferentfromskincutaneousmelanoma[85].
Onereasonunderlyingthedifficultyinassessingthesurvivalrateofmelanomaistherela-
tivelyhighheterogeneityofgeneticmutations,whichresultsinsubpopulations(clones)that
areresistanttotherapy[86–88].Recently,Tiroshetal.usedsingle-cellsequencingtoshowthat
distinctcloneswithinatumorcanhavedifferentgeneexpressionprofiles,andtherefore,can
interactwiththeirmicroenvironmentdifferently[89].Investigationofthegenesinourinter-
feronmoduleandthecorrespondingsignalingproteinsmayrevealhowtheinteractions
betweenmalignantcellsandtheirmicroenvironmentsismodulated.Forinstance,therela-
tivelyhighexpressionoftypeIinterferonsignalingpathwaygenesinlow–riskmelanomacases
canbeassociatedwithantitumorresponseoftheimmunesystem[90].Futureworkinthis
directioncanleverageavailabletechniquesfor1)clonaldecompositionbasedongeneticmuta-
tions[91–97],and2)deconvolutionofgeneexpressionprofilesintosignaturesthatarespecific
toacell-typeortumorclone[98–101].
Onelimitationofourpredictivemodelforclinicaluseistherelativelyhighnumberofcases
thatwereclassifiedasmedium–risk(312,77%).Thiscanbeaddressedbyimprovingthepre-
dictivemodelinafollow–upstudyor,alternatively,byusinganotherprognostictestforthe
medium–riskcases.Comparedtoothertests,ourpredictivemodelisbasedonarelatively
largenumberofgenes.Thisisadouble-edgedsword.Theinclusionofalargenumberofgenes
makesourmodelrobustwithrespecttorandomchangesintheexpressionofoneorseveral
genes,whicharecommonduetotechnicalorbiologicalnoise.Ontheotherhand,thelarge
numberofgenesmakesourtestdifficulttoapplyinclinicalsettings.Thiscanbeaddressedby
excludingthegenesthathavearelativelysmallercontributiontotheeigengenes.Specifically,a
greedyalgorithmcanbeusedtoexcludethegenesthathaveasmallerabsoluteweight(load-
ing)[102].Furtherfollow–upexperimentsondifferentdatasetswillbeneededtoshowthat
suchamodificationdoesnotaffecttheaccuracyofthemodelandtoensurethattoomany
geneswillnotbeexcluded.Nevertheless,follow-upstudiescanexaminethetherapeuticvalue
ofthegenesidentifiedinourstudythathavearelativelyhighcontributiontothemodel,as
wellastheirupstreamgenes(S3File).
Conclusion
Weidentifiedaspecificsetofgenesintheinterferonpathwaythatareunderexpressedinhigh–
riskmelanoma.Thisbiologicalsignature,togetherwiththeoverexpressionofothergenesinthe
mitoticcellcyclepathway,predictstheprognosisofmelanomawithrelativelyhighaccuracy.
PLOSONE|DOI:10.1371/journal.pone.0170025 January23,2017 8/15
IFNPathwayandMelanomaPrognosis
SupportingInformation
S1File.TheeigengenevaluesintheTCGAdatasetusefulforreproducingtheresults.
(XLS)
S2File.Thefrequencyofselectedmodulesby5–foldcross-validationofthepenalizedCox
model.
(XLS)
S3File.Thegenelistscorrespondingtothethreeselectedmodules.Thegenesymbol,Entrez
ID,andtheweightofeachgeneinthecorrespondingselectedmoduleisreported.
(XLS)
S4File.Supplementarycode.OurresultscanbeconvenientlyreproducedusingtheseR
scripts.Uncompressthetarballfile,installthepackagesmentionedinthesettings.R
script,andthensourcetherunall.Rscript.ThiscodeworksonUnix(LinuxorMacOSX).
DatawillbedownloadedfromTCGAandtheresultswillbesavedinthecurrentdirectory.A
desktopcomputerwitha2.8GHzCPUand8GBofmemorywillreproduceallresultsinless
thananhour.
(ZIP)
S1Fig.Thedistributionofthesizeofmodules.
(PNG)
S2Fig.Theoverrepresentationanalysisontheselectedmodules.
(PNG)
S3Fig.TheKaplan–Meiersurvivalcurvesforthevalidationdatasets.Colorsaresimilarto
(Fig1).IntheLUNDdataset,includingtheinterferonpathwaymoduleresultsinbetterpredic-
tionsofthesurvivaltime(a)withamoresignificantp-valueof2×10−3comparedtoanAFT
modelthatusesonlytwomodules(b).Similarly,intheLEEDSdataset,themodelpredictsthe
survivalratebetterwhentheinterferonpathwaymoduleisincluded(c)comparedtoamodel
thatusesonlytwomodules(d).
(PNG)
Acknowledgments
ThisstudywassupportedbytheNationalScienceFoundation(NSF)throughtheResearch
ExperiencesforUndergraduates(REU)program(CNS1358939),andalsothroughtheinfra-
structuregrant(CRI1305302).WeusedgeneexpressiondatageneratedbytheTCGA
ResearchNetwork:http://cancergenome.nih.gov/.
AuthorContributions
Conceptualization:HZ.
Datacuration:HZ.
Formalanalysis:AZPY.
Fundingacquisition:HZ.
Investigation:AZPY.
Methodology:HZ.
PLOSONE|DOI:10.1371/journal.pone.0170025 January23,2017 9/15
IFNPathwayandMelanomaPrognosis
Projectadministration:HZ.
Resources:HZ.
Software:AZPY.
Supervision:HZ.
Validation:HZAZPY.
Visualization:AZPY.
Writing–originaldraft:HZAZPY.
Writing–review&editing:HZAZPY.
References
1. SiegelRL,MillerKD,JemalA.Cancerstatistics,2015.CA:acancerjournalforclinicians.2015;65
(1):5–29.
2. BalchCM,GershenwaldJE,SoongSj,ThompsonJF,DingS,ByrdDR,etal.Multivariateanalysisof
prognosticfactorsamong2,313patientswithstageIIImelanoma:comparisonofnodalmicrometas-
tasesversusmacrometastases.JournalofClinicalOncology.2010;28(14):2452–2459.doi:10.1200/
JCO.2009.27.1627PMID:20368546
3. BalchCM,GershenwaldJE,SoongSj,ThompsonJF,AtkinsMB,ByrdDR,etal.Finalversionof2009
AJCCmelanomastagingandclassification.Journalofclinicaloncology.2009;27(36):6199–6206.
doi:10.1200/JCO.2009.23.4799PMID:19917835
4. GeramiP,CookRW,WilkinsonJ,RussellMC,DhillonN,AmariaRN,etal.Developmentofaprognos-
ticgeneticsignaturetopredictthemetastaticriskassociatedwithcutaneousmelanoma.ClinicalCan-
cerResearch.2015;21(1):175–183.doi:10.1158/1078-0432.CCR-13-3316PMID:25564571
5. SondakVK,MessinaJL.PredictionisDifficult,EspeciallyAbouttheFuture:ClinicalPrognosticTools
inMelanoma.Annalsofsurgicaloncology.2016;p.1–3.
6. GimottyPA,ElderDE,FrakerDL,BotbylJ,SellersK,ElenitsasR,etal.Identificationofhigh-risk
patientsamongthosediagnosedwiththincutaneousmelanomas.JournalofClinicalOncology.2007;
25(9):1129–1134.doi:10.1200/JCO.2006.08.1463PMID:17369575
7. MortonDL,ThompsonJF,CochranAJ,MozzilloN,ElashoffR,EssnerR,etal.Sentinel-nodebiopsy
ornodalobservationinmelanoma.NewEnglandJournalofMedicine.2006;355(13):1307–1317.doi:
10.1056/NEJMoa060992PMID:17005948
8. EkmekciogluS,DaviesMA,TaneseK,RoszikJ,Shin-SimM,BassettRL,etal.InflammatoryMarker
TestingIdentifiesCD74ExpressioninMelanomaTumorCells,anditsExpressionAssociateswith
FavorableSurvivalforStageIIIMelanoma.ClinicalCancerResearch.2016;p.clincanres–2226.doi:
10.1158/1078-0432.CCR-15-2226PMID:26783288
9. KimbroughCW,EggerME,McMastersKM,StrombergAJ,MartinRC,PhilipsP,etal.MolecularStag-
ingofSentinelLymphNodesIdentifiesMelanomaPatientsatIncreasedRiskofNodalRecurrence.
JournaloftheAmericanCollegeofSurgeons.2016;222(4):357–363.doi:10.1016/j.jamcollsurg.2015.
12.042PMID:26875070
10. WeissSA,HannifordD,HernandoE,OsmanI.Revisitingdeterminantsofprognosisincutaneousmel-
anoma.Cancer.2015;121(23):4108–4123.doi:10.1002/cncr.29634PMID:26308244
11. FrancisP,NamløsHM,Mu¨llerC,Ede´nP,FernebroJ,BernerJM,etal.Diagnosticandprognostic
geneexpressionsignaturesin177softtissuesarcomas:hypoxia-inducedtranscriptionprofilesignifies
metastaticpotential.BMCgenomics.2007;8(1):1.doi:10.1186/1471-2164-8-73
12. PaikS,ShakS,TangG,KimC,BakerJ,CroninM,etal.Amultigeneassaytopredictrecurrenceof
tamoxifen-treated,node-negativebreastcancer.NewEnglandJournalofMedicine.2004;351
(27):2817–2826.doi:10.1056/NEJMoa041588PMID:15591335
13. ColmanH,ZhangL,SulmanEP,McDonaldJM,ShooshtariNL,RiveraA,etal.Amultigenepredictor
ofoutcomeinglioblastoma.Neuro-oncology.2009;p.nop007.
14. GordonGJ,JensenRV,HsiaoLL,GullansSR,BlumenstockJE,RichardsWG,etal.Usinggene
expressionratiostopredictoutcomeamongpatientswithmesothelioma.JournaloftheNationalCan-
cerInstitute.2003;95(8):598–605.doi:10.1093/jnci/95.8.598PMID:12697852
PLOSONE|DOI:10.1371/journal.pone.0170025 January23,2017 10/15