Hi,

I am trying to compare two result lists from GOstats. The reason for doing this is that I did a proteomics experiment in which we established a new sample preparation protocol and I now want to kind of show that with this new method we do not enrich/deplete for certain proteins/protein functions. The obvious way, comparing the protein lists and check the overlap isn`t very telling as there is redundancy in the proteome, i.e. protein X does not neccessarily have the same ID in both sets, and further, both sets contain different amount of entries. So I thought doing a GO enrichment analysis to check whether certain catergories are enriched in one list versus the other would be appropriate.

So now I have the lists of the two sets and I wanted to do a simple statistical test to compare them GO term by GO term. As this is a count based list I wanted to use a GTest (or Chi square) . But as the input lists have different sizes, I would naturally expect less Counts per GO term in the smaller data set (GO.list.1) compared to the larger list (GO.list.2), meaning that the simple GTest:

>GO.list<-merge(GO.list.1, GO.list.2, by=2, all=F)
>colnames(GO.list)
 [1] "GOMFID"      "Pvalue.x"    "OddsRatio.x" "ExpCount.x"  "Count.x"    "Size.x"
 [7] "Term.x"     "test.x"      "ontology.x"  "subset.x"  "Pvalue.y"    "OddsRatio.y"  
[13] "ExpCount.y"  "Count.y"     "Size.y"   "Term.y"      "test.y"      "ontology.y"  
[19]  "subset.y"   

>GO.list$pvalue<-sapply(1:nrow(GO.list), 
 function(i) GTest(x=as.numeric(unlist(GO.list[i,c(5,14)]))))$pvalue)

will bias my result as it assumes equal proportions.

As the GOstats result file give you all kind of information (based on the input data) I thought I`d use these informations to calculate my expected proportions. My first thought now was to normalise the counts per GO term on the relative proportions of the ExpCounts given in the GOstats result file, i.e. columns 4  and 13 of the merged result file.

>GO.list<-merge(GO.list.1, GO.list.2, by=2, all=F)
>colnames(GO.list)
 [1] "GOMFID"      "Pvalue.x"    "OddsRatio.x" "ExpCount.x"  "Count.x"    "Size.x"
 [7] "Term.x"     "test.x"      "ontology.x"  "subset.x"  "Pvalue.y"    "OddsRatio.y"  
[13] "ExpCount.y"  "Count.y"     "Size.y"   "Term.y"      "test.y"      "ontology.y"  
[19]  "subset.y"   

>p1<-GO.list$ExpCount.x/(GO.list$ExpCount.y+GO.list$ExpCount.x)
>p2<-GO.list$ExpCount.y/(GO.list$ExpCount.y+GO.list$ExpCount.x)

>GO.list$pvalue<-sapply(1:nrow(GO.list), 
 function(i) GTest(x=as.numeric(unlist(GO.list[i,c(5,14)]))), p=c(p1,p2))$pvalue)

But I don`t think this is the correct way to do it either, because if I run GOstats on a merged list (i.e. merging my two protein in put lists) then the ExpCounts are (somewhat expectedly) not the sum of the two input lists.

So I was wondering if I could use the Count columns as a proxy for calculating the expected proportions:

>GO.list<-merge(GO.list.1, GO.list.2, by=2, all=F)
>colnames(GO.list)
 [1] "GOMFID"      "Pvalue.x"    "OddsRatio.x" "ExpCount.x"  "Count.x"    "Size.x"
 [7] "Term.x"     "test.x"      "ontology.x"  "subset.x"  "Pvalue.y"    "OddsRatio.y"  
[13] "ExpCount.y"  "Count.y"     "Size.y"   "Term.y"      "test.y"      "ontology.y"  
[19]  "subset.y"   

>p1<-((sum(GO.list$Count.x)/sum(GO.list$Count.y))*0.5)/(((sum(GO.list$Count.x)/sum(GO.list$Count.y))*0.5)+0.5) 
>p2<-((sum(GO.list$Count.y)/sum(GO.list$Count.x))*0.5)/(((sum(GO.list$Count.y)/sum(GO.list$Count.x))*0.5)+0.5)

#in my case: p1: 0.4176183
#p2: 0.5823817

>GO.list$pvalue<-sapply(1:nrow(GO.list), 
 function(i) GTest(x=as.numeric(unlist(GO.list[i,c(5,14)]))), p=c(p1,p2))$pvalue)

Does anybody have any suggestion if this is any way correct, or should I do something completely different?

Cheers

Sven

I wouldn't do it that way. In one set of data you have some evidence that a set of proteins that make up a group are being perturbed, and you want to test for evidence that the same group of proteins is (not) being affected in a different experiment. The gene set testing framework is the canonical way of doing that, and it has the added benefit that you don't necessarily need to measure the exact same things in each experiment (although you would like a pretty good overlap). The limma package has a bunch of different tests, but for your purposes I would imagine that roast would be applicable.