Θα ξεκινήσουμε τη σημερινή άσκηση με μια αναδρομή στην προηγούμενη. Θυμηθείτε πως στα πιο δύσκολα ερωτήματα ήταν να κάνετε μια καταμέτρηση των transcripts ανά χρωμόσωμα και ανα γονίδιο, ενώ γενικό ερώτημα ήταν να μετρήσετε το μέγεθος του αρχείου συνολικά.
Σε προγράμματα όπως το Excel ή άλλα αντίστοιχα προγράμματα επεξεργασίας λογιστικών φύλλων κάτι τέτοιο είναι εφικτό μέσα από συνδυασμό πράξεων, εφαρμογών συναρτήσεων, χρήση κελιών κλπ. Όμως υπάρχουν εναλλακτικές μέθοδοι που κάνουν τέτοιες πράξεις (και πολλές άλλες) πολύ πιο έυκολες και για το λόγο αυτό αξίζει να εξοικειωθεί κανείς μαζί τους. Μια τέτοια είναι η R, ένα υπολογιστικό "περιβάλλον" που επιτρέπει τεράστιο αριθμό αναλύσεων, στατιστικής επεξεργασίας και γραφικής αναπαράστασης. Μπορείτε να μάθετε λίγα περισσότερα για την R καθώς και να δείτε πως την εγκαθιστά κανείς σε PC εδώ.
Μπορείτε επίσης να δείτε μια εισαγωγή στην R στην παρουσίαση του Νίκου Φίκα εδώ: https://prezi.com/view/dnqH1E39FQ7Es6MNkHmb/
Στα πλαίσια μια απλής εισαγωγής ώρα θα δούμε πώς μπορούμε να απαντήσουμε στα ερωτήματα της προηγούμενης άσκησης με τη χρήση απλών εντολών της R. Κατ' αρχήν θα κατεβάσουμε το αρχείο με τα RefSeq γονίδια από τα συνημμένα αρχεία του Practical 1 και θα το αποθηκεύσουμε σε έναν φάκελο στον υπολογιστή μας. Στη συνέχεια θα ανοίξουμε την R από τον υπολογιστή μας και θα πλοηγηθούμε σε αυτόν τον φάκελο ως εξής:
Άνοιγμα R (ή R-studio)
setwd("C:/pathtofolder/")...Mε αυτόν τον τρόπο η R "βρίσκεται" στον φάκελο όπου περιέχονται τα δεδομένα τα οποία μπορούμε τώρα να διαβάσουμε σε μια μεταβλητή-πίνακα με την εντολή.
refseq<-read.delim("Human_genes_RefSeq.bed", header=T, sep="\t")Στην παραπάνω εντολή υποθέτουμε ότι υπάρχει ένα αρχείο που περιέχει σε μορφή bed τα δεδομένα σας. Αν δεν το έχετε μπορείτε να το κατεβάσετε εδώ.
δειτε τις πρώτες γραμμές του αρχείου ζητώντας να δείτε την αρχή της μεταβλητής refseq
head(refseq)Πώς θα υπολογίσουμε τώρα τον αριθμό των transcripts που περιέχει κάθε χρωμόσωμα; Πολύ απλά με τη χρήση της συνάρτησης table πάνω στη στήλη της μεταβλητής που μας ενδιαφέρει. Η στήλη με τα χρωμοσώματα συμβολίζεται ως refseq$chrom. Η εκτέλεση της εντολής:
table(refseq$chrom)Μας δίνει το αποτέλεσμα που επιθυμούμε, το οποίο μπορούμε να αποθηκεύσουμε σε μια μεταβλητή/πίνακα (με τη χρήση του "<-")
chromosomes<-table(refseq$chrom)Μπορούμε επιπλέον να κατατάξουμε αυτή τη λίστα από το μικρότερο (σε αριθμό) στο μεγαλύτερο χρωμόσωμα:
sort(chromosomes)->chromosomesνα βρουμε πόσα χρωμοσώματα έχουμε:
n<-length(chromosomes); nκαι να δούμε το αποτέλεσμα γραφικά σε μια πίτα:
pie(chromosomes, col=rainbow(51))ή σε ραβδόγραμμα (με μια λίγο πιο σύνθετη εντολή):
par(mar=c(5,10,3,3));barplot(chromosomes, col=rainbow(51), las=1, xlab="number of genes", main="chromosome name", cex.names=0.8, horiz=T)Έχοντας δει κάποια εισαγωγικά για τη χρήση της R, sτη συνέχεια θα δούμε πως μπορούμε να χρησιμοποιήσουμε την R για να αναλύσουμε όχι μόνο γονιδιωματικά δεδομένα συντεταγμένων αλλά και αλληλουχιών. Για το σκοπό αυτό θα χρησιμοποιήσουμε:
α. Αλληλουχίες που έχουμε αποκομίσει από βάσεις δεδομένων γονιδιωματικής
β. Συναρτήσεις γραμμένες στην R σε συνδυασμό με ενσωματωμένες συναρτήσεις.
Για τη σημερινή άσκηση θα πρέπει αρχικά να κατεβάσετε στον υπολογιστή σας τις συναρτήσεις readfastafile.R και gcContent.R που θα βρείτε στα παρακάτω link:
συνημμένες στο τέλος της σελίδας. Θα τις χρησιμοποιήσετε στη συνέχεια για διαβάσουμε τα αρχεία αλληλουχιών που θα χρησιμοποιήσουμε και να υπολογίσουμε το GC content τους. Η χρήση των συναρτήσεων στην R γίνεται με την απλή χρήση της συνάρτησης source(). Αφού κατεβάσετε τις συναρτήσεις στον υπολογιστή σας, δεν έχετε παρά να γράψετε στην R:
source("..../readfastafile.R") source("..../gcContent.R")Όπου αντί για "...." στις παραπάνω εντολές θα έχετε γράψει το πλήρες path για τον φάκελο στον οποίον έχετε σώσει τις συναρτήσεις, π.χ. :
source("/home/Christoforos/Documents/RPrograms/gcContent.R")Έχοντας ετοιμάσει το περιβάλλον και με τις κατάλληλες συναρτήσεις προς εκτέλεση μπορούμε να περάσουμε στο κυρίως πρόβλημα. Αυτό θα είναι να αναλύσουμε δύο κατηγορίες αλληλουχιών
Στη συνέχεια αυτό που χρειάζεται να κάνουμε είναι να διαβάσουμε κατ' αρχήν δύο συλλογές αλληλουχιών που μπορείτε να κατεβάσετε από το τέλος αυτής της σελίδας. Αφού κατεβάσετε τα δύο αρχεία SCGenome_GeneSeqs.fa και SCGenome_RegSeqs.fa να τα μεταφέρετε στο φάκελο από όπου εκτελείτε την R και στη συνέχεια να τα διαβάσετε σε δύο μεταβλητές με τις παρακάτω εντολές:
genes<-readfastafile("SCGenome_GeneSeqs.fa")regseq<-readfastafile("SCGenome_RegSeqs.fa")Τα αρχεία αυτά περιέχουν όλα τα γονίδια και όλες τις ρυθμιστικές αλληλουχίες του Σακχαρομύκητα έτσι όπως έχουν καταλογογραφηθεί από την βάση δεδομένων του S. cerevisiae (SGD). Mπορείτε να δείτε την αρχή κάθε αρχείου (καθώς είναι αρκετά μεγάλα) με τη χρήση της συναρτησης head().
head(genes)Eνώ για να δείτε πόσες αλληλουχίες περιέχονται σε κάθε συλλογή μπορείτε να τρέξετε τις εντολές:
ngenes<-length(genes) nregseq<-length(regseq)και έπειτα να καλέσετε τις μεταβλητές για προβολή στην οθόνη:
ngenesnregseqΣτη συνέχεια θα υπολογίσουμε τη σύσταση σε βάσεις γουανίνης+κυτοσίνης για κάθεμιά από τις δύο συλλογές. Αυτό θα γίνει με τη χρήση μια συνάρτησης του πακέτου seqinr που ονομάζεται (εύλογα) GC την οποία όμως θα χρησιμοποιήσουμε με επαναληπτικό τρόπο. Εκτελέστε τις παρακάτω εντολές για τις δύο συλλογές:
gcgenes<-gcContent(genes)gcreg<-gcContent(regseq)H συνάρτηση gcContent() παίρνει κάθε αλληλουχία της συλλογής και υπολογίζει το μήκος της και το ποσοστό βάσεων (G+C), παρουσιάζοντας τα δεδομένα αυτά σε μορφή πίνακα. Μπορείτε να δείτε τα περιεχόμενα των δύο αυτών πινάκων και παλι μέσα από τις πρώτες γραμμές τους:
head(gcgenes)ή με τη χρήση της συνάρτησης str() που σας δίνει πληροφορίες για το μέγεθος και το είδος του πίνακα
str(gcgenes)Από την παραπάνω βλέπουμε ότι ο gcgenes είναι ενας πίνακας της μορφής data.frame που περιέχει τρεις στήλες: όνομα, GC% και μήκος της αλληλουχίας για 6692 αλληλουχίες. Αν επιθυμούμε να δούμε ένα στοιχείο από τον κάθε πίνακα (δηλ. μόνο το GC ή μόνο τα μήκη των αλληλουχιών) μπορούμε να το πάρουμε δηλώνοντας το όνομά του, με το σύμβολο του δολλαρίου:
gcgenes$GCή βάζοντας σε αγκύλες τον αριθμό της στήλης μετά από ένα κόμμα π.χ.:
gcgenes[,2]Προσέξτε τη σύνταξη Χ[,2] που σημαίνει "δεύτερη στήλη του Χ" σε αντίθεση με την Χ[2,] που σημαίνει "δεύτερη γραμμή του Χ". Έχοντας στα χέρια μας δύο πίνακες που ονομάζονται gcgenes και gcreg τις κατανομές των τιμών GC για το σύνολο των αλληλουχιών της κάθε συλλογής. Πώς θα μπορούσαμε να τις συγκρίνουμε; Αρχικά με δύο ιστογράμματα το ένα κάτω από το άλλο:
par(mfrow=c(2,1)); # plotting two plots one below the otherhist(gcgenes[,2], breaks=100, xlim=c(0.2, 0.7), main="Genes", col="red")hist(gcreg[,2], breaks=100, xlim=c(0.2, 0.7), main="Regulatory", col="blue")Το ίδιο μπορεί να γίνει πιο εύκολα με μια απλή εφαρμογή μιας συνάρτησης "θηκογράμματος" (boxplot) που αναπαριστά τις κατανομές δύο ή περισσοτέρων συνόλων με τη μορφή "θηκών" (boxes) τα όρια των οποίων αντιστοιχούν σε χαρακτηριστικές τιμές της κάθε κατανομής:
par(mfrow=c(1,1)); # plotting back to 1x1boxplot(gcgenes[,2], gcreg[,2], notch=T, col=c("red","blue"), lwd=3, names=c("Genes", "Regulatory"))Πώς μπορούμε τέλος να δούμε τις χαρακτηριστικές τιμές των κατανομών; Μπορούμε να υπολογίσουμε τη μέση τιμή και την τυπική απόκλιση της κάθε συλλογής με τις αντίστοιχες συναρτήσεις mean(), sd():
mgcgenes<-mean(gcgenes$GC)sdgcgenes<-sd(gcgenes$GC)mgcreg<-mean(gcreg$GC)sdgcreg<-sd(gcreg$GC)Τις παραπάνω τιμές μπορούμε να χρησιμοποιήσουμε τώρα για να συγκρίνουμε με άγνωστες αλληλουχίες. Η εργασία που καλείστε να κάνετε βασίζεται στις παραπάνω αναλύσεις.
Ερωτήσεις:
1. Να κατεβάσετε τις πέντε αλληλουχίες (test1.fa-test5.fa) από εδώ
2. Να διαβάσετε καθεμία στην R, να την σώσετε σε μια μεταβλητή και στη συνέχεια να μετρήσετε το GC% με τη χρήση των συναρτήσεων readfastafile() και gcContent() όπως στα παραδείγματα που είδατε παραπάνω.
3. Nα υπολογίσετε τα z-scores για καθεμία από τις 5 αλληλουχίες σε σχέση με τις δύο κατανομές gcgenes, gcreg όπως τις υπολογίσατε παραπάνω. Θυμηθείτε ότι για καθεμία από τις δύο κατανομές οι μέσες τιμές και τυπικές αποκλίσεις διαφέρουν συνεπώς και οι τιμές GC-mean(GC)/sd(GC) θα διαφέρουν. Oι πράξεις μπορούν να γίνουν στην R με πολύ απλό τρόπο χρησιμοποιώντας το "/" ως τελεστή διαίρεσης. Για παράδειγμα:
z1<-(gc1$GC-mgcgenes)/sdgcgenes ή z1<-(gc1[,2]-mgcgenes)/sdgcgenes
H παραπάνω πράξη δίνει το z-score για την αλληλουχία test1.fa σε σχέση με το GC content των γονιδίων. Αφού υπολογίσετε και τα 5x2=10 z-scores να διακινδυνέψετε μια πρόβλεψη για το ποιες από τις πέντε αλληλουχίες είναι γονίδια και ποιες όχι και να δικαιολογήσετε την πρόβλεψή σας.
4. Bonus points αν βρειτε έναν τρόπο να αναπαραστήσετε γραφικά τις τιμές gc για τις άγνωστες αλληλουχίες πάνω στο ιστόγραμμα των gcgenes. Μπορείτε να δοκιμάσετε τις συναρτήσεις rug() και abline() της R.
H αναφορά σας δε θα πρέπει να ξεπερνάει τις τρεις σελίδες. Ζητούμε ένα αρχείο κειμένου στο οποίο θα περιγράφονται τα βήματα της διαδικασίας που ακολουθήσατε και τα αποτελέσματα στις παραπάνω ερωτήσεις σαφώς διατυπωμένα. Πίνακες/διαγράμματα που συνοδεύουν τις αναλύσεις σας θα πρέπει να βρίσκονται μέσα στην αναφορά και όχι σε ξεχωριστά αρχεία. Ονομάστε το αρχείο χρησιμοποιώντας το ονοματεπώνυμο και τον ΑΜ σας.
Ανεβάστε την αναφορά σας στο παρακάτω link:
https://www.dropbox.com/request/b8cu2iTmLrLSQgxC5GxB
μέχρι την Κυριακή 8/3/2020 και ώρα 23:59.