[Δείτε ένα ενημερωτικό βίντεο για την άσκηση εδώ]
Στη σημερινή άσκηση θα δούμε πώς αναλύουμε δεδομένα γονιδιακής έκφρασης που προέρχονται από ένα πείραμα μεγάλης κλίμακας. Σκοπός της άσκησης είναι να δούμε:
- Τη δομή δεδομένων γονιδιακής έκφρασης
- Τον υπολογισμό σχετικής έκφρασης και την στατιστική της αξιολόγηση
- Τρόπους αναπαράστασης των τιμών έκφρασης
- Εντοπισμό απορρυθμισμένων γονιδίων
Θα μελετήσουμε ένα πείραμα γονιδιακής έκφρασης σε μεγάλη κλίμακα και θα κάνουμε τους υπολογισμούς αλλά και τις γραφικές παραστάσεις με τη βοήθεια της R.
Δομή Δεδομένων γονιδιακής έκφρασης
Τα δεδομένα που θα μελετήσουμε προέρχονται από ένα πείραμα έκφρασης με τη χρήση μικροσυστοιχειών DNA (DNA microarrays). H ποσότητα του mRNA για 18703 γονίδια του ποντικού μελετήθηκε σε ένα στέλεχος αγρίου τύπου (Wt-control) σε σχέση με 5 διαφορετικές συνθήκες (Α, B, C, D, E). Kάθε πείραμα έγινε σε 4 βιολογικές επαναλήψεις. Το αποτέλεσμα είναι ένας πίνακας 18703x24 (τέσσερις επαναλήψεις για 5+1 συνθήκες). Οι κανονικοποιημένες τιμές έκφρασης περιέχονται στο αρχείο Expression_Data.tsv που θα βρείτε σε αυτό το link.
Ανάγνωση αρχείου στην R και επισκόπηση των δεδομένων
Θα διαβάσουμε το αρχείο dataset.tsv στην R (αφού το έχουμε αποθηκεύσει στο σωστό φάκελο και έχουμε θέσει αυτόν τον φάκελο, ως φάκελο εργασίας στην R κατά τα γνωστά) με την παρακάτω εντολή:
data<-read.delim("Expression_Data.tsv", header=T, sep="\t")που θα διαβάσει το αρχείο σε εναν πίνακα τύπου data.frame κρατώντας την πρώτη γραμμή για επικεφαλίδα
Για να πάρουμε μια ιδέα για τα περιεχόμενα του αρχείου μπορούμε απλώς να γράψουμε
head(data)οπότε και θα δούμε κάτι σαν το παρακάτω
Gene Wt_1 Wt_2 Wt_3 A_1 A_2 A_3 B_1 B_2 B_3 C_1 C_2 C_3 D_1 D_2 D_3 E_1 E_2 E_31 A1bg 3.9709 4.1174 4.0083 4.0271 4.2455 4.2449 3.7624 3.9783 3.9229 4.0293 3.8934 3.8316 3.8740 3.8785 3.8822 3.9112 3.9225 3.89132 A1cf 4.4204 4.4006 4.3663 4.2209 4.2776 4.2259 4.2970 4.3995 4.2479 4.3587 4.5550 4.4794 4.3744 4.4870 4.4653 4.4678 4.4799 4.50093 A2ld1 5.6756 5.9968 6.0226 5.4885 5.4918 5.3493 5.9686 6.1146 5.7160 7.7222 7.9057 7.7276 7.5761 7.7550 7.8009 7.6316 7.6726 7.61634 A2m 5.6274 6.8456 5.5570 4.9007 5.1312 5.5225 6.1936 5.7056 5.6869 5.7715 5.9767 5.6643 5.6800 5.3604 5.5496 5.1973 6.6096 6.01625 A3galt2 4.8495 4.8698 4.8714 4.9237 4.9776 5.2140 4.8181 4.8313 4.9395 4.9265 4.9678 4.8682 4.7075 4.8721 4.7357 4.7849 4.8205 4.79166 A4galt 8.1425 8.0071 8.2216 7.6556 7.5956 7.4403 7.9922 7.8596 7.9373 8.0794 8.2503 8.1693 8.2145 8.1684 8.1882 8.3621 8.4473 8.5088Eνώ αν θέλουμε να δούμε πιο αναλυτικά τι περιέχεται σε κάθε στήλη του data.frame μπορούμε να γράψουμε
str(data)που θα μας επιστρέψει αναλυτικά τα ονόματα κάθε στήλης και το είδος των δεδομένων που περιέχει η καθεμία.
Όπως συζητήσαμε και στο μάθημα, το πρώτο στάδιο της επεξεργασίας των δεδομένων είναι πάντα η κανονικοποίησή τους, η οποία γίνεται για να απαλειφθούν συστηματικές παρεκκλίσεις στις τιμές έκφρασης λόγω διαφορών στους χειρισμούς μεταξύ του κάθε δείγματος (απομόνωση περισσότερου υλικού (RNA), διαφορετικές ρυθμίσεις του μηχανήματος, τεχνικά σφάλματα κλπ). Στην περίπτωση του αρχείου μας αυτό έχει γίνει ήδη οπότε το βήμα αυτό θα παραλειφθεί. Για να βεβαιωθείτε, ωστόσο ότι οι τιμές είναι κανονικοποιημένες μπορείτε να δείτε τις κατανομές τους σε όλα τα δείγματα με τη χρήση απλών θηκογραμμάτων. Δοκιμάστε την παρακάτω εντολή:
boxplot(data[,2:19])H εντολή αυτή ζητά από την R να δημιουργήσει θηκογράμματα (boxplot) για κάθε στήλη από την 2 έως την 19 στο σετ δεδομένων data (αυτό λέει η σύνταξη data[,2:19])
Παρατήρηση των θηκογραμμάτων δείχνει πως σε κάθε δείγμα οι τιμές καταλαμβάνουν πολύ παραπλήσια εύρη, κάτι που σημαίνει πως δεν χρειάζεται κάποια περαιτέρω κανονικοποίηση.
Υπολογισμός Σχετικής Έκφρασης γονιδίων
Το επόμενο λογικό βήμα είναι να υπολογίσουμε τη σχετική έκφραση γονιδίων. Για να το κάνουμε αυτό θα πρέπει να σκεφτούμε:
α) έναντι ποιας συνθήκης θα υπολογίσουμε τη σχετική έκφραση, ποια θα είναι δηλαδή η συνθήκη/κατάσταση ελέγχου
β) με ποια ποσότητα θα ποσοτικοποιήσουμε τη σχετική έκφραση.
Στο παράδειγμά μας έχουμε μια κατάσταση που αφορά τα ζώα αγρίου τύπου χωρίς κανέναν επιπλέον χειρισμό. Είναι προφανές ότι θα επιλέξουμε αυτή τη συνθήκη (Wt) ως συνθήκη υποβάθρου. Πώς όμως θα υπολογίσουμε τη σχετική συχνότητα έκφρασης;
Όπως είδαμε στο μάθημα, ο πιο τυπικός τρόπος είναι η χρήση του λογαρίθμου λόγων έκφρασης (log-fold-change). Έτσι, για κάθε μία από τις άλλες συνθήκες Α, Β, C, D και Ε, θα υπολογίσουμε το δυαδικό λογάριθμο του λόγου της έκφρασης στη συνθήκη (Α, Β κλπ) προς την έκφραση στο Wt. H ποσότητα θα είναι:
log2(Έκφραση[Α]/Έκφραση[Wt])Στην περίπτωσή μας έχουμε 3 τιμές για την έκφραση σε κάθε κατάσταση οπότε θα πρέπει να υπολογίσουμε τη μέση τιμή της έκφρασης για κάθε γονίδιο στις δύο καταστάσεις. Αυτό θα το κάνουμε με τη χρήση μιας συνάρτησης της R που λέγεται rowMeans και που κάνει αυτό που λέει το όνομά της, υπολογίζει δηλαδή τη μέση τιμή σε κάθε σειρά ενός πίνακα. Αρκεί λοιπόν να επιλέξουμε τις σωστές στήλες για να υπολογίσουμε:
wt<-rowMeans(data[,2:4])τη μέση έκφραση στο Wt (στήλες 2 έως 4 του πίνακα data) και αντίστοιχα:
A<-rowMeans(data[,5:7])για την κατάσταση Α
Η σχετική έκφραση τώρα θα είναι:
logFC_A<-A-wtΣτην προκειμένη περίπτωση εφαρμόζουμε αφαίρεση, καθώς οι τιμές μας είναι ήδη λογαριθμημένες.
Η μεταβλητή Α_wt είναι τώρα μια λίστα με 18703 τιμές κάθεμιά από τις οποίες αντιστοιχεί στη σχετική έκφραση του κάθε γονιδίου του αρχικού σετ δεδομένων μεταξύ της κατάστασης Α και του στελέχους αγρίου τύπου (κατάσταση υποβάθρου). Με απολύτως αντίστοιχο τρόπο να υπολογίσετε τις ποσότητες logFC_B, logFC_C, logFC_D, logFC_E.
Στατιστική σημασία των δεδομένων
O υπολογισμός της σχετικής έκφρασης είναι το ένα μόνο μέρος της ανάλυσης. Το δεύτερο, εξίσου σημαντικό αν όχι σημαντικότερο, είναι η αξιολόγηση του κατά πόσο η μεταβολή της έκφρασης μεταξύ δύο καταστάσεων είναι στατιστικά σημαντικη, αντανακλά δηλαδή ένα φυσικό γεγονός και δεν είναι αποτέλεσμα μιας τυχαίας διακύμανσης. Στο παράδειγμά μας είμαστε τυχεροί να έχουμε 4 πειραματικές επαναλήψεις για κάθε κατάσταση (συνθήκη). Αυτό μας επιτρέπει να υπολογίσουμε πόσο συστηματική είναι η διαφορά μεταξύ των 4 τιμών Wt και των 4 τιμών για κάθε κατάσταση. Αν υποθέσουμε ότι οι τιμές προέρχονται από μια κανονική κατανομή (δεν έχουμε λόγο να μην το πιστεύουμε) μπορούμε να εφαρμόσουμε το t.test για να ελέγξουμε κατά πόσο οι μέσες τιμές των Wt, A διαφέρουν συστηματικά ή όχι. Για την πρώτη τιμή του συνόλου data αυτό μπορεί να γίνει ως εξής:
t.test(data[1,2:4], data[1,5:7])η οποία συγκρίνει τις τέσσερις τιμές Wt και τις τρεις τιμές Α για το γονίδιο A1bg (που είναι το πρώτο στον πίνακα data και το οποίο μπορείτε να δείτε ανα γράψετε data[1,1]). To αποτέλεσμα που βλέπετε στην οθόνη περιέχει μια σειρά από πληροφορίες, οι βασικότερες εκ των οποίων είναι οι δύο μέσες τιμές (κάτω-κάτω στο μήνυμα) και οι τιμές t (το στατιστικό τεστ) και κυρίως p-value. Το p-value δίνει την πιθανότητα να ισχύει η υπόθεσή ότι τα δύο δείγματα προέρχονται από κατανομές με την ίδια μέση τιμή. Αυτό που θέλουμε λοιπόν για να μπορέσουμε να πούμε πως η διαφορά είναι συστηματική είναι μικρές τιμές p-value. Γενικώς θεωρούμε ότι τιμές p-value μικρότερες από 0.05 είναι ικανοποιητικές για να πούμε ότι η διαφορά που παρατηρούμε είναι στατιστικά σημαντική.
Η επόμενη σειρά εντολών υπολογίζει το p-value για καθένα από τα 18702 γονίδια του data για τις διαφορές τιμών μεταξύ Α και Wt, και τα αποθηκεύει σε μια λίστα τιμών που ονομάζουμε pval_A. (θα πάρει λίγο χρόνο).
pval_A<-vector(mode="numeric", length=length(data[,1]));for(i in 1:length(data[,1])){pval_A[i]<-t.test(data[i,2:4], data[i,5:7])$p.value}H λίστα pval_A τώρα περιέχει την πληροφορία σχετικά με το πόσο στατιστικά σημαντική είναι κάθε τιμή logFC_A. Ας δούμε πώς μπορούμε να την χρησιμοποιήσουμε.
Εξαγωγή απορρυθμισμένων γονιδίων. Γραφική αναπαράσταση σχετικής έκφρασης
Έχουμε πλέον στα χέρια μας τις τιμές αλλαγής της έκφρασης και του πόσο σημαντική είναι αυτή για όλα τα γονίδια. Αν θέλουμε να εντοπίσουμε ποια γονίδια είναι απορρυθμισμένα δεν έχουμε παρά να ζητήσουμε από την R να βρει ποια έχουν α) πολύ μεγάλη ή πολύ μικρή τιμή logFC_A και β) τιμή pval_A μικρότερη από ένα όριο. Ένα ζεύγος συνηθισμένων τιμών-κατωφλίων είναι: logFC > 1 η logFC < -1 και p-value<0.05. Mπορούμε να εντοπίσουμε τα γονίδια που πληρούν αυτές τις προϋποθέσεις με τις παρακάτω εντολές:
which(logFC_A>=1 & pval_A<=0.05)->upwhich(logFC_A<=-1 & pval_A<=0.05)->downτο σύμβολο "&" σημαίνει το λογικό "ΚΑΙ", θέλουμε δηλαδή να ισχύουν και τα δύο προαπαιτούμενα
ενώ αντίστοιχα τα γονίδια που δεν αλλάζουν σημαντικά την εκφρασή τους θα είναι όλα τα υπόλοιπα
which(abs(logFC_A)<1 | pval_A>0.05)->nochange Τα γονίδια αυτά καθαυτά μπορούμε να τα βρούμε αν αντικαταστήσουμε τα up, down στον αρχικό πίνακα data για τη στήλη των γονιδίων. Έτσι η εντολή:
upA<-data[up,1]θα αποθηκεύσει σε μια λίστα upA τα ονόματα των γονιδίων που υπερ-εκφράζονται στην κατάσταση Α.
Τελευταίο στάδιο της Άσκησής μας είναι η γραφική παράσταση. Θα δημιουργήσουμε ένα volcano plot όπως το περιγράψαμε στο μάθημα με την παρακάτω σειρά εντολών:
plot(logFC_A[nochange], -log10(pval_A[nochange]), type="p", pch=19, col="grey", xlab="log2(FC)", ylab="-log10(p-value)", main="Deregulation A vs Wt");lines(logFC_A[up], -log10(pval_A[up]), type="p", pch=19, col="dark red");lines(logFC_A[down], -log10(pval_A[down]), type="p", pch=19, col="dark blue");legend("topleft", c("no change", "Up", "Down"), fill=c("grey","dark red", "dark blue"), bty="n")Άσκηση
1. Να υπολογίσετε τα logFC και pvalue για όλες τις συνθήκες.
2. Nα δημιουργήσετε volcano plots για κάθεμια από τις 5 καταστάσεις.
3. Να αναφέρετε τον αριθμό των κοινά απορυθμισμένων γονιδίων για κάθε κατάσταση. Χρησιμοποιείστε για αυτό το σκοπό τη συνάρτηση intersect() της R με τον ακόλουθο τρόπο:
upA<-data[x,1]upB<-data[y,1]commonAB<-intersect(upA,upB)noofcommonABgenes<-length(commonAB)Να καταγράψετε τους αριθμούς των κοινών απορρυθμισμένων για κάθε κατάσταση σε δύο πίνακες 5Χ5 (ξεχωριστά για υπερ- και υπο-εκφραζόμενα):
H αναφορά σας δε θα πρέπει να ξεπερνάει τις τρεις σελίδες. Ζητούμε ένα αρχείο κειμένου στο οποίο θα περιγράφονται τα βήματα της διαδικασίας που ακολουθήσατε και τα αποτελέσματα στις παραπάνω ερωτήσεις σαφώς διατυπωμένα. Πίνακες/διαγράμματα που συνοδεύουν τις αναλύσεις σας θα πρέπει να βρίσκονται μέσα στην αναφορά και όχι σε ξεχωριστά αρχεία. Ονομάστε το αρχείο χρησιμοποιώντας το ονοματεπώνυμο και τον ΑΜ σας.
Ανεβάστε την αναφορά σας στο παρακάτω link
https://www.dropbox.com/request/7iUw2iIxD5ok71ahf5yt
μέχρι τη Κυριακή 10/5/2020 στις 23.59 π.μ.