Συγκρίνοντας λίστες γονιδίων

Συχνά θελουμε να συγκρίνουμε τα αποτελέσματα από δύο πειράματα ή παραπάνω ώστε να δούμε αν δύο ασθένειες για παράδειγμα έχουν ως βάση τα ίδια γονίδια.

Μεθοδολογία

Εφαρμόζουμε την μεθοδολογία που ήδη έχουμε μάθει σε κάθε dataset με στόχο την ευρεση γονιδίων που διαφοροποιούνται σε κάθε ένα από τα dataset. Ας υποθέσουμε λοιπόν ότι έχουμε τα ακόλουθα datasets που αποθηκευονται στις μεταβλητές a και b.

## !!! the next datasets contains many null values, specified by 'null'. They are NAs actuallys
a = read.table("GDS4193.clean", h=TRUE, na.strings = 'null')
b = read.table("GDS5272.clean", h=TRUE, na.strings = 'null')

ids = complete.cases(a) & complete.cases(b)

d1 = a[ids,3:ncol(a)]
d2 = b[ids,3:ncol(b)]

genes = a[ids,2]

d1.group = factor(rep(c("patient", "control"), c(17, 9)), levels = c("control", "patient"))
d2.group = factor(rep(c("patient", "control"), c(5,5)), levels = c("control", "patient"))

DIfferential Expression Analysis

library(limma)
dge = function(d, f){
  mm = model.matrix(~0+f)
  #print(mm)
  colnames(mm) = c("control", "patient")
  linModel = lmFit(object = d, design = mm)
  contrs = makeContrasts(patient - control, levels = mm)
  fit.contrs = contrasts.fit(linModel, contrs)
  efit = eBayes(fit.contrs)
  dge.genes = topTable(efit, adjust.method = "fdr", number = nrow(d))
  return(dge.genes)
}


dge2 = function(d, f){
  pvalues = apply(d, 1, function(x){
    lmod = lm(x~f)
    pvalue = summary(lmod)$coefficients[2,4]
    return(pvalue)
  })
  return(pvalues)
}

d1.table = dge(d1, d1.group)
d2.table = dge(d2, d2.group)

Παρατηρήστε ότι σε διορθωμένα pvalues, δεν υπάρχει κανένα γονίδιο στο πρώτο dataset που να έχει p-value μικρότερο από 0.05. Αρα δεν έχει νοημα να βρούμε κοινά γονίδια αφού στο πρωτο dataset δεν υπάρχουν γονίδια που διαφοροποιούν την έκφραση τους.

Όμως, για το μάθημα, ας αλλάξουμε το κατώφλι και να έχουμε μερικά γονιδια. Έστω λοιπόν ότι θέτουμε το επίπεδο σημαντικότητας με βάση το p-value (κι οχι το διορθωμένο p-value sto 0.01). Τότε θα είχαμε την εξής ανάλυση

thr = 0.01
col = "P.Value"
head(d1.table)
##            logFC  AveExpr         t      P.Value adj.P.Val         B
## 21478  0.3469212 4.634833  4.601141 9.159436e-05 0.7907169 -1.198435
## 45354 -0.3175116 3.593118 -4.333349 1.869210e-04 0.7907169 -1.484166
## 24706 -0.3700350 7.210280 -4.327807 1.896942e-04 0.7907169 -1.490144
## 19181 -0.4619599 7.576772 -4.316609 1.954233e-04 0.7907169 -1.502229
## 37311 -0.3433600 9.524445 -4.217453 2.542434e-04 0.7907169 -1.609651
## 27461 -0.3874054 8.467904 -4.207200 2.612459e-04 0.7907169 -1.620799
setgenes1 = row.names(d1.table) [d1.table[,col] < thr]
setgenes2 = row.names(d2.table) [d2.table[,col] < thr]


commonGenes = intersect(setgenes1, setgenes2)
k = length(commonGenes)
k1 = length(setgenes1)
k2 = length(setgenes2)
n = nrow(d1.table)
hyper.p.value = phyper(q=k-1, m = k1, n=n-k1, k = k2, lower.tail=FALSE)
p1 = hyper.p.value

## as fisher's exact test
cont.table = matrix(c(k, k2-k, k1-k, n-k1-k2+k), nrow=2)
p2 = fisher.test(x=cont.table, alternative = "greater")$p.value
## binomial test
p.success = k1/n
p3 = binom.test(x=k, n=k2, p=p.success, alternative = "greater")$p.value
## binomial test with pbinom
p4 = pbinom(q = k-1, size = k2, prob = p.success, lower.tail = FALSE)
print(c(hyper=p1, fisher=p2, binTest=p3, pbinom=p4))
##      hyper     fisher    binTest     pbinom 
## 0.03333899 0.03333899 0.03784243 0.03784243

Υπάρχει και μία εναλλακτική προσσέγγιση με βάση την οποία δεν χρειάζεται να θέσουμε πιο κομμάτι της λίστας των γονιδίων θεωρούμε ως διαφοροποιημένο. Το παίρνουμε όλο από την αρχή και το ελέγχουμε τμηματικά

## it's important that gene names are sorted based on the p-value
setgenes1 = row.names(d1.table) 
setgenes2 = row.names(d2.table) 

getAllpvalues = function(s1, s2){
  pvalues = vector(mode = 'numeric', length = length(s1))
  fcs = vector(mode = 'numeric', length = length(s1))
  for(i in 1:length(s1)){
    ts1 = s1[1:i]
    ts2 = s2[1:i]
    commonGenes = intersect(ts1, ts2)
    k = length(commonGenes)
    k1 = length(ts1)
    k2 = length(ts2)
    n = nrow(d1.table)
    hyper.p.value = phyper(q=k-1, m = k1, n=n-k1, k = k2, lower.tail=FALSE)
    pvalues[i] = hyper.p.value
    fc = (k1/n)/(k/k2) 
    fcs[i] = fc
  }
  return(list(pvalues = pvalues, fcs = fcs))
}
res = getAllpvalues(setgenes1, setgenes2)
all.size.pvalues = res$pvalues
fc = res$fcs

{xlim=c(1,3000)
plot(1:length(all.size.pvalues), all.size.pvalues, type='l', axes=FALSE, xlab="Size of top", ylab="", xlim=xlim)
axis(side = 1, labels = TRUE)
axis(side = 2)
mtext(side=2, "P-value", 2)
par(new=TRUE)
plot(1:length(all.size.pvalues), fc,  type='l', axes=FALSE, ylab="", col='red', xlim=xlim, xlab="")
axis(side = 4, labels = TRUE)
mtext(side=4, "Fold Change", 2)}

## the same plot in log-scale
##plot(1:length(all.size.pvalues), log10(all.size.pvalues), type='l')

Φαίνεται να υπάρχει ένα ενδιαφέρον (τοπικό) ελάχιστο στην αρχή που ειναι ενδιαφέρον. Αυτό ειναι το εξής:

min(all.size.pvalues[1:500])
## [1] 0.008160177
int.id = which.min(all.size.pvalues[1:500])

common.genes = intersect(setgenes1[1:int.id], setgenes2[1:int.id])
as.numeric(common.genes)
##  [1] 11454 19659 37741 35053 19088 33266 35052 10474 13887 28544 29205
## [12] 36128 44478
my.genes = a[as.numeric(common.genes), 2]
my.genes
##  [1] ST14    RGS13   SLC44A1 SMDT1   SMAGP   GALNT2  SMDT1   EIF5B  
##  [9] IFI44L  SEMA4A  AMIGO3  RFX2    CCDC85A
## 30806 Levels: 1552563_a_at 1552829_at 1552867_at 1552961_at ... ZZZ3
write.table(my.genes, file="", quote=FALSE, row.names = FALSE, col.names=FALSE)
## ST14
## RGS13
## SLC44A1
## SMDT1
## SMAGP
## GALNT2
## SMDT1
## EIF5B
## IFI44L
## SEMA4A
## AMIGO3
## RFX2
## CCDC85A