Differential Gene Expression
Pavlos Pavlidis
4/8/2021
Εισαγωγή
Στο σημερινό μάθημα θα εμβαθύνουμε λίγο στην χρήση γραμμικών μοντέλων για DGE analysis. (Differential Gene Expression analysis).
Ας θεωρήσουμε το dataset που χρησιμοποιούμε στα τελευταία μαθήματα
a <- read.table("GDS3709.soft.clean", h=TRUE)
des <- read.table("GDS3709.desc")
##b <- a[, 3:(39+2)]
data = a[,3:ncol(a)]
gender = rep(c("F", "M"), c(39, 40))
smoking = rep( rep( c("s", "n"), 2), c(19, 20, 20, 20))Βλέπουμε ότι έχουμε δύο παράγοντες το κάπνισμα (smoking) και το φύλο (gender). Οι παράγοντες αυτοι επηρεάζουν την έκφραση γονιδίων και επομένως μπορεί να θέλουμε να βρούμε τα γονίδια των οποίων η έκφραση επηρεάζεται από τους παράγοντες αυτούς.
Το πρώτο βήμα ειναι να προσδιορίσουμε το model.matrix το οποίο στις γραμμές έχει τα δείγματα και στις στήλες τα treatments. Ουσιαστικά το model.matrix μας λέει πως σχετίζονται τα δείγματα μας με τα treatments, ποιο ειναι το μοντέλο που στήνουμε για να μοντελοποιήσουμε την έκφραση γονιδίων. Πολύ σημαντικές πληροφορίες βρίσκονται εδώ:
Πάμε λοιπόν να δούμε πως εφαρμόζεται αυτό στο δικό μας dataset.
mm1 = model.matrix(~gender+smoking)
mm1## (Intercept) genderM smokings
## 1 1 0 1
## 2 1 0 1
## 3 1 0 1
## 4 1 0 1
## 5 1 0 1
## 6 1 0 1
## 7 1 0 1
## 8 1 0 1
## 9 1 0 1
## 10 1 0 1
## 11 1 0 1
## 12 1 0 1
## 13 1 0 1
## 14 1 0 1
## 15 1 0 1
## 16 1 0 1
## 17 1 0 1
## 18 1 0 1
## 19 1 0 1
## 20 1 0 0
## 21 1 0 0
## 22 1 0 0
## 23 1 0 0
## 24 1 0 0
## 25 1 0 0
## 26 1 0 0
## 27 1 0 0
## 28 1 0 0
## 29 1 0 0
## 30 1 0 0
## 31 1 0 0
## 32 1 0 0
## 33 1 0 0
## 34 1 0 0
## 35 1 0 0
## 36 1 0 0
## 37 1 0 0
## 38 1 0 0
## 39 1 0 0
## 40 1 1 1
## 41 1 1 1
## 42 1 1 1
## 43 1 1 1
## 44 1 1 1
## 45 1 1 1
## 46 1 1 1
## 47 1 1 1
## 48 1 1 1
## 49 1 1 1
## 50 1 1 1
## 51 1 1 1
## 52 1 1 1
## 53 1 1 1
## 54 1 1 1
## 55 1 1 1
## 56 1 1 1
## 57 1 1 1
## 58 1 1 1
## 59 1 1 1
## 60 1 1 0
## 61 1 1 0
## 62 1 1 0
## 63 1 1 0
## 64 1 1 0
## 65 1 1 0
## 66 1 1 0
## 67 1 1 0
## 68 1 1 0
## 69 1 1 0
## 70 1 1 0
## 71 1 1 0
## 72 1 1 0
## 73 1 1 0
## 74 1 1 0
## 75 1 1 0
## 76 1 1 0
## 77 1 1 0
## 78 1 1 0
## 79 1 1 0
## attr(,"assign")
## [1] 0 1 2
## attr(,"contrasts")
## attr(,"contrasts")$gender
## [1] "contr.treatment"
##
## attr(,"contrasts")$smoking
## [1] "contr.treatment"Παρατήρηση 1
Δεν φαίνονται όλες οι πιθανές στήλες στο model.matrix
Αυτό γίνεται επειδή οι στήλες στο model.matrix πρέπει να ειναι ανεξάρτητες μεταξύ τους. Στο συγκεκριμένο παράδειγμα το Intercept περιλαμβάνει τα άτομα που ειναι F,n-s, δηλαδή Θηλυκά που δεν καπνίζουν. Αυτή η ομάδα χρησιμοποιείται σαν CONTROL group, και αυτό ειναι ‘κατά σύμβαση’. Μπορούμε φυσικά να το τροποποιήσουμε αυτό.
Τώρα μπορούμε να βρούμε την (μέση) επίδραση του να ειναι κάποιο άτομο θηλυκό στην έκφραση του γονιδίου 1
expr = as.numeric(data[1,])
expr## [1] 1124.0 1196.0 982.8 1075.0 1114.0 1302.0 1279.0 1210.0 1199.0 1172.0
## [11] 1245.0 1307.0 1263.0 1282.0 1191.0 1189.0 1202.0 1372.0 1343.0 1151.0
## [21] 1447.0 1362.0 1033.0 1307.0 1295.0 1198.0 1042.0 1105.0 884.7 1267.0
## [31] 1358.0 1270.0 1375.0 1181.0 1176.0 1196.0 1177.0 1227.0 1391.0 1172.0
## [41] 1271.0 1416.0 1123.0 1120.0 981.3 1115.0 1272.0 1234.0 1252.0 1202.0
## [51] 1199.0 1009.0 1385.0 1230.0 1192.0 1114.0 1218.0 1235.0 1037.0 1350.0
## [61] 1341.0 1189.0 1209.0 1521.0 1669.0 1205.0 1111.0 1105.0 1019.0 1234.0
## [71] 1080.0 1402.0 1182.0 1326.0 1105.0 1361.0 1309.0 1315.0 1110.0lm(formula = expr~gender)##
## Call:
## lm(formula = expr ~ gender)
##
## Coefficients:
## (Intercept) genderM
## 1217.705 5.302#αν ειχαμε
gender = relevel(factor(gender), "F")
linmod = lm(formula = expr~gender)
summary(linmod)##
## Call:
## lm(formula = expr ~ gender)
##
## Residuals:
## Min 1Q Median 3Q Max
## -333.01 -96.86 -15.71 85.14 445.99
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1217.705 20.820 58.487 <2e-16 ***
## genderM 5.302 29.260 0.181 0.857
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 130 on 77 degrees of freedom
## Multiple R-squared: 0.0004263, Adjusted R-squared: -0.01256
## F-statistic: 0.03284 on 1 and 77 DF, p-value: 0.8567Ο τροπος με το intercept (θέτοντας δηλαδή ενα group σαν βάση) ειναι πολλές φορές άνετος για να βρούμε gene expression differences. Όμως μερικές φορές ειναι καλύτερα να ορίζουμε την σύγκριση που θέλουμε να μελετήσουμε.
library(limma)
mm = model.matrix(~0+gender+smoking)
mm## genderF genderM smokings
## 1 1 0 1
## 2 1 0 1
## 3 1 0 1
## 4 1 0 1
## 5 1 0 1
## 6 1 0 1
## 7 1 0 1
## 8 1 0 1
## 9 1 0 1
## 10 1 0 1
## 11 1 0 1
## 12 1 0 1
## 13 1 0 1
## 14 1 0 1
## 15 1 0 1
## 16 1 0 1
## 17 1 0 1
## 18 1 0 1
## 19 1 0 1
## 20 1 0 0
## 21 1 0 0
## 22 1 0 0
## 23 1 0 0
## 24 1 0 0
## 25 1 0 0
## 26 1 0 0
## 27 1 0 0
## 28 1 0 0
## 29 1 0 0
## 30 1 0 0
## 31 1 0 0
## 32 1 0 0
## 33 1 0 0
## 34 1 0 0
## 35 1 0 0
## 36 1 0 0
## 37 1 0 0
## 38 1 0 0
## 39 1 0 0
## 40 0 1 1
## 41 0 1 1
## 42 0 1 1
## 43 0 1 1
## 44 0 1 1
## 45 0 1 1
## 46 0 1 1
## 47 0 1 1
## 48 0 1 1
## 49 0 1 1
## 50 0 1 1
## 51 0 1 1
## 52 0 1 1
## 53 0 1 1
## 54 0 1 1
## 55 0 1 1
## 56 0 1 1
## 57 0 1 1
## 58 0 1 1
## 59 0 1 1
## 60 0 1 0
## 61 0 1 0
## 62 0 1 0
## 63 0 1 0
## 64 0 1 0
## 65 0 1 0
## 66 0 1 0
## 67 0 1 0
## 68 0 1 0
## 69 0 1 0
## 70 0 1 0
## 71 0 1 0
## 72 0 1 0
## 73 0 1 0
## 74 0 1 0
## 75 0 1 0
## 76 0 1 0
## 77 0 1 0
## 78 0 1 0
## 79 0 1 0
## attr(,"assign")
## [1] 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$gender
## [1] "contr.treatment"
##
## attr(,"contrasts")$smoking
## [1] "contr.treatment"cm = makeContrasts(genderF-genderM, levels=mm)
linFit = lmFit(data, mm)
fitc = contrasts.fit(linFit, cm)
fitc = eBayes(fitc)
topTable(fitc)## logFC AveExpr t P.Value adj.P.Val B
## 16147 -168.812110 125.27684 -12.68127 1.670634e-20 5.349121e-16 29.68518
## 11358 -2790.748987 1903.62709 -12.64290 1.956697e-20 5.349121e-16 29.57238
## 31009 672.599117 425.56608 12.36609 6.150566e-20 1.120941e-15 28.75114
## 13857 -189.458827 122.85527 -12.10697 1.811345e-19 2.475882e-15 27.97046
## 23518 298.884077 191.36570 11.86982 4.900215e-19 4.637628e-15 27.24590
## 13858 -122.842575 86.18557 -11.86083 5.089304e-19 4.637628e-15 27.21823
## 20513 -9.937403 13.74244 -11.63110 1.342739e-18 1.048775e-14 26.50697
## 16510 -82.364623 70.56278 -11.20342 8.293826e-18 5.143349e-14 25.15963
## 37747 -116.008617 85.39362 -11.18581 8.943253e-18 5.143349e-14 25.10350
## 33848 2250.830727 1443.55747 11.17399 9.407131e-18 5.143349e-14 25.06584Παραδείγματα
Παράδειγμα 1 - Προσομοιωμένα δεδομένα
Στο παράδειγμα αυτό δεν θα χρησιμοποιήσουμε αληθινά δεδομένα αλλά προσομοιωμένα, δηλαδή θα τα κατασκευάσουμε στον υπολογιστή.
## Data generation
nrow=10
ncol=100
grp1 = factor(rep(c(0,1), each=ncol/2))
grp2 = factor(rep(c(0,1), ncol/2))
simdata = matrix(rnorm(nrow*ncol, 100, 30), nrow=nrow, ncol=ncol)
simdata[,((ncol/2+1):ncol)] = matrix(rnorm(nrow*ncol/2, 150, 30), nrow=nrow, ncol=ncol/2)
simdata[, seq(from=2, to=ncol, by=2)] = simdata[, seq(from=2, to=ncol, by=2)] + matrix(rnorm(nrow*ncol/2, 80, 30), nrow=nrow, ncol=ncol/2)
## Visualization with a heatmap
library(gplots)##
## Attaching package: 'gplots'## The following object is masked from 'package:stats':
##
## lowessheatmap.2(simdata, Colv=F, Rowv=F, dendrogram='none', trace = 'none')
## Let's build some models
## model 1 -- we try to explain all the variability of the data using just the grp1
sim.linmod.1 = lm(simdata[1,]~grp1)
summary(sim.linmod.1)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 139.56233 6.430580 21.702914 3.268323e-39
## grp11 60.48126 9.094214 6.650521 1.672253e-09## model 2 -- we try to explain all the variability of the data using just the grp2
sim.linmod.2 = lm(simdata[1,]~grp2)
summary(sim.linmod.2)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 140.54894 6.523057 21.54648 5.880243e-39
## grp21 58.50806 9.224996 6.34234 6.991041e-09## model 3 -- we try to explain the variability of the data using ADDITIVELY the grp1 and grp2
## this means that the effect of grp2 is the same for both levels of grp1
sim.linmod.3 = lm(simdata[1,]~0+grp1*grp2)
summary(sim.linmod.3)$coefficients## Estimate Std. Error t value Pr(>|t|)
## grp10 116.06294 6.884015 16.859774 1.892171e-30
## grp11 165.03493 6.884015 23.973646 2.612338e-42
## grp21 46.99879 9.735467 4.827584 5.210780e-06
## grp11:grp21 23.01853 13.768030 1.671883 9.780393e-02## model 4 -- we try to explain the variability of the data using both groups but also assuming that
## the effect of grp2 might be different for each level of grp1
sim.linmod.4 = lm(simdata[1,]~0+grp1*grp2)
summary(sim.linmod.4)$coefficients## Estimate Std. Error t value Pr(>|t|)
## grp10 116.06294 6.884015 16.859774 1.892171e-30
## grp11 165.03493 6.884015 23.973646 2.612338e-42
## grp21 46.99879 9.735467 4.827584 5.210780e-06
## grp11:grp21 23.01853 13.768030 1.671883 9.780393e-02## Next, we will try to see how the coefficients in the models above, are connected to the mean values of each of the groups below. We will do this only for the first row of the data matrix, but anyway for all of the rows it's the same
## Για να δούμε λίγο τους μέσους όρους
grp1.0 = grp1 == 0
grp1.1 = grp1 == 1
grp2.0 = grp2 == 0
grp2.1 = grp2 == 1
## just for the first row
i = 1
mean1.0 = mean(simdata[i,grp1.0])
mean1.1 = mean(simdata[i,grp1.1])
mean2.0 = mean(simdata[i,grp2.0])
mean2.1 = mean(simdata[i,grp2.1])
mean1.0_2.0 = mean(simdata[i,grp1.0&grp2.0])
mean1.0_2.1 = mean(simdata[i,grp1.0&grp2.1])
mean1.1_2.0 = mean(simdata[i,grp1.1&grp2.0])
mean1.1_2.1 = mean(simdata[i,grp1.1&grp2.1])
sim.means = c(m1.0 = mean1.0, m1.1=mean1.1, m2.0 = mean2.0, m2.1 = mean2.1,
m1.0_2.0=mean1.0_2.0,
m1.0_2.1=mean1.0_2.1,
m1.1_2.0=mean1.1_2.0,
m1.1_2.1=mean1.1_2.1)
sim.means## m1.0 m1.1 m2.0 m2.1 m1.0_2.0 m1.0_2.1 m1.1_2.0 m1.1_2.1
## 139.5623 200.0436 140.5489 199.0570 116.0629 163.0617 165.0349 235.0523## Let's plot the data with different colors for each of the group
pointcolors = rep("white", ncol)
pointcolors[grp1.0&grp2.0] = "red"
pointcolors[grp1.0&grp2.1] = "black"
pointcolors[grp1.1&grp2.0] = "blue"
pointcolors[grp1.1&grp2.1] = "green"
## rn adds just some noise on the x-axis, so that the points are not one over the other
## only for visualization reasons
rn = rnorm(ncol, 0, .03)
{plot(as.numeric(grp1)+rn-1, simdata[i,], col=pointcolors)
m1 = lm(simdata[i,]~grp1)
abline(a=m1$coefficients[1], b=m1$coefficients[2], col="lightblue", lw=2)}
sim.means## m1.0 m1.1 m2.0 m2.1 m1.0_2.0 m1.0_2.1 m1.1_2.0 m1.1_2.1
## 139.5623 200.0436 140.5489 199.0570 116.0629 163.0617 165.0349 235.0523m1##
## Call:
## lm(formula = simdata[i, ] ~ grp1)
##
## Coefficients:
## (Intercept) grp11
## 139.56 60.48Στην παραπάνω περίπτωση, θεωρήσαμε τα δεδομένα και προσπαθήσαμε να εξηγήσουμε όλο τις διαφορετικές τιμές έκφρασης με βάση μόνο το grp1, δηλάδή τον πρώτο τρόπο ομαδοποίησης. Λέμε λοιπόν ότι όλες οι τιμές έκφρασης κυμαίνονται γύρω από δύο διαφορετικές τιμές. Η μία υπολογίζεται για τα άτομα που ανήκουν στην ομάδα grp1 == 1 και η άλλη στην ομάδα grp1 == 2
Πράγματι, οι δύο αυτές τιμές ειναι οι μέσοι όροι των ατόμων που ανήκουν στο grp1 == 1 και grp1 == 2. Το βλέπουμε ως εξής αυτό.
print(sim.means)## m1.0 m1.1 m2.0 m2.1 m1.0_2.0 m1.0_2.1 m1.1_2.0 m1.1_2.1
## 139.5623 200.0436 140.5489 199.0570 116.0629 163.0617 165.0349 235.0523print(m1$coefficients)## (Intercept) grp11
## 139.56233 60.48126Παρατηρούμε οτι m1.0 == Intercept και Intercept + grp11 = 135.19376 + 56.26633 == m1.1
Είναι δηλαδή σαν να έχουμε βρει την μέση επίδραση του παράγοντα grp1 στην έκφραση. Και γράφω “μέση επίδραση”, επειδή η επίδραση αυτή αναφέρεται και στα άτομα που έχουν grp2 == 1 και στα άτομα που έχουν grp2 == 2. Με άλλα λόγια, αν ένα έχουμε δύο άτομα που το ένα ανήκει στο grp1 == 1 και το άλλο στο grp1 == 2, προβλέπουμε ότι το δευτερο άτομο θα έχει ανεβασμένη έκφραση κατά 56.26633, άσχετα με το τι γίνεται με το grp2.
Ας βάλουμε τώρα στο παιχνίδι και το grp2. Τότε το μοντέλο γίνεται ως εξής:
m2 = lm(simdata[i,]~grp1+grp2)
summary(m2)##
## Call:
## lm(formula = simdata[i, ] ~ grp1 + grp2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -96.643 -20.303 1.836 24.064 87.870
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 110.308 6.017 18.334 < 2e-16 ***
## grp11 60.481 6.947 8.706 8.32e-14 ***
## grp21 58.508 6.947 8.422 3.38e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 34.74 on 97 degrees of freedom
## Multiple R-squared: 0.602, Adjusted R-squared: 0.5938
## F-statistic: 73.35 on 2 and 97 DF, p-value: < 2.2e-16Αυτό σημαίνει οτι θέλουμε να εξηγήσουμε την έκφραση ως συνάρτηση και των δύο παραγόντων, grp1, και grp2. Όμως θεωρούμε ότι η επίδραση του grp2 ειναι η ΙΔΙΑ και για τα επίπεδα του grp1 και αντίστροφα. Αυτό θα γίνει κατανοητό με το παρακάτω διάγραμμα
## note: the function uname removes the names from a vector
{plot(as.numeric(grp1)+rn-1, simdata[i,], col=pointcolors)
abline(a=unname(m2$coefficients[1]), b=m2$coefficients[2], col="lightblue", lw=2)
abline(a=m2$coefficients[1]+m2$coefficients[3], b=m2$coefficients[2], col="pink", lw=2)
legend("top", legend=c("grp1, grp2 == 0,0", "grp1,grp2 == 0,1", "grp1,grp2 == 1,0", "grp1,grp2 == 1,1"), col=c("red", "black", "blue", "green"), pch=19)
}
Τέλος στο πιο πολύπλοκο γραμμικό μοντέλο επιτρέπουμε και τις αλληλεπιδράσεις μεταξύ των παραγόντων. Αυτό σημαίνει ότι μπορεί η επίδραση του grp2 στην έκφραση να ειναι διαφορετική αν το grp1 == 1 ή το grp1 == 2.
#NOTE: the symbol * between grp1 and grp2
m3 = lm(simdata[i,]~grp1*grp2)
m3$coefficients## (Intercept) grp11 grp21 grp11:grp21
## 116.06294 48.97200 46.99879 23.01853summary(m3)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 116.06294 6.884015 16.859774 1.892171e-30
## grp11 48.97200 9.735467 5.030267 2.278489e-06
## grp21 46.99879 9.735467 4.827584 5.210780e-06
## grp11:grp21 23.01853 13.768030 1.671883 9.780393e-02Αυτό σημαίνει οτι θέλουμε να εξηγήσουμε την έκφραση ως συνάρτηση και των δύο παραγόντων, grp1, και grp2. Όμως θεωρούμε ότι η επίδραση του grp2 ειναι η ΜΠΟΡΕΙ ΝΑ ΜΗΝ ΕΙΝΑΙ Η ΙΔΙΑ για τα επίπεδα του grp1 και αντίστροφα. Αυτό θα γίνει κατανοητό με το παρακάτω διάγραμμα
## note: the function uname removes the names from a vector
{plot(as.numeric(grp1)+rn-1, simdata[i,], col=pointcolors)
abline(a=unname(m3$coefficients[1]), b=m3$coefficients[2], col="lightblue", lw=2)
abline(a=m3$coefficients[1]+m3$coefficients[3], b=m3$coefficients[2]+m3$coefficients[4], col="pink", lw=2)
legend("top", legend=c("grp1, grp2 == 0,0", "grp1,grp2 == 0,1", "grp1,grp2 == 1,0", "grp1,grp2 == 1,1"), col=c("red", "black", "blue", "green"), pch=19)
}
Ουσιαστικά εδώ η αλλαγή δεν φαίνεται, επειδή όπως έχουν κατασκευαστεί τα δεδομένα δεν υπάρχει στην ουσία αλληλεπίδραση μεταξύ των grp1 και grp2. Η επίδραση του grp2 στην έκφραση ειναι ανεξάρτητη της τιμής του grp1.
Ας φτιάξουμε ένα άλλο σύνολο δεδομένων, όπου θα υπάρχει αλληλεπίδραση μεεταξύ των δύο παραγόντων.
## Data generation
nrow=10
ncol=100
grp1 = factor(rep(c(0,1), each=ncol/2))
grp2 = factor(rep(c(0,1), ncol/2))
simdata = matrix(rnorm(nrow*ncol, 100, 30), nrow=nrow, ncol=ncol)
simdata[,((ncol/2+1):ncol)] = matrix(rnorm(nrow*ncol/2, 150, 30), nrow=nrow, ncol=ncol/2)
simdata[, seq(from=2, to=ncol/2, by=2)] = simdata[, seq(from=2, to=ncol/2, by=2)] + matrix(rnorm(nrow*ncol/4, -30, 30), nrow=nrow, ncol=ncol/4)
simdata[, seq(from=ncol/2+2, to=ncol, by=2)] = simdata[, seq(from=ncol/2+2, to=ncol, by=2)] + matrix(rnorm(nrow*ncol/4, 70, 30), nrow=nrow, ncol=ncol/4)
## Visualization with a heatmap
heatmap.2(simdata, Colv=F, Rowv=F, dendrogram='none', trace = 'none')
#NOTE: the symbol * between grp1 and grp2
m3 = lm(simdata[i,]~grp1*grp2)
m3$coefficients## (Intercept) grp11 grp21 grp11:grp21
## 104.75430 48.29763 -40.36663 101.63603summary(m3)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 104.75430 7.795597 13.437624 8.756354e-24
## grp11 48.29763 11.024639 4.380881 3.017455e-05
## grp21 -40.36663 11.024639 -3.661493 4.102615e-04
## grp11:grp21 101.63603 15.591194 6.518810 3.287917e-09Αυτό σημαίνει οτι θέλουμε να εξηγήσουμε την έκφραση ως συνάρτηση και των δύο παραγόντων, grp1, και grp2. Όμως θεωρούμε ότι η επίδραση του grp2 ειναι η ΜΠΟΡΕΙ ΝΑ ΜΗΝ ΕΙΝΑΙ Η ΙΔΙΑ για τα επίπεδα του grp1 και αντίστροφα. Αυτό θα γίνει κατανοητό με το παρακάτω διάγραμμα
## note: the function uname removes the names from a vector
{plot(as.numeric(grp1)+rn-1, simdata[i,], col=pointcolors)
abline(a=unname(m3$coefficients[1]), b=m3$coefficients[2], col="lightblue", lw=2)
abline(a=m3$coefficients[1]+m3$coefficients[3], b=m3$coefficients[2]+m3$coefficients[4], col="pink", lw=2)
legend("top", legend=c("grp1, grp2 == 0,0", "grp1,grp2 == 0,1", "grp1,grp2 == 1,0", "grp1,grp2 == 1,1"), col=c("red", "black", "blue", "green"), pch=19)
}
i = 1
mean1.0 = mean(simdata[i,grp1.0])
mean1.1 = mean(simdata[i,grp1.1])
mean2.0 = mean(simdata[i,grp2.0])
mean2.1 = mean(simdata[i,grp2.1])
mean1.0_2.0 = mean(simdata[i,grp1.0&grp2.0])
mean1.0_2.1 = mean(simdata[i,grp1.0&grp2.1])
mean1.1_2.0 = mean(simdata[i,grp1.1&grp2.0])
mean1.1_2.1 = mean(simdata[i,grp1.1&grp2.1])
sim.means = c(m1.0 = mean1.0, m1.1=mean1.1, m2.0 = mean2.0, m2.1 = mean2.1,
m1.0_2.0=mean1.0_2.0,
m1.0_2.1=mean1.0_2.1,
m1.1_2.0=mean1.1_2.0,
m1.1_2.1=mean1.1_2.1)
sim.means## m1.0 m1.1 m2.0 m2.1 m1.0_2.0 m1.0_2.1 m1.1_2.0
## 84.57099 183.68663 128.90312 139.35450 104.75430 64.38767 153.05193
## m1.1_2.1
## 214.32133Προσοχή !!! Αν υπάρχει αλληλεπίδραση, τότε δεν έχει πολύ νόημα να μιλάμε για παράγοντες χωρίς να θεωρήσουμε την αλληλεπίδραση
#NOTE: the symbol + between grp1 and grp2
m4 = lm(simdata[i,]~grp1+grp2)
m4$coefficients## (Intercept) grp11 grp21
## 79.34530 99.11564 10.45138summary(m4)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 79.34530 8.066981 9.835811 3.033134e-16
## grp11 99.11564 9.314947 10.640495 5.582823e-18
## grp21 10.45138 9.314947 1.122001 2.646317e-01## note: the function uname removes the names from a vector
{plot(as.numeric(grp1)+rn-1, simdata[i,], col=pointcolors)
abline(a=unname(m4$coefficients[1]), b=m4$coefficients[2], col="lightblue", lw=2)
abline(a=m4$coefficients[1]+m4$coefficients[3], b=m4$coefficients[2], col="pink", lw=2)
legend("top", legend=c("grp1, grp2 == 0,0", "grp1,grp2 == 0,1", "grp1,grp2 == 1,0", "grp1,grp2 == 1,1"), col=c("red", "black", "blue", "green"), pch=19)
}
ΠΡΟΣΟΧΗ!!! Αυτό προφανώς ειναι λάθος. Και το λάθος δημιουργήθηκε επειδή υποθέσαμε οτι δεν υπάρχει αλληλεπίδραση, ενώ στην πραγματικότητα υπάρχει.*
Το ερώτημα μπορεί να διατυπωθεί πιο γενικά: Αν πολλαπλοί παράγοντες επιδρούν στην έκφραση γονιδίων τότε μπορούμε να αγνοήσουμε τους υπόλοιπους και να μελετήσουμε μόνον αυτό που μας αφορά;
Θα δούμε ότι αυτό δεν ειναι πολύ καλή ιδέα. Γενικά, ότι πληροφορία έχουμε για το δείγμα που αναλύουμε ειναι καλό να τις χρησιμοποιούμε στο μοντέλο της έκφρασης που έχουμε φτιάξει.
Μελέτη έκφρασης αγνοώντας κάποιον παράγοντα (δεν συστήνεται να το κάνετε αυτό)
Θεωρήστε την έκφραση του παρακάτω γονιδίου
## Η έκφραση ενός γονιδίου σαν "βάση" ειναι 200 για τα αρσενικά ατομα
males_g1 = rnorm(100, 200, 30)
## Η έκφραση ενός γονιδίου σαν "βάση" ειναι 250 για τα θηλυκά ατομα
females_g1 = rnorm(100, 250, 30)
## Ας υποθέσουμε οτι καπνίζουν τα "ζυγά" αρσενικά άτομα και τα "ζυγά θηλυκά άτομα"
## Επίσης ας υποθέσουμε οτι στα αρσενικά η επίδραση του καπνίσματος ειναι +50, ενώ στα θηλυκά -50. Τότε
zyga=seq(from=2, to=100, by=2)
allzyga = seq(from=2, to=200, by=2)
males_g1[zyga] = males_g1[zyga] + rnorm(length(zyga), 50, 10)
females_g1[zyga] = females_g1[zyga] + rnorm(length(zyga), -50, 10)
g1 = c(males_g1, females_g1)
# for plotting
col=rep(c("blue", "red"), each=100)
pch=rep(c(18, 1), 100)
xlabel = rep(c(0,1), 100) + rnorm(200, 0, 0.03)
plot(xlabel, g1, col=col, pch=pch)
## Differential gene expression without considering the GENDER
## t-test
## allzyga are the indexes for smokers
## -allzyga denotes the complementary indexes of allzyga, i.e., non-smokers
boxplot(list(nonsmokers=g1[-allzyga], smokers=g1[allzyga]))
Το πρόβλημα που προκύπτει ειναι οτι τελικά το t-test δεν θα βγει σημαντικό στο σύνολο των δειγμάτων, αν και ειναι σημαντικό τοσο για τα θηλυκά όσο και για τα αρσενικά άτομα.
mytest = t.test(g1[allzyga], g1[-allzyga])
print(mytest)##
## Welch Two Sample t-test
##
## data: g1[allzyga] and g1[-allzyga]
## t = -1.1404, df = 194.58, p-value = 0.2555
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -16.436222 4.392173
## sample estimates:
## mean of x mean of y
## 224.9757 230.9977Άρα συνοψίζοντας το κομμάτι αυτό: Αν υπάρχει κάποιος παράγοντας ακόμη και αν δεν μας ενδιαφέρει πρέπει να μπαίνει στο στατιστικό μοντέλο. Στην περίπτωση που δείξαμε παραπάνω, αυτό συμβαίνει επειδή υπάρχει αλληλεπίδραση του SMOKING με το GENDER (Αυτό φαίνεται επειδή στα αρσενικά άτομα το κάπνισμα αυξάνει την τιμή έκφραης, ενω στα θηλυκά την μειώνει). Όμως πρόβλημα μπορέι να προκύψει και όταν δεν υπάρχει τέτοια αλληλεπίδραση
## Η έκφραση ενός γονιδίου σαν "βάση" ειναι 200 για τα αρσενικά ατομα
males_g1 = rnorm(100, 200, 30)
## Η έκφραση ενός γονιδίου σαν "βάση" ειναι 250 για τα θηλυκά ατομα
females_g1 = rnorm(100, 500, 30)
## Ας υποθέσουμε οτι καπνίζουν τα "ζυγά" αρσενικά άτομα και τα "ζυγά θηλυκά άτομα"
## Επίσης ας υποθέσουμε οτι στα αρσενικά η επίδραση του καπνίσματος ειναι +50, ενώ στα θηλυκά -50. Τότε
zyga=seq(from=2, to=100, by=2)
allzyga = seq(from=2, to=200, by=2)
males_g1[zyga] = males_g1[zyga] + rnorm(length(zyga), 30, 10)
females_g1[zyga] = females_g1[zyga] + rnorm(length(zyga), 30, 10)
g1 = c(males_g1, females_g1)
col=rep(c("blue", "red"), each=100)
pch=rep(c(18, 1), 100)
xlabel = rep(c(0,1), 100) + rnorm(200, 0, 0.03)
plot(xlabel, g1, col=col, pch=pch)
#As a boxplot
boxplot(list(nonsmokers=g1[-allzyga], smokers=g1[allzyga]))
Για να πάροουμε το t-test
mytest = t.test(g1[allzyga], g1[-allzyga])
print(mytest)##
## Welch Two Sample t-test
##
## data: g1[allzyga] and g1[-allzyga]
## t = 1.1639, df = 197.78, p-value = 0.2459
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -17.38170 67.44722
## sample estimates:
## mean of x mean of y
## 374.1051 349.0723Ενώ το t-test στα αρσενικά άτομα είναι στατιστικά σημαντικό
zyga## [1] 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34
## [18] 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68
## [35] 70 72 74 76 78 80 82 84 86 88 90 92 94 96 98 100mytest = t.test(males_g1[zyga], males_g1[-zyga])
print(mytest)##
## Welch Two Sample t-test
##
## data: males_g1[zyga] and males_g1[-zyga]
## t = 5.4028, df = 96.792, p-value = 4.7e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 18.59492 40.19048
## sample estimates:
## mean of x mean of y
## 227.5618 198.1691Άρα το πρόβλημα που προκύπτει αν δεν λάβουμε υπόψιν τους επιπλέον παράγοντες του πειράματος ειναι οτι μπορεί να μην βρούμε διαφορική έκφραση ενω στην πραγματικότητα υπάρχει. Δεν έχω υπόψιν μου κάποια περίπτωση για το αντίθετο: δηλαδή να βρούμε διαφορική έκφραση στο σύνολο του δείγματος για έναν παράγοντα ενω στην πραγματικότητα σε κανένα από τα υποσύνολα που ορίζουν οι υπόλοιποι παράγοντες να μην υπάρχει διαφορική έκφραση.
Στο παραπάνω παράδειγμα, αν ειχαμε λάβει υπόψιν το γεγονός οτι υπάρχουν δύο φύλα (ακόμη και χωρίς αλληλεπίδραση)
gender = factor(rep(c(0,1), each=100))
smoking = factor(rep(c(0,1), 100))
lm.test = lm(g1~smoking+gender)
print(summary(lm.test))##
## Call:
## lm(formula = g1 ~ smoking + gender)
##
## Residuals:
## Min 1Q Median 3Q Max
## -65.762 -19.647 -4.203 18.382 78.405
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 200.349 3.446 58.137 <2e-16 ***
## smoking1 25.033 3.979 6.291 2e-09 ***
## gender1 297.446 3.979 74.749 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.14 on 197 degrees of freedom
## Multiple R-squared: 0.9662, Adjusted R-squared: 0.9658
## F-statistic: 2814 on 2 and 197 DF, p-value: < 2.2e-16Εδω βλέπουμε προφανώς, οτι “conditional on the gender”, δηλαδή δεδομένου ότι ξέρουμε το φύλο το κάπνισμα έχει επίδραση στην έκφραση γονιδίων. Αν δεν το ξέρουμε, τοτε ειναι τοσο μεγάλο το variance (επειδή το φύλο προσδίδει μεγάλο διαφορετικότητα στην έκφραση σε σχέση με το κάπνισμα) που δεν μπορούμε να βρούμε με το υπάρχον δείγμα οτι υπάρχει επίδραση του καπνίσματος.
Παράδειγμα 2 - Ευρεση γονιδίων που διαφοροποιούνται εξαιτίας του καπνίσματος
Το Dataset GDS3709 που χρησιμοποιούμε σαν παράδειγμα μελετάει δύο παράγοντες. Το κάπνισμα και το φύλο. Θα μπορούσαμε να μελετήσουμε την επίδραση του καπνίσματος αγνοώντας το φύλο. Όμως όπως είδαμε στο παραπάνω παράδειγμα, με τα προσομοιωμένα datasets αν υπάρχει και άλλος παράγοντας καλό ειναι να τον λάβουμε υπόψιν.
Ας υποθέσουμε λοιπόν ότι μας ενδιαφέρει ο παράγοντας κάπνισμα. Δεν θα θεωρήσουμε οτι υπάρχει αλληλεπίδραση μεταξύ φύλου και καπνίσματος. Το φύλο θα έχει το ρόλο μιας “αδιάφορης” για εμάς μεταβλητής που όμως πρέπει να “ελέγξουμε” κατά’τον στατιστικό έλεγχο.
## expr is the vector with expression values
## f1 and f2, are the two factors for the test, respectively
## f1 represents smoking and f2 represents gender
mytest = function(expr, f1, f2){
f1 = factor(f1)
f2 = factor(f2)
mylm = lm(expr~f1+f2)
#print(summary(mylm)$coefficients)
## first row is intercept, thus 2nd row is the row that we are intrested in (f1)
pvalue = summary(mylm)$coefficients[2,"Pr(>|t|)"]
return(pvalue)
}
## load the data
a <- read.table("GDS3709.soft.clean", h=TRUE)
#des <- read.table("GDS3709.desc")
##b <- a[, 3:(39+2)]
data = a[,3:ncol(a)]
## notice: the first factor is the basal factor
gender = factor(rep(c("F", "M"), c(39, 40)), levels = c("F", "M"))
smoking = factor( rep( rep( c("s", "n"), 2), c(19, 20, 20, 20)), levels=c("n", "s"))
mypvalues = apply(data, 1, mytest, smoking, gender)
length(mypvalues)## [1] 54675head(mypvalues)## [1] 0.18352009 0.05255009 0.52163062 0.15044807 0.31345358 0.03541250Επειδή έχουμε κάνει πολλά τεστ, πρέπει να διορθώσουμε για έλεγχο πολλαπλών υποθέσεων
## correct the p-values
p.corrected.bonf = p.adjust(mypvalues, method = "bonferroni")
p.corrected.fdr = p.adjust(mypvalues, method="fdr")
## the gene names are in the second column of the data frame
gene.names = a[,2]
genes.bonf = gene.names[p.corrected.bonf < 0.05]
genes.bonf## [1] AHCYL1 IMPDH2 CYP1B1 CYP1B1 CYP1B1 BAMBI SURF2
## [8] ALDH3A1 ITGB7 CYP1A1 S100A12 ITGAL CYTL1 DNLZ
## [15] MFSD12 CD207 CORO2A LINC00673 LGR6 JAML AHRR
## [22] CYSLTR1 TIFAB
## 30806 Levels: 1552563_a_at 1552829_at 1552867_at 1552961_at ... ZZZ3genes.fdr = gene.names[p.corrected.fdr < 0.05]
genes.fdr## [1] MAP3K6 UMODL1 LYNX1 FBXO7
## [5] MIER3 L3MBTL4 PER3 USP25
## [9] MS4A4A C10orf90 CASC22 UBE2G2
## [13] CHIC1 TMEM72 FCN1 LOC283482
## [17] AK097628 MAF LOC101927820 PER3
## [21] PLD1 YWHAZ PTP4A1 AHCYL1
## [25] AHCYL1 MIR1292 IPO7 PRDX2
## [29] EIF5B EIF5B TUFM CHD4
## [33] PSMD13 HSD17B4 NQO1 NQO1
## [37] PFDN1 SRM ANXA2 KIAA0100
## [41] ELAC2 CEACAM5 IMPDH2 LOC101928623
## [45] DKK3 SLC6A8 DDAH2 SLC35B1
## [49] CYP1B1 CYP1B1 CYP1B1 GSTM3
## [53] LYN SERTAD2 SNRPD1 ABHD14A-ACY1
## [57] R3HDM1 DNTTIP2 STK39 TCEB1
## [61] GPX2 CNPY2 DNAJB12 SOX9
## [65] GSTA4 LTBR CCDC86 RABEPK
## [69] KATNB1 ABCB6 LGALS9 FCGBP
## [73] NPRL2 BAMBI POLG PRAF2
## [77] GGH SKIV2L VSNL1 NQO2
## [81] IGF2BP3 IGF2BP3 PIK3CD MAPT
## [85] FZD6 PHACTR2 PHACTR2 TOE1
## [89] AKR1C1 RUNX3 RUNX3 MRC1
## [93] UGT1A1 CCL5 EPHB6 COL9A3
## [97] MIR8085 NCF4 ACAP1 SURF2
## [101] MPPED2 CORO2A PCSK5 PCSK5
## [105] ALDH3A1 ITGB7 STAC CYP1A1
## [109] ITGAM JAK2 S100A12 CASP1
## [113] UGT1A3 CEACAM7 CEACAM7 CA4
## [117] CFP MAPT UGT1A1 PTGES
## [121] PIR OGT NCF4 PPARG
## [125] UGT1A3 GORASP2 ATRX DDX18
## [129] HNRNPUL2-BSCL2 TXN2 PRPF4 POLR2H
## [133] MKNK1 POP7 BMP7 BMP7
## [137] CYB561D2 CD24 LRRN3 CDC42EP2
## [141] ABCB1 CST7 CHRNA3 CD1A
## [145] MTX1 ANXA2 ABLIM1 EIF3K
## [149] NQO1 AKR1C4 TTC39A NPRL3
## [153] ITPK1 GSTM4 CSRP2 BMP7
## [157] LEPROT LEPROT LEPROT SLC4A4
## [161] EIF3B VEGFA LOC101930400 FCER1A
## [165] NOLC1 RSL1D1 SERPINE2 OGT
## [169] CHMP7 SULF1 TCF4 DNAJC8
## [173] CAST SLC25A44 EIF3K CBS
## [177] ARHGAP45 MCF2L PTPN2 PCBP2
## [181] ATG2A POM121C MFSD9 ITGAL
## [185] ANXA2 RPL14 PCSK5 SECTM1
## [189] MET CD6 GGCX MAP4K1
## [193] DKK3 NOP2 SPRR1A ZNF394
## [197] ZNF75D DDAH2 UGT1A3 DDAH2
## [201] DDX27 CYFIP2 IGF2BP3 MS4A1
## [205] GALNT2 ARHGAP5 MSRB1 C9orf78
## [209] NOL6 TOMM6 IPO4 TIMM9
## [213] TMEM14A THAP7 ERO1A TEX264
## [217] NIT2 MRPS33 LRRC8D GRAMD3
## [221] CCDC51 RNF25 LAGE3 TTI2
## [225] ADGRL1 TRIM62 CEP76 OTUB2
## [229] WFDC1 CLIC3 FBXO22 GDPD3
## [233] SLC27A5 PRR7 RASAL1 SLC13A4
## [237] CYTL1 LYVE1 DNLZ MFSD12
## [241] CFAP45 USP25 CD207 UGT1A3
## [245] EIF3K LSG1 LSG1 PPAN-P2RY11
## [249] ASB6 CLEC7A FAAP100 SDHAF1
## [253] ARRB1 MECOM WLS GSTA3
## [257] TBC1D2 ACE2 PNPO ARRB1
## [261] ZNF703 NOB1 TMEM14B TMEM14C
## [265] SLC25A28 MAD2L2 SLC4A11 OXR1
## [269] PINX1 GDA AF116671 UBQLN4
## [273] CYB561A3 DNTTIP1 GPT2 MRPS25
## [277] KRTCAP2 CPSF3 MRPS26 MAPT
## [281] TSPAN5 PET117 ERO1A MTMR10
## [285] LOC101930349 CCDC137 NLN ALKBH6
## [289] PLEKHG2 MANEAL MCEE PTPN14
## [293] MRPS15 MECOM CYGB CORO2A
## [297] FNDC10 KREMEN1 TMEM119 LINC00673
## [301] ACE PIK3R5 C10orf99 C10orf99
## [305] LGR6 LAYN JAML NADK2
## [309] SORCS2 NEDD8-MDP1 NEDD8-MDP1 GJC1
## [313] POLR2J2 SNX20 WLS ADAMTS15
## [317] CCT5 RASSF6 ADRB1 AHRR
## [321] FAM118B ARHGAP42 MBNL3 KDM5A
## [325] FUOM LINC00673 CYSLTR1 YIF1B
## [329] NMRAL1 LOC102724842 CYSLTR1 USF1
## [333] PRLHR OIP5-AS1 TUBD1 ANKRD33B
## [337] PIGN HSPA9 LOC101927254 LTBR
## [341] TEX264 LOC100129884 234954_at SAMHD1
## [345] CSRNP3 CYYR1 ZNF542P CCDC117
## [349] PELI3 VPS37D RASSF6 BF115793
## [353] HS1BP3 TIFAB BF938956 KIAA1217
## [357] TFRC BF056769 SYTL3 BMPR2
## [361] OSBPL6 TNFRSF11A RHOJ KCNA2
## [365] TIFAB AA588092 T26531 TLCD2
## [369] NMRAL1P1 GPR82 MYO1G TSSK2
## [373] SLC25A44 CHST3 NEFH KCTD2
## [377] MADD PRDX2 DDX23 ARRB1
## [381] GAPDH
## 30806 Levels: 1552563_a_at 1552829_at 1552867_at 1552961_at ... ZZZ3Παρόμοια ανάλυση χρησιμοποιώντας το limma
Η παραπάνω ανάλυση μπορεί να αυτοματοποιηθεί χρησιμοποιώντας το limma πακέτο, το model.matrix και τα contrasts
library(limma)
## load the data
a <- read.table("GDS3709.soft.clean", h=TRUE)
#des <- read.table("GDS3709.desc")
##b <- a[, 3:(39+2)]
data = a[,3:ncol(a)]
## notice: the first factor is the basal factor
gender = factor(rep(c("F", "M"), c(39, 40)), levels = c("F", "M"))
smoking = factor( rep( rep( c("s", "n"), 2), c(19, 20, 20, 20)), levels=c("n", "s"))
## I use the 0 in order to show both of the levels of smoking as reference
modmat = model.matrix(~0+smoking+gender)
head(modmat)## smokingn smokings genderM
## 1 0 1 0
## 2 0 1 0
## 3 0 1 0
## 4 0 1 0
## 5 0 1 0
## 6 0 1 0colnames(modmat) = c("nosmoking", "smoking", "males")
modmat## nosmoking smoking males
## 1 0 1 0
## 2 0 1 0
## 3 0 1 0
## 4 0 1 0
## 5 0 1 0
## 6 0 1 0
## 7 0 1 0
## 8 0 1 0
## 9 0 1 0
## 10 0 1 0
## 11 0 1 0
## 12 0 1 0
## 13 0 1 0
## 14 0 1 0
## 15 0 1 0
## 16 0 1 0
## 17 0 1 0
## 18 0 1 0
## 19 0 1 0
## 20 1 0 0
## 21 1 0 0
## 22 1 0 0
## 23 1 0 0
## 24 1 0 0
## 25 1 0 0
## 26 1 0 0
## 27 1 0 0
## 28 1 0 0
## 29 1 0 0
## 30 1 0 0
## 31 1 0 0
## 32 1 0 0
## 33 1 0 0
## 34 1 0 0
## 35 1 0 0
## 36 1 0 0
## 37 1 0 0
## 38 1 0 0
## 39 1 0 0
## 40 0 1 1
## 41 0 1 1
## 42 0 1 1
## 43 0 1 1
## 44 0 1 1
## 45 0 1 1
## 46 0 1 1
## 47 0 1 1
## 48 0 1 1
## 49 0 1 1
## 50 0 1 1
## 51 0 1 1
## 52 0 1 1
## 53 0 1 1
## 54 0 1 1
## 55 0 1 1
## 56 0 1 1
## 57 0 1 1
## 58 0 1 1
## 59 0 1 1
## 60 1 0 1
## 61 1 0 1
## 62 1 0 1
## 63 1 0 1
## 64 1 0 1
## 65 1 0 1
## 66 1 0 1
## 67 1 0 1
## 68 1 0 1
## 69 1 0 1
## 70 1 0 1
## 71 1 0 1
## 72 1 0 1
## 73 1 0 1
## 74 1 0 1
## 75 1 0 1
## 76 1 0 1
## 77 1 0 1
## 78 1 0 1
## 79 1 0 1
## attr(,"assign")
## [1] 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$smoking
## [1] "contr.treatment"
##
## attr(,"contrasts")$gender
## [1] "contr.treatment"myfit = lmFit(object=data, design = modmat )
mycontr = makeContrasts(smoking-nosmoking, levels = modmat)
fitc = contrasts.fit(myfit, mycontr)
fitc = eBayes(fitc)
difTable = topTable(fitc)
difTable## logFC AveExpr t P.Value adj.P.Val B
## 36708 91.20237 138.52975 8.600944 7.318853e-13 4.001583e-08 17.869145
## 38609 36.71535 87.34291 7.780183 2.769455e-11 7.530571e-07 14.665652
## 29122 24.85955 49.10924 7.689495 4.132001e-11 7.530571e-07 14.311334
## 11885 222.41465 211.12380 7.155318 4.312740e-10 5.658587e-06 12.229046
## 37349 31.87156 86.03557 7.113581 5.174748e-10 5.658587e-06 12.066920
## 11886 80.83565 65.19885 6.695668 3.173051e-09 2.891442e-05 10.451036
## 15197 42.40659 49.40544 6.413110 1.066616e-08 8.331031e-05 9.368711
## 15071 1616.27922 3660.62025 6.273696 1.930472e-08 1.272879e-04 8.838606
## 15311 264.98513 894.35316 6.254379 2.095274e-08 1.272879e-04 8.765388
## 29713 98.50650 165.20063 6.069341 4.575549e-08 2.499695e-04 8.067091rownames(difTable)[difTable$adj.P.Val < 0.05]## [1] "36708" "38609" "29122" "11885" "37349" "11886" "15197" "15071"
## [9] "15311" "29713"sigIds = as.numeric(rownames(difTable)[difTable$adj.P.Val < 0.05])
sigGenes = gene.names[sigIds]
sigGenes## [1] LINC00673 AHRR CYTL1 CYP1B1 JAML CYP1B1 CYP1A1
## [8] ALDH3A1 S100A12 CD207
## 30806 Levels: 1552563_a_at 1552829_at 1552867_at 1552961_at ... ZZZ3genes.fdr## [1] MAP3K6 UMODL1 LYNX1 FBXO7
## [5] MIER3 L3MBTL4 PER3 USP25
## [9] MS4A4A C10orf90 CASC22 UBE2G2
## [13] CHIC1 TMEM72 FCN1 LOC283482
## [17] AK097628 MAF LOC101927820 PER3
## [21] PLD1 YWHAZ PTP4A1 AHCYL1
## [25] AHCYL1 MIR1292 IPO7 PRDX2
## [29] EIF5B EIF5B TUFM CHD4
## [33] PSMD13 HSD17B4 NQO1 NQO1
## [37] PFDN1 SRM ANXA2 KIAA0100
## [41] ELAC2 CEACAM5 IMPDH2 LOC101928623
## [45] DKK3 SLC6A8 DDAH2 SLC35B1
## [49] CYP1B1 CYP1B1 CYP1B1 GSTM3
## [53] LYN SERTAD2 SNRPD1 ABHD14A-ACY1
## [57] R3HDM1 DNTTIP2 STK39 TCEB1
## [61] GPX2 CNPY2 DNAJB12 SOX9
## [65] GSTA4 LTBR CCDC86 RABEPK
## [69] KATNB1 ABCB6 LGALS9 FCGBP
## [73] NPRL2 BAMBI POLG PRAF2
## [77] GGH SKIV2L VSNL1 NQO2
## [81] IGF2BP3 IGF2BP3 PIK3CD MAPT
## [85] FZD6 PHACTR2 PHACTR2 TOE1
## [89] AKR1C1 RUNX3 RUNX3 MRC1
## [93] UGT1A1 CCL5 EPHB6 COL9A3
## [97] MIR8085 NCF4 ACAP1 SURF2
## [101] MPPED2 CORO2A PCSK5 PCSK5
## [105] ALDH3A1 ITGB7 STAC CYP1A1
## [109] ITGAM JAK2 S100A12 CASP1
## [113] UGT1A3 CEACAM7 CEACAM7 CA4
## [117] CFP MAPT UGT1A1 PTGES
## [121] PIR OGT NCF4 PPARG
## [125] UGT1A3 GORASP2 ATRX DDX18
## [129] HNRNPUL2-BSCL2 TXN2 PRPF4 POLR2H
## [133] MKNK1 POP7 BMP7 BMP7
## [137] CYB561D2 CD24 LRRN3 CDC42EP2
## [141] ABCB1 CST7 CHRNA3 CD1A
## [145] MTX1 ANXA2 ABLIM1 EIF3K
## [149] NQO1 AKR1C4 TTC39A NPRL3
## [153] ITPK1 GSTM4 CSRP2 BMP7
## [157] LEPROT LEPROT LEPROT SLC4A4
## [161] EIF3B VEGFA LOC101930400 FCER1A
## [165] NOLC1 RSL1D1 SERPINE2 OGT
## [169] CHMP7 SULF1 TCF4 DNAJC8
## [173] CAST SLC25A44 EIF3K CBS
## [177] ARHGAP45 MCF2L PTPN2 PCBP2
## [181] ATG2A POM121C MFSD9 ITGAL
## [185] ANXA2 RPL14 PCSK5 SECTM1
## [189] MET CD6 GGCX MAP4K1
## [193] DKK3 NOP2 SPRR1A ZNF394
## [197] ZNF75D DDAH2 UGT1A3 DDAH2
## [201] DDX27 CYFIP2 IGF2BP3 MS4A1
## [205] GALNT2 ARHGAP5 MSRB1 C9orf78
## [209] NOL6 TOMM6 IPO4 TIMM9
## [213] TMEM14A THAP7 ERO1A TEX264
## [217] NIT2 MRPS33 LRRC8D GRAMD3
## [221] CCDC51 RNF25 LAGE3 TTI2
## [225] ADGRL1 TRIM62 CEP76 OTUB2
## [229] WFDC1 CLIC3 FBXO22 GDPD3
## [233] SLC27A5 PRR7 RASAL1 SLC13A4
## [237] CYTL1 LYVE1 DNLZ MFSD12
## [241] CFAP45 USP25 CD207 UGT1A3
## [245] EIF3K LSG1 LSG1 PPAN-P2RY11
## [249] ASB6 CLEC7A FAAP100 SDHAF1
## [253] ARRB1 MECOM WLS GSTA3
## [257] TBC1D2 ACE2 PNPO ARRB1
## [261] ZNF703 NOB1 TMEM14B TMEM14C
## [265] SLC25A28 MAD2L2 SLC4A11 OXR1
## [269] PINX1 GDA AF116671 UBQLN4
## [273] CYB561A3 DNTTIP1 GPT2 MRPS25
## [277] KRTCAP2 CPSF3 MRPS26 MAPT
## [281] TSPAN5 PET117 ERO1A MTMR10
## [285] LOC101930349 CCDC137 NLN ALKBH6
## [289] PLEKHG2 MANEAL MCEE PTPN14
## [293] MRPS15 MECOM CYGB CORO2A
## [297] FNDC10 KREMEN1 TMEM119 LINC00673
## [301] ACE PIK3R5 C10orf99 C10orf99
## [305] LGR6 LAYN JAML NADK2
## [309] SORCS2 NEDD8-MDP1 NEDD8-MDP1 GJC1
## [313] POLR2J2 SNX20 WLS ADAMTS15
## [317] CCT5 RASSF6 ADRB1 AHRR
## [321] FAM118B ARHGAP42 MBNL3 KDM5A
## [325] FUOM LINC00673 CYSLTR1 YIF1B
## [329] NMRAL1 LOC102724842 CYSLTR1 USF1
## [333] PRLHR OIP5-AS1 TUBD1 ANKRD33B
## [337] PIGN HSPA9 LOC101927254 LTBR
## [341] TEX264 LOC100129884 234954_at SAMHD1
## [345] CSRNP3 CYYR1 ZNF542P CCDC117
## [349] PELI3 VPS37D RASSF6 BF115793
## [353] HS1BP3 TIFAB BF938956 KIAA1217
## [357] TFRC BF056769 SYTL3 BMPR2
## [361] OSBPL6 TNFRSF11A RHOJ KCNA2
## [365] TIFAB AA588092 T26531 TLCD2
## [369] NMRAL1P1 GPR82 MYO1G TSSK2
## [373] SLC25A44 CHST3 NEFH KCTD2
## [377] MADD PRDX2 DDX23 ARRB1
## [381] GAPDH
## 30806 Levels: 1552563_a_at 1552829_at 1552867_at 1552961_at ... ZZZ3Υπάρχουν μερικές διαφορές, όμως ειναι μικρές και ειναι εξαιτίας του λίγο διαφορετικού τρόπου που το lm χειρίζεται το variance σε σχέση με το limma. Το limma χρησιμοποιεί το eBayes ώστε να έχει μια καλύτερη εκτίμηση του variance σε κάθε γονίδιο.
Σε κάθε περίπτωση το ερώτημα ειναι:
Ωραία! βρήκαμε τα γονίδια που ειναι διαφοροποιημένα. Τι κάνουν τα γονίδια αυτά;
Για να απαντήσουμε το ερώτημα αυτό, θα χρησιμοποιήσουμε μια βιβλιοθήκη που ονομάζεται gprofiler. Το ίδιο εργαλείο υπάρχει και σε online version https://biit.cs.ut.ee/gprofiler/gost
gprofiler ή αλλιώς Gene set enrichment analysis ή GO analysis
Η ιδέα ειναι απλή. Θέλουμε να δούμε αν τα γονίδια που βρίσκουμε σαν αποτέλεσμα σχετίζονται με κάποια βιολογική διεργασία περισσότερο από οτι θα περιμέναμε στην τύχη. Έτσι, Αν σε μια βιολογική κατηγορία (πχ γονιδια που σχετίζονται με την καρδιά) υπάρχουν γενικά Χ γονίδια (πχ 1000) από τα συνολικά 10000 γονίδια του ανθρώπου (δηλαδή, έστω οτι με την καρδιά σχετίζονται περίπου το 10% των γονιδίων μας), και στην δική μας λίστα γονιδίων που ειναι διαφοροποιημένα βρίσκω 6 στα 10 οτι σχετίζονται με την καρδιά, τότε προφανώς υπάρχει μια υπεραντιπροσώπευση της λειτουργίας “σχετίζεται με την καρδιά” στην δική μου λίστα σε σχέση με το σύνολο των γονιδίων. Τέτοιες υπερ-αντιπροσωπευσεις ψάχνουμε να βρούμε.
##install.packages('gprofiler2')
library(gprofiler2)
gostres <- gost(query = sigGenes, organism = "hsapiens")
## To get the results in a table format
gostres$result[1,]## query significant p_value term_size query_size intersection_size
## 1 query_1 TRUE 7.493605e-05 122 7 4
## precision recall term_id source term_name
## 1 0.5714286 0.03278689 GO:0006805 GO:BP xenobiotic metabolic process
## effective_domain_size source_order parents
## 1 18117 2774 GO:0044237, GO:0071466## and plot a summary of the results
p <- gostplot(gostres, capped = FALSE, interactive = FALSE)
p
phyper(q=0, m = 122, n = 18117-122, k=7, lower.tail = FALSE)## [1] 0.04620388Η στατιστική ανάλυση πίσω από το gProfiler
Αν σκεφτούμε λίγο το gProfiler, τότε αρκεί να σκεφτούμε τα εξής:
Σύνολο W: Τα γονίδια που αντιστοιχούν σε μια βιολογική κατηγορία Χ (υποσύνολο από το Τ) Σύνολο Q: Τα γονίδια που παρέχουμε σαν είσοδο (πχ τα Differentially Expressed Genes) Σύνολο T: Το σύνολο των γονιδίων του οργανισμού Σύνολο S: To Υποσύνολο των γονιδίων του Q που ανήκουν στο Χ
Αυτό που εξετάζουμε ειναι αν το Q έχει περισσότερα Χ γονίδια από ότι θα περιμέναμε αν απλά το Q το παίρναμε τυχαία από το T. Δηλαδή για παράδειγμα, αν στο Τ υπάρχουν 10% γονίδια που αντιστοιχούν στην βιολογική κατήγορία Χ, ενώ στο Q το αντίστοιχο ποσοστό ειναι 40% τοτε αυτό σημαίνει ότι το Q δεν ειναι τυχαία κατηγορία από το Q αλλά “έχουμε επιλέξει” (ή το πείραμα δίνει τέτοια αποτελέσματα) γονίδια που έχουν περισσότερο από το X. Επομένως μπορούμε να πούμε ότι υπάρχει υπερ-αντιπροσώπευση του Χ στο Q σε σύγκριση με αυτό που θα περιμέναμε στην τύχη.
Για να εκτιμήσουμε την πιθανότητα να βρούμε τόσα στοιχεία του Q που ανήκουν στο Χ (δηλαδή το S), από αυτά που θα περιμέναμε αν το Q ήταν “απλά ένα τ υχαίο υποσύνολο το T”, θα χρησιμοποιήσουμε είτε την διωνυμική κατανομή είτε την υπεργεωμετρική κατανομή. Ειναι σχεδόν οι ίδιες, μόνο που η διωνυμική έπιτρέπει επανάθεση, ενώ η υπεργεωμετρική όχι. Η υπεργεωμετρική ειναι μάλλον πιο ακριβής, όμως και οι δύο κατανομές θα δώσουν παρόμοιο αποτέλεσμα τουλάχιστον όταν αναφερόμαστε σε βιολογικές κατηγορίες που έχουν αρκετό μεγάλο πλήθος γονιδίων.
Αν γράφαμε στην R τον κώδικα για την εκτίμηση της πιθανότητας να υπάρχει υπεραντιπροσώπευση της κατηγορίας Χ, τότε ο κώδικας θα ήταν ως εξής:
phyper(q = |S|-1, m = |W|, n = |T|-|W|, k = |Q|, lower.tail = F)
όπου, |S| είναι το μέγεθος του συνόλου S κοκ.