Hi,

Being new to R, I am following the protocol for GAGE from DESeq in the GAGE vignette, but I don't think it fits the setup of my experiment.

My read count table is as follow:

                   REF SAMP1 SAMP2 SAMP3
ENSG00000000001.6    0     1     1     4
ENSG00000000008.3     1     6     3     1
ENSG00000000434.1    15    14     3     4
ENSG00000000417.5    80   279   108   349
ENSG00000000230.1  279    69    23   373
ENSG00000000876.8     2     4     2     0​

I need to compare the fold changes independently for REF vs SAMP1, then REF vs SAMP2, then REF vs SAMP3.

But would need to map them on the same map further with GAGE.

Starting from the DESeq fork:

# add the group info
grp.idx <- factor(c("Reference","Sample","Sample","Sample"))
cds <- newCountDataSet(cnts, grp.idx)

# normalise library size
cds <- estimateSizeFactors(cds)

# variance estimation
cds <- estimateDispersions(cds,method="blind",sharingMode="fit-only")

At this point should I insert a for loop to process the 3 pairs of REF/SAMP* ?

And then, keep going with the workflow:

##### BEGGINING FOR LOOP HERE ??? #####

# Calling differential expression
deseq.res <- nbinomTest(cds, "Reference","Sample")
deseq.fc=deseq.res$log2FoldChange

# Remove the '-Inf' an 'Inf' fold change values
deseq.fc[deseq.fc>10]=10
deseq.fc[deseq.fc< -10]=-10

exp.fc=deseq.fc

# Plot dispersion
# Plot the log2 fold change against mean normalised counts
# Plot Histogram p-values
# Table of differential expressed genes

##### END FOR LOOP HERE ??? #####

### GAGE fork ###
##### 'Re-assemble' the 3 fold change data frames HERE ??? #####
data(kegg.gs)
cn=colnames(cnts)
ref.idx=grep('REF', cn)
samp.idx=grep('SAMP',cn)

# to see both stats / up- / down-expressed genes
fc.kegg.p <- gage(exp.fc_reassembled, gsets=kegg.gs, ref=ref.idx, samp=samp.idx, compare="unpaired", same.dir=TRUE)

### use 'exp.fc_reassembled' as data frame for PATHVIEW and GO Enrichment forks ###

QUESTION #1: is it the right logic to insert a for loop for doing these pairwise comparisons? If yes, how and where do I insert the loop?

QUESTION #2: is it the right way to 're-assemble' the 3 fold change data frames together? If yes, how to do that?

Sorry if they are basic questions.

Thanks !