Εισαγωγή

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

Χρήση παραγόντων (factors)

Γενικά για τους παράγοντες

Οι παράγοντες είναι στοιχεία των dataframes που αντιστοιχούν σε κατηγορικές μεταβλητές. Οι παράγοντες είναι στην ουσία διανύσματα χαρακτήρων που όμως μεταφράζονται πάντα σε κατηγορικά δεδομένα όταν χειριζόμαστε πλαίσια δεδομένων. Στο παρακάτω παράδειγμα βλέπουμε τη λειτουργία των παραγόντων στο ενσωματωμένο στην R πλαίσιο δεδομένων iris:

str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
is.factor(iris$Species)
## [1] TRUE

Στην τελευταία γραμμή του output της \(str()\) βλέπουμε ότι το διάνυσμα Species που περιέχει τα ονόματα των ειδών ίριδος αναγνωρίζεται ως factor με 3 επίπεδα (levels) τα οποία αντιστοιχούν στις τρείς διαφορετικές τιμές που παίρνει το διάνυσμα Species. Επί της ουσίας όλα τα διανύσματα που είναι factors είναι διανύσματα χαρακτήρων ωστόσο εφόσον χαρακτηριστούν ως factors η συνάρτηση \(is.character()\) δεν επαληθεύεται:

is.character(iris$Species)
## [1] FALSE

Για ποιο λόγο η R μετατρέπει τα διανύσματα χαρακτήρων σε factors; Αυτό συμβαίνει επειδή οι παράγοντες είναι ιδιαίτερα χρήσιμοι για μια σειρά από χειρισμούς σε πλαίσια δεδομένων.

Παρά το γεγονός ότι η R θα αποδώσει την ιδιότητα παράγοντα σε κάθε διάνυσμα χαρακτήρων μέσα σε ένα dataframe, αυτό έχει νόημα μόνο στην περίπτωση που ο αριθμός των μοναδικών τιμών μέσα στο διάνυσμα είναι περιορισμένος ώστε να επιτρέπει την ομαδοποίηση. Το πλήθος των μοναδικών τιμών μπορεί κανείς να δει με τη χρήση της συνάρτησης \(levels()\) η οποία τις επιστρέφει σε ένα διάνυσμα χαρακτήρων:

levels(iris$Species)
## [1] "setosa"     "versicolor" "virginica"

Στην πράξη οποιοδήποτε διάνυσμα χαρακτήρων μπορεί να μετατραπεί σε παράγοντα με τη χρήση της συνάρτησης \(factor()\):

a_vector<-c(rep("X",5),rep("Y",3), rep("Z",1))
a_factor<-factor(a_vector)
levels(a_factor)
## [1] "X" "Y" "Z"

Χρήση επιπέδων (levels) σε παράγοντες

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

# Data in vectors
income <- c(32,51,12,23,26,27,24,26,25,18,30,22,24,21,25,27,18)
age <- c(48,59,26,23,37,32,20,25,45,55,44,42,51,30,35,29,19)
education <- c("middle","high","basic","high","high","middle","high","high","middle","basic", "high", "basic", "basic", "middle", "middle", "middle", "basic")
ageclass<-c("higher","higher","lower","lower","higher","lower","lower", "lower","higher","higher","higher","higher","higher","lower","higher","lower",  "lower") 
# Create the dataframe.
income_data <- data.frame(income,age,education,ageclass)
levels(income_data$education)
## [1] "basic"  "high"   "middle"

Σε αρκετές περιπτώσεις θέλουμε μια διαφορετική σειρά στα επίπεδα του παράγοντα. Αυτό γίνεται με την αντικατάσταση του παράγοντα από έναν παράγοντα όπου τα επίπεδα ορίζονται αναλυτικά όπως στο παράδειγμα:

income_data$education<-factor(income_data$education, levels=c("basic", "middle", "high"))
levels(income_data$education)
## [1] "basic"  "middle" "high"

Χρήση Παραγόντων σε Πλαίσια Δεδομένων

Μια πρώτη συνάρτηση που δρα σε παράγοντες είναι η \(table()\) που επενεργεί αποκλειστικά σε διανύσματα παραγόντων για να μας δώσει το πλήθος του κάθε level στο διάνυσμα:

table(income_data$education)
## 
##  basic middle   high 
##      5      6      6

Η R επιτρέπει την κατηγοριοποίηση μέσω πολλαπλών παραγόντων με τη χρήση της συνάρτησης \(table()\). Η συγκεκριμένη συνάρτηση μπορεί να δεχτεί ως είσοδο δύο ή περισσότερα διανύσματα παραγόντων και να δημιουργήσει έναν πίνακα “σύμπτωσης” (contingency table) που περιέχει τον αριθμό των στοιχείων των συνδυασμών τους.

table(income_data$education, income_data$ageclass)
##         
##          higher lower
##   basic       3     2
##   middle      3     3
##   high        3     3

Η βασική χρησιμότητα των παραγόντων είναι η κατηγοριοποίηση των δεδομένων σε ένα dataframe. Η R αναγνωρίζει τις μεταβλητές που μπορούν να χρησιμοποιηθούν ως παράγοντες και μπορεί να κατατάξει τις υπόλοιπες μεταβλητές με βάση τα επίπεδά τους πριν κάνει υπολογισμούς σε διανυσματικό επίπεδο. Χαρακτηριστικές συναρτήσεις που επιτρέπουν υπολογισμούς αυτού του είδους είναι η \(by()\) και η \(aggregate()\). Ας δούμε πώς μπορούμε να εκμεταλλευτούμε τους παράγοντες με αυτές στο παράδειγμα του income_data που δημιουργήσαμε προηγουμένως:

by(income_data$income, income_data$education, mean)
## income_data$education: basic
## [1] 18.8
## ------------------------------------------------------------ 
## income_data$education: middle
## [1] 26.16667
## ------------------------------------------------------------ 
## income_data$education: high
## [1] 30

Στο πρώτο παράδειγμα η \(by(a,b,c)\) ορίζει αρχικά το διάνυσμα στο οποίο θέλει να πραγματοποιήσει τον υπολογισμό, στη δεύτερη θέση το διάνυσμα παραγόντων με το οποίο θα κάνει την κατηγοριοποίηση και στην τρίτη θέση τη συνάρτηση που θα εφαρμοστεί για τους υπολογισμούς. Το αποτέλεσμα εδώ είναι η μέση τιμή του εισοδήματος ανά μορφωτικό επίπεδο. Αντίστοιχη λειτουργία πραγματοποιεί η \(aggregate()\) με τη βασική διαφορά ότι τα διανύσματα παραγόντων που θα χρησιμοποιηθούν δίνονται με τη μορφή λίστας. Στο παρακάτω παράδειγμα υπολογίζουμε τη μέση τιμή των ηλικιών ανά μορφωτικό επίπεδο:

aggregate(income_data$age, by=list(Education=education), mean)
##   Education        x
## 1     basic 38.60000
## 2      high 34.66667
## 3    middle 36.50000

Ένα σημαντικό πλεονέκτημα της \(aggregate()\) είναι ότι μπορεί να εφαρμοστεί σε συνδυασμούς παραγόντων που οδηγούν σε πιο πολύπλοκες κατηγοριοποιήσεις. Στη συνέχεια βλέπουμε πώς μπορούμε να χρησιμοποιήσουμε την \(aggregate()\) για να υπολογίσουμε χαρακτηριστικές τιμές για αυτούς τους συνδυασμούς, κάτι που δεν μπορούμε να κάνουμε με την \(by()\).

aggregate(income_data$income, by=list(Education=education, Age=ageclass), mean)
##   Education    Age        x
## 1     basic higher 21.33333
## 2      high higher 35.66667
## 3    middle higher 27.33333
## 4     basic  lower 15.00000
## 5      high  lower 24.33333
## 6    middle  lower 25.00000

Από το οποίο προκύπτει ότι υψηλότερο μέσο εισόδημα έχουν τα μεγαλύτερης ηλικίας και μορφωτικού επιπέδου άτομα.

Λίστες στην R

Οι λίστες (lists) είναι η πιο πολύπλοκη κατηγορία δεδομένων στην R. Αποτελούνται από σύνολα διανυσμάτων που δεν είναι υποχρεωτικό να είναι του ίδιου τύπου αλλά ούτε και του ίδιου μήκους. Στην πιο διευρυμένη μορφή τους μπορούν να περιέχουν ακόμα και πίνακες ή και συναρτήσεις ως στοιχεία. Μια λίστα δημιουργείται με τη χρήση της ομώνυμης συνάρτησης, η οποία δέχεται μια παράθεση στοιχείων, όπως στο παράδειγμα:

l<-list(c(1,3,5), c("Steven", "Robbie", "Kenny", "Ian"), TRUE, c(0.9, 0.8))
l
## [[1]]
## [1] 1 3 5
## 
## [[2]]
## [1] "Steven" "Robbie" "Kenny"  "Ian"   
## 
## [[3]]
## [1] TRUE
## 
## [[4]]
## [1] 0.9 0.8

Βλέπουμε ότι η λίστα απαρτίζεται από τις παραμέτρους που περάσαμε στην \(list()\) με τη σειρά που το κάναμε. Για να δημιουργήσουμε μια πιο καλά οργανωμένη εκδοχή της θα δώσουμε ονόματα στα στοιχεία της \(l\) με τη χρήση της \(names\):

names(l)<-c("odds", "names", "logical", "decimals")
l
## $odds
## [1] 1 3 5
## 
## $names
## [1] "Steven" "Robbie" "Kenny"  "Ian"   
## 
## $logical
## [1] TRUE
## 
## $decimals
## [1] 0.9 0.8

Προσέξτε ότι τώρα τα ονόματα των στοιχείων της λίστας βρίσκονται στη θέση που προηγουμένων υπήρχε αρίθμηση σε διπλές αγκύλες. Ίσως το πιο σημαντικό σημείο σε ό,τι αφορά τις λίστες είναι η χρήση αυτών των διπλών αγκυλών. Είτε φέρουν ονόματα, είτε όχι, η θέση των στοιχείων μιας λίστας δίνεται από έναν αριθμό μέσα σε διπλές αγκύλες αντί για μονές όπως ισχύει για τα διανύσματα και τους πίνακες. Η διαφορά στη χρήση μονών και διπλών αγκυλών είναι αρκετά λεπτή και μπορεί να διαφύγει σε πρώτο επίπεδο. Προσέξτε τη διαφορά μεταξύ:

l[1]
## $odds
## [1] 1 3 5

και:

l[[1]]
## [1] 1 3 5

Πρακτικά δε φαίνεται να υπάρχει καμία. Δείτε όμως τι είδους αντικείμενο κρύβεται πίσω από τις δύο αυτές δομές:

class(l[1])
## [1] "list"
class(l[[1]])
## [1] "numeric"

Χρήση δηλαδή των μονών αγκυλών επιστρέφει ένα αντικείμενο τύπου λίστας ενώ οι διπλές αγκύλες επιστρέφουν το διάνυσμα του πρώτου στοιχείου. Αυτό είναι πολύ σημαντικό καθώς όπως έχουμε συζητήσει και παραπάνω συγκεκριμένες συναρτήσεις μπορούν να εφαρμοστούν μόνο σε συγκεκριμένους τύπους αντικειμένων. Πώς μπορούμε τώρα να εντοπίσουμε ένα στοιχείο εντός των διανυσμάτων που απαρτίζουν μια λίστα; Ακολουθώντας τη λογική που έχουμε δει ως τώρα, αρκεί να προσδιορίσουμε αρχικά το διάνυσμα που θέλουμε με διπλές αγκύλες και στη συνέχεια να ορίσουμε σε αυτό την επιθυμητή θέση με μονές. Έτσι π.χ. το τρίτο όνομα στο στοιχείο \(names\) της \(l\) δίνεται ως εξής:

l[[2]][3]
## [1] "Kenny"

Εναλλακτικά, η πρόσβαση σε αυτό μπορεί να γίνει και με τη χρήση του ονόματος του διανύσματος όπως είδαμε στα dataframes:

l$names[3]
## [1] "Kenny"

Στην περίπτωση που θέλουμε να πάρουμε πίσω όλη τη λίστα σε ένα διάνυσμα μπορούμε να την “καταρρεύσουμε” σε διάνυσμα με την συνάρτηση \(unlist()\):

unlist(l)
##     odds1     odds2     odds3    names1    names2    names3    names4   logical 
##       "1"       "3"       "5"  "Steven"  "Robbie"   "Kenny"     "Ian"    "TRUE" 
## decimals1 decimals2 
##     "0.9"     "0.8"

Όπως φαίνεται, στην περίπτωση ονοματισμένης λίστας, τα ονόματα διατηρούνται και στο διάνυσμα, κάτι που είναι αρκετά βολικό.

Χειρισμός Λιστών με συναρτήσεις τύπου apply()

Οι χειρισμοί των λιστών μοιάζουν πολύ με τους αντίστοιχους των διανυσμάτων. Έτσι μπορούμε να ενώσουμε δύο ή περισσότερες λίστες με την συνάρτηση \(c()\):

l1<-list(c("A","B"), c(1,2,3))
l2<-list(c("C","D","E"), c(T,F,F,F))
l3<-c(l1,l2)
l3
## [[1]]
## [1] "A" "B"
## 
## [[2]]
## [1] 1 2 3
## 
## [[3]]
## [1] "C" "D" "E"
## 
## [[4]]
## [1]  TRUE FALSE FALSE FALSE

Η βασικότερη χρησιμότητα των λιστών είναι η δυνατότητα εφαρμογής συναρτήσεων στα στοιχεία τους με τη χρήση των συναρτήσεων \(apply()\) που λειτουργούν καλώντας άλλες συναρτήσεις και εφαρμόζοντάς τις απευθείας, με αποτέλεσμα την ταχύτατη εκτέλεσή τους. Απαραίτητη προϋπόθεση είναι τα στοιχεία της λίστας να μπορούν να αποτελέσουν όρισμα της συνάρτησης. Στην περίπτωση των λιστών, το πλεονέκτημα που έχουμε σε σχέση με τους πίνακες είναι ότι τα διανύσματα που μπορούμε να αναλύσουμε δεν είναι υποχρεωτικά ίδιου μήκους. Έστω για παράδειγμα ότι έχουμε μια σειρά από διανύσματα των οποίων θέλουμε να υπολογίσουμε τις μέγιστες τιμές. Μπορούμε να αποφύγουμε μια επαναληπτική κλήση της αντίστοιχης συνάρτησης \(max()\) ενσωματώνοντας τα διανύσματα σε μια λίστα και καλώντας την \(max()\) σε μία κλήση της \(sapply()\):

x<-c(1, 8, 9, 12, 5)
y<-c(6, 7, 0.8, 2, 5.1)
z<-c(-10, 5, 2, 3.8)
my_list<-list(x,y,z)
sapply(my_list, FUN=max)
## [1] 12  7  5

Το αποτέλεσμα είναι ένα διάνυσμα τριών τιμών που αντιστοιχεί στις μέγιστες τιμές των \(x\), \(y\) και \(z\). Το ίδιο αποτέλεσμα καταχωρημένο σε λίστα προκύπτει αν αντί για την \(sapply()\) καλέσουμε την \(lapply()\):

lapply(my_list, FUN=max)
## [[1]]
## [1] 12
## 
## [[2]]
## [1] 7
## 
## [[3]]
## [1] 5

Analysis of Variance ANOVA

Με τον όρο Ανάλυση Διακύμανσης ή Ανάλυση Διασποράς (Analysis of Variance) αναφερόμαστε σε μια δέσμη από στατιστικές αναλύσεις μοντελοποίησης που σκοπό έχουν την εκτίμηση της διαφοράς των μέσων τιμών μεταξύ ομάδων σε ευρύτερα δείγματα. Όπως και πολλά άλλα μεθοδολογικά εργαλεία, η ANOVA αναπτύχθηκε από τον Ronald Fisher και στηρίζεται στην διαισθητικά απλή αλλά τεχνικά μάλλον πολύπλοκη βασική αρχή της ανάλυσης της παρατηρούμενης διακύμανσης μιας μεταβλητής σε διαφορετικές συνιστώσες που σχετίζονται με εγγενείς τάσεις και τυχαία σφάλματα. Η ANOVA είναι η μεθοδολογική προσέγγιση επιλογής για την σύγκριση περισσότερων από δύο δείγματα. Εννοιολογικά δε διαφέρει πολύ από μια προσέγγιση πολλαπλών ζευγαρωτών t-test, ωστόσο είναι γενικά πιο αυστηρή και έτσι οδηγεί σε πιο συντηρητικά αποτελέσματα με περιορισμένα σφάλματα τύπου Ι, οδηγεί δηλαδή σπανιότερα στην απόρριψη μηδενικών υποθέσεων που είναι αληθείς.

Η βάση της ANOVA είναι η εκτίμηση της διαφοράς στη μέση τιμή μιας παραμετρικής μεταβλητής σε σχέση με μια ή περισσότερες κατηγορικές μεταβλητές. Ανάλογα με το σχεδιασμό του πειράματος, των αριθμών των μεταβλητών και των ομάδων δειγμάτων υπάρχουν διάφορες παραλλαγές της ANOVA. Έτσι υπάρχει μονο- (one-way) και πολυ-παραγοντική ανάλυση (n-way ANOVA) ανάλογα με το αν η διακύμανση μελετάται σε σχέση με μία ή περισσότερες κατηγορικές μεταβλητές, η πολυμεταβλητή (mulitvariate ANOVA ή MANOVA) όταν η μετρούμενες μεταβλητές απόκρισεις είναι παραπάνω από μία, ενώ ανάλογα με το αν υπάρχουν επιπλέον παραμετρικές μεταβλητές που μπορούν να επηρεάσουν την μετρούμενη (target variable) χρησιμοποιείται ένας συνδυασμός ΑNOVA και παλινδρόμησης (regression) που ονομάζεται Analysis of Co-Variance.

Σε πρακτικό επίπεδο, παρότι χρησιμοποιείται στον έλεγχο υποθέσεων, η ANOVA αποτελεί στην ουσία μια κατηγορία μοντέλων που ονομάζονται γενικευμένα γραμμικά μοντέλα (generalized linear models). Σημείο εκκίνησης είναι η μηδενική υπόθεση ότι οι τιμές όλων των ομάδων/κατηγοριών που εξετάζονται είναι τυχαία δείγματα που προέρχονται από την ίδια κατανομή. Η διαφορά με άλλους (ζευγαρωτούς) ελέγχους υποθέσεων είναι ότι η ANOVA λειτουργεί προσπαθώντας να μοντελοποιήσει την σχέση μεταξύ της μεταβλητής που μετρούμε, που ονομάζεται και μεταβλητή απόκρισης (response variable) και των κατηγορικών μεταβλητών που ονομάζονται και επεξηγηματικές (explanatory variables). Ο στόχος του μοντέλου είναι να εκτιμηθεί η σχέση αυτή σε επίπεδο διασποράς.

H R διαθέτει διάφορες συναρτήσεις για την διενέργεια της ANOVA. Σε αυτή την φάση θα εξετάσουμε την \(aov()\), η οποία ακολουθεί την γενική σύνταξη των γραμμικών μοντέλων που θα δούμε πιο αναλυτικά σε επόμενα κεφάλαια. Συνοπτικά η σύνταξή της είναι:

aov(ResponseVariable ~ A + B + C …, data=data.frame)

όπου ResponseVariable είναι η παραμετρική μεταβλητή που μας ενδιαφέρει να ελέγξουμε και Α, Β, C, κλπ είναι οι κατηγορικές μεταβλητές στη βάση των οποίων θέλουμε να κάνουμε τον έλεγχο. Σε μια προσέγγιση, στην οποία οι κατηγορικές μεταβλητές είναι ανεξάρτητες, αυτές δίνονται με τον τελεστή “+” ανάμεσά τους, ενώ σε περίπτωση που έχουμε συσχετιζόμενες κατηγορικές μεταβλητές (factorial ANOVA) μπορούμε να το δηλώσουμε με "*".

aov(ResponseVariable ~ A * B * C …, data=data.frame)

Ο πιο σημαντικός τελεστής και στις δύο περιπτώσεις είναι ο τελεστής σύνδεσης (~) ο οποίος στην R υποκαθιστά τη σχέση f() των μαθηματικών μοντέλων, υποδηλώνει δηλαδή την σχέση εξαρτημένων (αριστερά) και ανεξάρτητων (δεξιά) μεταβλητών.

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

Στην περίπτωση που θα δούμε παρακάτω μελετάμε τη μάζα ζώων (πουλερικών) στη βάση διαφορετικών διατροφικών συμπληρωμάτων. Εδώ έχουμε την απλούστερη περίπτωση μιας μεταβλητής απόκρισης (βάρος) και μιας κατηγορικής μεταβλητής (διατροφικό συμπλήρωμα). Με βάση τα παραπάνω η κλήση της συνάρτησης \(aov()\) θα είναι:

aov(weight~feed, data=chickwts)
## Call:
##    aov(formula = weight ~ feed, data = chickwts)
## 
## Terms:
##                     feed Residuals
## Sum of Squares  231129.2  195556.0
## Deg. of Freedom        5        65
## 
## Residual standard error: 54.85029
## Estimated effects may be unbalanced

Το output από μόνο του δε φαίνεται να μας λέει πολλά καθώς αναγράφει μόνο μέρος των στατιστικών που έχει υπολογίσει η ANOVA χωρς να γίνεται λόγος για στατιστική σημαντικότητα. Τις πληροφορίες αυτές παίρνουμε αν αποθηκεύσουμε το αποτέλεσμα της ανάλυσης σε μια μεταβλητή και εκτελέσουμε πάνω της την \(summary()\):

aov(weight~feed, data=chickwts)->fit
summary(fit)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## feed         5 231129   46226   15.37 5.94e-10 ***
## Residuals   65 195556    3009                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

To αποτέλεσμα σε αυτήν την περίπτωση είναι αρκετά πιο πλούσιο σε πληροφορία. Οι βαθμοί ελευθερίας υπήρχαν και πριν και είναι 5 για την κατηγορική μεταβλητή \(feed\) (N-1, με N=6 διατροφικά συμπληρώματα) ενώ η αντίστοιχη τιμή df για τα σφάλματα είναι 65, όπως προκύπτει αν αφαιρέσουμε από το συνολικό πλήθος των τιμών (Ν-1, για Ν=71) τους βαθμούς ελευθερίας που έχουν ενσωματωθεί στην κατηγορική μεταβλητή (Ν=5).

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

Προετομασία αρχείων για διενέργεια ANOVA

Οι συναρτήσεις που περιγράφηκαν παραπάνω χρειάζονται σαν είσοδο μια συγκεκριμένη μορφή data-frame στην οποία κάθε κατηγορική μεταβλητή θα είναι ξεχωριστό factor. Σε γενικές γραμμές τα αρχεία ήδη έχουν αυτήν την δομή, ωστόσο αν δεν συμβαίνει κάτι τέτοιο μπορείτε εύκολα να προσθέσετε την κατηγορική μεταβλητή μόνοι/ες σας.

Χειρισμός αρχείων fasta (+ περισσότερα για τις λίστες)

Tα αρχεία fasta είναι αρχεία ειδικού μορφότυπου με τον οποίο αναπαριστούμε πρωτοταγείς βιολογικές αλληλουχίες (νουκλεϊκών οξέων ή πρωτεϊνών). Αποτελούνται από ένα (single fasta) ή περισσότερα (multi-fasta) block γραμμών που αναλύονται πάντα σε: 1. μια πρώτη γραμμή που ξεκινά με “>” και περιέχει στοιχεία για την αλληλουχία που ακολουθεί 2. ένα block γραμμών που περιέχουν την αλληλουχία με συγκεκριμένο αριθμό χαρακτήρων ανά γραμμή (ακολουθούμενο από χαρακτήρα “αλλαγής γραμμής” στο τέλος)

Παρότι μπορούμε να γράψουμε τις δικές μας συναρτήσεις για να διαβάσουμε αρχεία fasta, υπάρχουν έτοιμες συναρτήσεις σε ειδικά πακέτα ανάλυσης αλληλουχιών όπως το seqinr

library(seqinr)

στο οποίο η συνάρτηση \(read.fasta()\) κάνει ακριβώς αυτό. Στο παράδειγμα βλέπουμε πώς εκτελουμε την \(read.fasta()\) σε ένα multi-fasta αρχείο με περισσότερες από μία αλληλουχίες

fastaseqs<-read.fasta("~/Dropbox/Teaching/EAP_BNP54/Class_Winter_2020/Exercises_Fall_2020/Exercise2/histone4.fa")

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

str(fastaseqs)
## List of 10
##  $ sp|P62805|H4_HUMAN: 'SeqFastadna' chr [1:103] "m" "s" "g" "r" ...
##   ..- attr(*, "name")= chr "sp|P62805|H4_HUMAN"
##   ..- attr(*, "Annot")= chr ">sp|P62805|H4_HUMAN Histone H4 OS=Homo sapiens GN=HIST1H4A PE=1 SV=2"
##  $ sp|P62806|H4_MOUSE: 'SeqFastadna' chr [1:103] "m" "s" "g" "r" ...
##   ..- attr(*, "name")= chr "sp|P62806|H4_MOUSE"
##   ..- attr(*, "Annot")= chr ">sp|P62806|H4_MOUSE Histone H4 OS=Mus musculus GN=Hist1h4a PE=1 SV=2"
##  $ sp|P84040|H4_DROME: 'SeqFastadna' chr [1:103] "m" "t" "g" "r" ...
##   ..- attr(*, "name")= chr "sp|P84040|H4_DROME"
##   ..- attr(*, "Annot")= chr ">sp|P84040|H4_DROME Histone H4 OS=Drosophila melanogaster GN=His4 PE=1 SV=2"
##  $ sp|P59259|H4_ARATH: 'SeqFastadna' chr [1:103] "m" "s" "g" "r" ...
##   ..- attr(*, "name")= chr "sp|P59259|H4_ARATH"
##   ..- attr(*, "Annot")= chr ">sp|P59259|H4_ARATH Histone H4 OS=Arabidopsis thaliana GN=At1g07660 PE=1 SV=2"
##  $ sp|P62803|H4_BOVIN: 'SeqFastadna' chr [1:103] "m" "s" "g" "r" ...
##   ..- attr(*, "name")= chr "sp|P62803|H4_BOVIN"
##   ..- attr(*, "Annot")= chr ">sp|P62803|H4_BOVIN Histone H4 OS=Bos taurus PE=1 SV=2"
##  $ sp|P62784|H4_CAEEL: 'SeqFastadna' chr [1:103] "m" "s" "g" "r" ...
##   ..- attr(*, "name")= chr "sp|P62784|H4_CAEEL"
##   ..- attr(*, "Annot")= chr ">sp|P62784|H4_CAEEL Histone H4 OS=Caenorhabditis elegans GN=his-1 PE=1 SV=2"
##  $ sp|P62801|H4_CHICK: 'SeqFastadna' chr [1:103] "m" "s" "g" "r" ...
##   ..- attr(*, "name")= chr "sp|P62801|H4_CHICK"
##   ..- attr(*, "Annot")= chr ">sp|P62801|H4_CHICK Histone H4 OS=Gallus gallus GN=H4-I PE=1 SV=2"
##  $ sp|P62804|H4_RAT  : 'SeqFastadna' chr [1:103] "m" "s" "g" "r" ...
##   ..- attr(*, "name")= chr "sp|P62804|H4_RAT"
##   ..- attr(*, "Annot")= chr ">sp|P62804|H4_RAT Histone H4 OS=Rattus norvegicus GN=Hist1h4b PE=1 SV=2"
##  $ sp|P02309|H4_YEAST: 'SeqFastadna' chr [1:103] "m" "s" "g" "r" ...
##   ..- attr(*, "name")= chr "sp|P02309|H4_YEAST"
##   ..- attr(*, "Annot")= chr ">sp|P02309|H4_YEAST Histone H4 OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) GN=HHF1 PE=1 SV=2"
##  $ sp|P62799|H4_XENLA: 'SeqFastadna' chr [1:103] "m" "s" "g" "r" ...
##   ..- attr(*, "name")= chr "sp|P62799|H4_XENLA"
##   ..- attr(*, "Annot")= chr ">sp|P62799|H4_XENLA Histone H4 OS=Xenopus laevis PE=1 SV=2"

για να δούμε καλύτερα την πληροφορία ανά αλληλουχία μπορούμε να εξετάσουμε το πρώτο στοιχείο της λίστας

fastaseqs[1]
## $`sp|P62805|H4_HUMAN`
##   [1] "m" "s" "g" "r" "g" "k" "g" "g" "k" "g" "l" "g" "k" "g" "g" "a" "k" "r"
##  [19] "h" "r" "k" "v" "l" "r" "d" "n" "i" "q" "g" "i" "t" "k" "p" "a" "i" "r"
##  [37] "r" "l" "a" "r" "r" "g" "g" "v" "k" "r" "i" "s" "g" "l" "i" "y" "e" "e"
##  [55] "t" "r" "g" "v" "l" "k" "v" "f" "l" "e" "n" "v" "i" "r" "d" "a" "v" "t"
##  [73] "y" "t" "e" "h" "a" "k" "r" "k" "t" "v" "t" "a" "m" "d" "v" "v" "y" "a"
##  [91] "l" "k" "r" "q" "g" "r" "t" "l" "y" "g" "f" "g" "g"
## attr(,"name")
## [1] "sp|P62805|H4_HUMAN"
## attr(,"Annot")
## [1] ">sp|P62805|H4_HUMAN Histone H4 OS=Homo sapiens GN=HIST1H4A PE=1 SV=2"
## attr(,"class")
## [1] "SeqFastadna"

το οποίο είναι περιέχει την αλληλουχία μαζί με κάποια επιπλέον στοιχεία (attributes). Αν θέλουμε να χειριστούμε την αλληλουχία μόνο, χωρίς τα επιπλέον θα δράσουμε πάνω σε αυτό με την unlist:

unlist(fastaseqs[1])
##   sp|P62805|H4_HUMAN1   sp|P62805|H4_HUMAN2   sp|P62805|H4_HUMAN3 
##                   "m"                   "s"                   "g" 
##   sp|P62805|H4_HUMAN4   sp|P62805|H4_HUMAN5   sp|P62805|H4_HUMAN6 
##                   "r"                   "g"                   "k" 
##   sp|P62805|H4_HUMAN7   sp|P62805|H4_HUMAN8   sp|P62805|H4_HUMAN9 
##                   "g"                   "g"                   "k" 
##  sp|P62805|H4_HUMAN10  sp|P62805|H4_HUMAN11  sp|P62805|H4_HUMAN12 
##                   "g"                   "l"                   "g" 
##  sp|P62805|H4_HUMAN13  sp|P62805|H4_HUMAN14  sp|P62805|H4_HUMAN15 
##                   "k"                   "g"                   "g" 
##  sp|P62805|H4_HUMAN16  sp|P62805|H4_HUMAN17  sp|P62805|H4_HUMAN18 
##                   "a"                   "k"                   "r" 
##  sp|P62805|H4_HUMAN19  sp|P62805|H4_HUMAN20  sp|P62805|H4_HUMAN21 
##                   "h"                   "r"                   "k" 
##  sp|P62805|H4_HUMAN22  sp|P62805|H4_HUMAN23  sp|P62805|H4_HUMAN24 
##                   "v"                   "l"                   "r" 
##  sp|P62805|H4_HUMAN25  sp|P62805|H4_HUMAN26  sp|P62805|H4_HUMAN27 
##                   "d"                   "n"                   "i" 
##  sp|P62805|H4_HUMAN28  sp|P62805|H4_HUMAN29  sp|P62805|H4_HUMAN30 
##                   "q"                   "g"                   "i" 
##  sp|P62805|H4_HUMAN31  sp|P62805|H4_HUMAN32  sp|P62805|H4_HUMAN33 
##                   "t"                   "k"                   "p" 
##  sp|P62805|H4_HUMAN34  sp|P62805|H4_HUMAN35  sp|P62805|H4_HUMAN36 
##                   "a"                   "i"                   "r" 
##  sp|P62805|H4_HUMAN37  sp|P62805|H4_HUMAN38  sp|P62805|H4_HUMAN39 
##                   "r"                   "l"                   "a" 
##  sp|P62805|H4_HUMAN40  sp|P62805|H4_HUMAN41  sp|P62805|H4_HUMAN42 
##                   "r"                   "r"                   "g" 
##  sp|P62805|H4_HUMAN43  sp|P62805|H4_HUMAN44  sp|P62805|H4_HUMAN45 
##                   "g"                   "v"                   "k" 
##  sp|P62805|H4_HUMAN46  sp|P62805|H4_HUMAN47  sp|P62805|H4_HUMAN48 
##                   "r"                   "i"                   "s" 
##  sp|P62805|H4_HUMAN49  sp|P62805|H4_HUMAN50  sp|P62805|H4_HUMAN51 
##                   "g"                   "l"                   "i" 
##  sp|P62805|H4_HUMAN52  sp|P62805|H4_HUMAN53  sp|P62805|H4_HUMAN54 
##                   "y"                   "e"                   "e" 
##  sp|P62805|H4_HUMAN55  sp|P62805|H4_HUMAN56  sp|P62805|H4_HUMAN57 
##                   "t"                   "r"                   "g" 
##  sp|P62805|H4_HUMAN58  sp|P62805|H4_HUMAN59  sp|P62805|H4_HUMAN60 
##                   "v"                   "l"                   "k" 
##  sp|P62805|H4_HUMAN61  sp|P62805|H4_HUMAN62  sp|P62805|H4_HUMAN63 
##                   "v"                   "f"                   "l" 
##  sp|P62805|H4_HUMAN64  sp|P62805|H4_HUMAN65  sp|P62805|H4_HUMAN66 
##                   "e"                   "n"                   "v" 
##  sp|P62805|H4_HUMAN67  sp|P62805|H4_HUMAN68  sp|P62805|H4_HUMAN69 
##                   "i"                   "r"                   "d" 
##  sp|P62805|H4_HUMAN70  sp|P62805|H4_HUMAN71  sp|P62805|H4_HUMAN72 
##                   "a"                   "v"                   "t" 
##  sp|P62805|H4_HUMAN73  sp|P62805|H4_HUMAN74  sp|P62805|H4_HUMAN75 
##                   "y"                   "t"                   "e" 
##  sp|P62805|H4_HUMAN76  sp|P62805|H4_HUMAN77  sp|P62805|H4_HUMAN78 
##                   "h"                   "a"                   "k" 
##  sp|P62805|H4_HUMAN79  sp|P62805|H4_HUMAN80  sp|P62805|H4_HUMAN81 
##                   "r"                   "k"                   "t" 
##  sp|P62805|H4_HUMAN82  sp|P62805|H4_HUMAN83  sp|P62805|H4_HUMAN84 
##                   "v"                   "t"                   "a" 
##  sp|P62805|H4_HUMAN85  sp|P62805|H4_HUMAN86  sp|P62805|H4_HUMAN87 
##                   "m"                   "d"                   "v" 
##  sp|P62805|H4_HUMAN88  sp|P62805|H4_HUMAN89  sp|P62805|H4_HUMAN90 
##                   "v"                   "y"                   "a" 
##  sp|P62805|H4_HUMAN91  sp|P62805|H4_HUMAN92  sp|P62805|H4_HUMAN93 
##                   "l"                   "k"                   "r" 
##  sp|P62805|H4_HUMAN94  sp|P62805|H4_HUMAN95  sp|P62805|H4_HUMAN96 
##                   "q"                   "g"                   "r" 
##  sp|P62805|H4_HUMAN97  sp|P62805|H4_HUMAN98  sp|P62805|H4_HUMAN99 
##                   "t"                   "l"                   "y" 
## sp|P62805|H4_HUMAN100 sp|P62805|H4_HUMAN101 sp|P62805|H4_HUMAN102 
##                   "g"                   "f"                   "g" 
## sp|P62805|H4_HUMAN103 
##                   "g"

στην προκειμένη περίπτωση η λιστα “κληροδοτεί” ονόματα στο διάνυσμα κάτι που δεν θέλουμε γιατί δημιουργεί σύγχυση. Μπορούμε να τα αφαιρέσουμε με την \(unname()\):

myseq<-unlist(fastaseqs[1])
myseq<-unname(myseq)
myseq
##   [1] "m" "s" "g" "r" "g" "k" "g" "g" "k" "g" "l" "g" "k" "g" "g" "a" "k" "r"
##  [19] "h" "r" "k" "v" "l" "r" "d" "n" "i" "q" "g" "i" "t" "k" "p" "a" "i" "r"
##  [37] "r" "l" "a" "r" "r" "g" "g" "v" "k" "r" "i" "s" "g" "l" "i" "y" "e" "e"
##  [55] "t" "r" "g" "v" "l" "k" "v" "f" "l" "e" "n" "v" "i" "r" "d" "a" "v" "t"
##  [73] "y" "t" "e" "h" "a" "k" "r" "k" "t" "v" "t" "a" "m" "d" "v" "v" "y" "a"
##  [91] "l" "k" "r" "q" "g" "r" "t" "l" "y" "g" "f" "g" "g"

που αποδίδει την αλληλουχία σε διάνυσμα με κάθε κατάλοιπο να είναι ξεχωριστό στοιχείο.

Μετατροπή διανυσμάτων σε ενιαία strings και αντιστρόφως

Κάποιες συναρτήσεις επιδρούν πάνω σε ενιαία strings. Αν θέλουμε να μετατρέψουμε το παραπάνω διάνυσμα χαρακτήρων σε ενιαίο string μπορούμε να δράσουμε με την paste, θέτοντας την παράμετρο collapse ίση με τον κενό χαρακτήρα "":

oneseq<-paste(myseq, collapse = "")
oneseq
## [1] "msgrgkggkglgkggakrhrkvlrdniqgitkpairrlarrggvkrisgliyeetrgvlkvflenvirdavtytehakrktvtamdvvyalkrqgrtlygfgg"

H αντίστροφη διαδικασία είναι κάπως ποιο πολυπλοκη. Αν δηλαδή θέλουμε να “ανοίξουμε” μια αλληλουχία από ενιαίο string σε διάνυσμα, χρησιμοποιούμε αρχικά την \(strsplit()\)

strsplit(oneseq, "")
## [[1]]
##   [1] "m" "s" "g" "r" "g" "k" "g" "g" "k" "g" "l" "g" "k" "g" "g" "a" "k" "r"
##  [19] "h" "r" "k" "v" "l" "r" "d" "n" "i" "q" "g" "i" "t" "k" "p" "a" "i" "r"
##  [37] "r" "l" "a" "r" "r" "g" "g" "v" "k" "r" "i" "s" "g" "l" "i" "y" "e" "e"
##  [55] "t" "r" "g" "v" "l" "k" "v" "f" "l" "e" "n" "v" "i" "r" "d" "a" "v" "t"
##  [73] "y" "t" "e" "h" "a" "k" "r" "k" "t" "v" "t" "a" "m" "d" "v" "v" "y" "a"
##  [91] "l" "k" "r" "q" "g" "r" "t" "l" "y" "g" "f" "g" "g"

που δημιουργεί μια λίστα, την οποία μπορούμε να πάρουμε ως διάνυσμα κατά τα γνωστά είτε με διπλή αγκύλη:

strsplit(oneseq, "")[[1]]
##   [1] "m" "s" "g" "r" "g" "k" "g" "g" "k" "g" "l" "g" "k" "g" "g" "a" "k" "r"
##  [19] "h" "r" "k" "v" "l" "r" "d" "n" "i" "q" "g" "i" "t" "k" "p" "a" "i" "r"
##  [37] "r" "l" "a" "r" "r" "g" "g" "v" "k" "r" "i" "s" "g" "l" "i" "y" "e" "e"
##  [55] "t" "r" "g" "v" "l" "k" "v" "f" "l" "e" "n" "v" "i" "r" "d" "a" "v" "t"
##  [73] "y" "t" "e" "h" "a" "k" "r" "k" "t" "v" "t" "a" "m" "d" "v" "v" "y" "a"
##  [91] "l" "k" "r" "q" "g" "r" "t" "l" "y" "g" "f" "g" "g"

ή με την \(unlist()\)

unlist(strsplit(oneseq, ""))
##   [1] "m" "s" "g" "r" "g" "k" "g" "g" "k" "g" "l" "g" "k" "g" "g" "a" "k" "r"
##  [19] "h" "r" "k" "v" "l" "r" "d" "n" "i" "q" "g" "i" "t" "k" "p" "a" "i" "r"
##  [37] "r" "l" "a" "r" "r" "g" "g" "v" "k" "r" "i" "s" "g" "l" "i" "y" "e" "e"
##  [55] "t" "r" "g" "v" "l" "k" "v" "f" "l" "e" "n" "v" "i" "r" "d" "a" "v" "t"
##  [73] "y" "t" "e" "h" "a" "k" "r" "k" "t" "v" "t" "a" "m" "d" "v" "v" "y" "a"
##  [91] "l" "k" "r" "q" "g" "r" "t" "l" "y" "g" "f" "g" "g"

Ας δούμε και μια τελευταία περίπτωση ανάλυσης. Με τη συνάρτηση που θα βρείτε εδώ μπορούμε να διαβάσουμε το ίδιο αρχείο χωρίς τα attributes.

source("~/Dropbox/Programs/R_functions/readfastafile.R")
newseqs<-readfastafile("~/Dropbox/Teaching/EAP_BNP54/Class_Winter_2020/Exercises_Fall_2020/Exercise2/histone4.fa")
str(newseqs)
## List of 10
##  $ sp|P62805|H4_HUMAN Histone H4 OS=Homo sapiens GN=HIST1H4A PE=1 SV=2                                     : chr "MSGRGKGGKGLGKGGAKRHRKVLRDNIQGITKPAIRRLARRGGVKRISGLIYEETRGVLKVFLENVIRDAVTYTEHAKRKTVTAMDVVYALKRQGRTLYGFGG"
##  $ sp|P62806|H4_MOUSE Histone H4 OS=Mus musculus GN=Hist1h4a PE=1 SV=2                                     : chr "MSGRGKGGKGLGKGGAKRHRKVLRDNIQGITKPAIRRLARRGGVKRISGLIYEETRGVLKVFLENVIRDAVTYTEHAKRKTVTAMDVVYALKRQGRTLYGFGG"
##  $ sp|P84040|H4_DROME Histone H4 OS=Drosophila melanogaster GN=His4 PE=1 SV=2                              : chr "MTGRGKGGKGLGKGGAKRHRKVLRDNIQGITKPAIRRLARRGGVKRISGLIYEETRGVLKVFLENVIRDAVTYTEHAKRKTVTAMDVVYALKRQGRTLYGFGG"
##  $ sp|P59259|H4_ARATH Histone H4 OS=Arabidopsis thaliana GN=At1g07660 PE=1 SV=2                            : chr "MSGRGKGGKGLGKGGAKRHRKVLRDNIQGITKPAIRRLARRGGVKRISGLIYEETRGVLKIFLENVIRDAVTYTEHARRKTVTAMDVVYALKRQGRTLYGFGG"
##  $ sp|P62803|H4_BOVIN Histone H4 OS=Bos taurus PE=1 SV=2                                                   : chr "MSGRGKGGKGLGKGGAKRHRKVLRDNIQGITKPAIRRLARRGGVKRISGLIYEETRGVLKVFLENVIRDAVTYTEHAKRKTVTAMDVVYALKRQGRTLYGFGG"
##  $ sp|P62784|H4_CAEEL Histone H4 OS=Caenorhabditis elegans GN=his-1 PE=1 SV=2                              : chr "MSGRGKGGKGLGKGGAKRHRKVLRDNIQGITKPAIRRLARRGGVKRISGLIYEETRGVLKVFLENVIRDAVTYCEHAKRKTVTAMDVVYALKRQGRTLYGFGG"
##  $ sp|P62801|H4_CHICK Histone H4 OS=Gallus gallus GN=H4-I PE=1 SV=2                                        : chr "MSGRGKGGKGLGKGGAKRHRKVLRDNIQGITKPAIRRLARRGGVKRISGLIYEETRGVLKVFLENVIRDAVTYTEHAKRKTVTAMDVVYALKRQGRTLYGFGG"
##  $ sp|P62804|H4_RAT Histone H4 OS=Rattus norvegicus GN=Hist1h4b PE=1 SV=2                                  : chr "MSGRGKGGKGLGKGGAKRHRKVLRDNIQGITKPAIRRLARRGGVKRISGLIYEETRGVLKVFLENVIRDAVTYTEHAKRKTVTAMDVVYALKRQGRTLYGFGG"
##  $ sp|P02309|H4_YEAST Histone H4 OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) GN=HHF1 PE=1 SV=2: chr "MSGRGKGGKGLGKGGAKRHRKILRDNIQGITKPAIRRLARRGGVKRISGLIYEEVRAVLKSFLESVIRDSVTYTEHAKRKTVTSLDVVYALKRQGRTLYGFGG"
##  $ sp|P62799|H4_XENLA Histone H4 OS=Xenopus laevis PE=1 SV=2                                               : chr "MSGRGKGGKGLGKGGAKRHRKVLRDNIQGITKPAIRRLARRGGVKRISGLIYEETRGVLKVFLENVIRDAVTYTEHAKRKTVTAMDVVYALKRQGRTLYGFGG"

Η συγκεκριμένη συνάρτηση κρατάει κάθε αλληλουχία σαν ενιαίο string. Μπορούμε να δημιουργησουμε διανύσματα για όλες τις αλληλουχίες ταυτόχρονα με τη χρήση της \(lapply()\) ως εξής:

newlistseq<-lapply(newseqs, function(x) {strsplit(x, "")})
head(newlistseq)
## $`sp|P62805|H4_HUMAN Histone H4 OS=Homo sapiens GN=HIST1H4A PE=1 SV=2`
## $`sp|P62805|H4_HUMAN Histone H4 OS=Homo sapiens GN=HIST1H4A PE=1 SV=2`[[1]]
##   [1] "M" "S" "G" "R" "G" "K" "G" "G" "K" "G" "L" "G" "K" "G" "G" "A" "K" "R"
##  [19] "H" "R" "K" "V" "L" "R" "D" "N" "I" "Q" "G" "I" "T" "K" "P" "A" "I" "R"
##  [37] "R" "L" "A" "R" "R" "G" "G" "V" "K" "R" "I" "S" "G" "L" "I" "Y" "E" "E"
##  [55] "T" "R" "G" "V" "L" "K" "V" "F" "L" "E" "N" "V" "I" "R" "D" "A" "V" "T"
##  [73] "Y" "T" "E" "H" "A" "K" "R" "K" "T" "V" "T" "A" "M" "D" "V" "V" "Y" "A"
##  [91] "L" "K" "R" "Q" "G" "R" "T" "L" "Y" "G" "F" "G" "G"
## 
## 
## $`sp|P62806|H4_MOUSE Histone H4 OS=Mus musculus GN=Hist1h4a PE=1 SV=2`
## $`sp|P62806|H4_MOUSE Histone H4 OS=Mus musculus GN=Hist1h4a PE=1 SV=2`[[1]]
##   [1] "M" "S" "G" "R" "G" "K" "G" "G" "K" "G" "L" "G" "K" "G" "G" "A" "K" "R"
##  [19] "H" "R" "K" "V" "L" "R" "D" "N" "I" "Q" "G" "I" "T" "K" "P" "A" "I" "R"
##  [37] "R" "L" "A" "R" "R" "G" "G" "V" "K" "R" "I" "S" "G" "L" "I" "Y" "E" "E"
##  [55] "T" "R" "G" "V" "L" "K" "V" "F" "L" "E" "N" "V" "I" "R" "D" "A" "V" "T"
##  [73] "Y" "T" "E" "H" "A" "K" "R" "K" "T" "V" "T" "A" "M" "D" "V" "V" "Y" "A"
##  [91] "L" "K" "R" "Q" "G" "R" "T" "L" "Y" "G" "F" "G" "G"
## 
## 
## $`sp|P84040|H4_DROME Histone H4 OS=Drosophila melanogaster GN=His4 PE=1 SV=2`
## $`sp|P84040|H4_DROME Histone H4 OS=Drosophila melanogaster GN=His4 PE=1 SV=2`[[1]]
##   [1] "M" "T" "G" "R" "G" "K" "G" "G" "K" "G" "L" "G" "K" "G" "G" "A" "K" "R"
##  [19] "H" "R" "K" "V" "L" "R" "D" "N" "I" "Q" "G" "I" "T" "K" "P" "A" "I" "R"
##  [37] "R" "L" "A" "R" "R" "G" "G" "V" "K" "R" "I" "S" "G" "L" "I" "Y" "E" "E"
##  [55] "T" "R" "G" "V" "L" "K" "V" "F" "L" "E" "N" "V" "I" "R" "D" "A" "V" "T"
##  [73] "Y" "T" "E" "H" "A" "K" "R" "K" "T" "V" "T" "A" "M" "D" "V" "V" "Y" "A"
##  [91] "L" "K" "R" "Q" "G" "R" "T" "L" "Y" "G" "F" "G" "G"
## 
## 
## $`sp|P59259|H4_ARATH Histone H4 OS=Arabidopsis thaliana GN=At1g07660 PE=1 SV=2`
## $`sp|P59259|H4_ARATH Histone H4 OS=Arabidopsis thaliana GN=At1g07660 PE=1 SV=2`[[1]]
##   [1] "M" "S" "G" "R" "G" "K" "G" "G" "K" "G" "L" "G" "K" "G" "G" "A" "K" "R"
##  [19] "H" "R" "K" "V" "L" "R" "D" "N" "I" "Q" "G" "I" "T" "K" "P" "A" "I" "R"
##  [37] "R" "L" "A" "R" "R" "G" "G" "V" "K" "R" "I" "S" "G" "L" "I" "Y" "E" "E"
##  [55] "T" "R" "G" "V" "L" "K" "I" "F" "L" "E" "N" "V" "I" "R" "D" "A" "V" "T"
##  [73] "Y" "T" "E" "H" "A" "R" "R" "K" "T" "V" "T" "A" "M" "D" "V" "V" "Y" "A"
##  [91] "L" "K" "R" "Q" "G" "R" "T" "L" "Y" "G" "F" "G" "G"
## 
## 
## $`sp|P62803|H4_BOVIN Histone H4 OS=Bos taurus PE=1 SV=2`
## $`sp|P62803|H4_BOVIN Histone H4 OS=Bos taurus PE=1 SV=2`[[1]]
##   [1] "M" "S" "G" "R" "G" "K" "G" "G" "K" "G" "L" "G" "K" "G" "G" "A" "K" "R"
##  [19] "H" "R" "K" "V" "L" "R" "D" "N" "I" "Q" "G" "I" "T" "K" "P" "A" "I" "R"
##  [37] "R" "L" "A" "R" "R" "G" "G" "V" "K" "R" "I" "S" "G" "L" "I" "Y" "E" "E"
##  [55] "T" "R" "G" "V" "L" "K" "V" "F" "L" "E" "N" "V" "I" "R" "D" "A" "V" "T"
##  [73] "Y" "T" "E" "H" "A" "K" "R" "K" "T" "V" "T" "A" "M" "D" "V" "V" "Y" "A"
##  [91] "L" "K" "R" "Q" "G" "R" "T" "L" "Y" "G" "F" "G" "G"
## 
## 
## $`sp|P62784|H4_CAEEL Histone H4 OS=Caenorhabditis elegans GN=his-1 PE=1 SV=2`
## $`sp|P62784|H4_CAEEL Histone H4 OS=Caenorhabditis elegans GN=his-1 PE=1 SV=2`[[1]]
##   [1] "M" "S" "G" "R" "G" "K" "G" "G" "K" "G" "L" "G" "K" "G" "G" "A" "K" "R"
##  [19] "H" "R" "K" "V" "L" "R" "D" "N" "I" "Q" "G" "I" "T" "K" "P" "A" "I" "R"
##  [37] "R" "L" "A" "R" "R" "G" "G" "V" "K" "R" "I" "S" "G" "L" "I" "Y" "E" "E"
##  [55] "T" "R" "G" "V" "L" "K" "V" "F" "L" "E" "N" "V" "I" "R" "D" "A" "V" "T"
##  [73] "Y" "C" "E" "H" "A" "K" "R" "K" "T" "V" "T" "A" "M" "D" "V" "V" "Y" "A"
##  [91] "L" "K" "R" "Q" "G" "R" "T" "L" "Y" "G" "F" "G" "G"

Ανάλυση καταλοίπων σε fasta

Βασικά χρησιμοποιούμε τις γενικές συναρτήσεις ανάλυσης σειρών χαρακτήρων όπως η grep(), gsub() κλπ. Αν μας ενδιαφέρει να βρούμε το πλήθος εμφανίσεων ενός χαρακτήρα σε ένα string μια πρώτη λύση είναι: α) να τον αντικαταστήσουμε με τον κενό χαρακτήρα παντού στο string και β) να υπολογίσουμε τη διαφορά στα μήκη των δύο strings. To μήκος ενός string δίνεται από την \(nchar()\):

nchar(oneseq)
## [1] 103

η αντικατάσταση π.χ. του καταλοίπου “a” παντού στο string γίνεται με την \(gsub()\):

anotherseq<-gsub("a", "", oneseq)
nchar(anotherseq)
## [1] 96

H διαφορά των δύο είναι προφανώς το πλήθος των “a” στην αρχική αλληλουχία:

nchar(oneseq)-nchar(anotherseq)
## [1] 7

Mια πιο δομημένη λύση που υπολογίζει όλα τα πλήθη των καταλοίπων απευθείας είναι αρχικά να μετατραπεί η αλληλουχία σε διάνυσμα όπως έχουμε δει:

vectorseq<-unlist(strsplit(oneseq, ""))

και να εφαρμόσουμε την \(table()\) στο διάνυσμα που προκύπτει:

table(vectorseq)
## vectorseq
##  a  d  e  f  g  h  i  k  l  m  n  p  q  r  s  t  v  y 
##  7  3  4  2 17  2  6 11  8  2  2  1  2 14  2  7  9  4

H \(table()\) δημιουργεί ένα αντικείμενο τύπου table που μπορούμε να προσπελάσουμε με τη χρήση της συνάρτησης names για να βρούμε την τιμή για ένα συγκεκριμένο κατάλοιπο:

tab<-table(vectorseq)
tab[which(names(tab)=="m")]
## m 
## 2

Βρόχοι επανάληψης for/while

Παρότι η χρήση βρόχων επανάληψης γενικώς αποθαρρύνεται στην R, προς όφελος διανυσματικών πράξεων, εντούτοις υποστηρίζονται κανονικά με σχετικά απλή σύνταξη.

Η γενική σύνταξη είναι η συνθήκη επανάληψης σε παρένθεση και οι πράξεις μέσα σε άγκιστρα. Η διαφορά μεταξύ των for/while είναι ότι για τους πρώτους πρέπει να δοθεί αναλυτικά ο αριθμός των επαναλήψεων με αρχή και τέλος:

x<-seq(1:10)
y<-0
for(i in -1:10){
  y[i]<-x[i]**2
}
plot(x,y, type="o")

ενώ για τους while βρόχους χρειάζεται να έχει ήδη δηλωθεί μια συνθήκη έναρξης, η οποία να επικαιροποιείται μέσα στα άγκιστρα

x<-seq(1:10)
y<-0
i<-1
while(i<=10){
  y[i]<-x[i]**2
  i<-i+1
}
plot(x,y, type="o")