Question

Analysis of pathways using edgeR

0

Entering edit mode

fawazfebin ▴ 60

@fawazfebin-14053

Last seen 3.7 years ago

Hi

I was analysing public RNA-Seq datasets from GEO for arriving at the most affected pathways /GO terms using the edgeR program. My objective was to determine whether Ca signalling pathway had a considerable effect in the corresponding phenotype. I did a DE analysis with edgeR and arrived at result in which the Ca signalling pathway was one among the top (topKEGG) hundred pathways :

(No.80) Pathway N Up Down P.Up P.Down

path:hsa04020 Calcium signaling pathway 79 8 5 3.534894e-02 3.626250e-01 (at 2hr)

(No.72) Pathway N Up Down P.Up P.Down

path:hsa04020 Calcium signaling pathway 79 14 12 1.795998e-02 3.278448e-01 (at 4hr)

What can I conclude on this pathway regarding its effect ? Thanks in advance.

rnaseq differential expression edgeR kegg • 2.5k views

ADD COMMENT • link updated 6.1 years ago by Gordon Smyth 50k • written 6.1 years ago by fawazfebin ▴ 60

score 3 · Answer 1 · 2018-03-09

3

Entering edit mode

Gordon Smyth 50k

@gordon-smyth

Last seen 37 minutes ago

WEHI, Melbourne, Australia

Well, the analysis shows that your up-regulated genes are moderately enriched for Ca signally pathway genes (P = 0.035 at 2 hours and P = 0.018 at 4 hours).

Down-regulated genes are not enriched (P = 0.36 at 2 hours and P = 0.33 at 4 hours).

Further interpretation is a biological question, and is up to you!

If you are only interested in the KEGG Calcium Signalling pathway, you might consider doing a fry test for this pathway:

> KEGG <- getGeneKEGGLinks(species.KEGG = "hsa")
> HSA0420 <- KEGG$GeneID[ KEGG$PathwayID=="path:hsa04020" ]
> fry(y, index=HSA0420, design=design, contrast=yourcontrast)

ADD COMMENT • link 6.1 years ago Gordon Smyth 50k

0

Entering edit mode

Thank you for the valuable guidance ,Sir. Is 'y' the DGElist object in the last command of fry test ? I get the following error while running it :

> fry(y.DEX, index=HSA0420, design=design.DEX, contrast="ResistantvsControl.2hr")
Error in .zscoreDGE(y = y, design = design, contrast = contrast) :
  contrast doesn't match any column of design

ADD REPLY • link 6.1 years ago fawazfebin ▴ 60

1

Entering edit mode

Well, the error message is self-explanatory, is it not?

Yes, 'y' is the DGEList that you used for your edgeR analysis.

'contrast' is just the same contrast that you used for your edgeR DE analysis. It has the same meaning in fry() as it has in glmLRT() or glmQLFTest(). You must already know how to specify the contrast in edgeR, otherwise you couldn't have done a DE analysis in the first place.

ADD REPLY • link 6.1 years ago Gordon Smyth 50k

0

Entering edit mode

Yes Sir, I had used the same DGElist object and the same contrast used in DE analysis here..but couldnt understand why the error came still..Thanks

ADD REPLY • link 6.1 years ago fawazfebin ▴ 60

1

Entering edit mode

Type

colnames(design.DEX)

to see what the columns of the design matrix are.

You must have changed at least one of the design matrix or the contrast when running fry.

ADD REPLY • link 6.1 years ago Gordon Smyth 50k

0

Entering edit mode

> colnames(design.DEX)
[1] "control"       "resistant.2hr" "resistant.4hr"

> DEX.contrasts
               Contrasts
   Levels          resistant.2hrvscontrol   resistant.4hrvscontrol
  control                           -1                     -1
  resistant.2hr                      1                      0
  resistant.4hr                      0                      1

>  fry(y.DEX, index=HSA0420, design=design.DEX, contrast="resistant.2hrvscontrol")
Error in .zscoreDGE(y = y, design = design, contrast = contrast) :
  contrast doesn't match any column of design

I am confused at the error. What could be wrong here? Thanks in advance.

ADD REPLY • link 6.1 years ago fawazfebin ▴ 60

1

Entering edit mode

Perhaps you intended to do this:

fry(y.DEX, index=HSA0420, design=design.DEX,
    contrast=DEX.contrasts[,"resistant.2hrvscontrol"])

ADD REPLY • link 6.1 years ago Gordon Smyth 50k

0

Entering edit mode

Yes, Sir. Great thanks.

ADD REPLY • link 6.1 years ago fawazfebin ▴ 60

0

Entering edit mode

> HSA0420 <- KEGG$GeneID[ KEGG$PathwayID=="path:hsa04020" ]

> fry(y.DEX, index=HSA0420, design=design.DEX,
+     contrast=DEX.contrasts[,"resistant.2hrvscontrol"])
     NGenes Direction PValue PValue.Mixed
set1      0        Up    NaN            0

Can the output of the above fry function be interpreted as the calcium signalling pathway is over represented ? Thanks in advance.

ADD REPLY • link 6.0 years ago fawazfebin ▴ 60

1

Entering edit mode

Obviously not. fry is telling you that there are no genes in the set, and so it's returning a nonsense p-value of NaN. Are you sure that the row names of y.DEX are Entrez IDs?