Fasta sequences

Οι βιολογικές αλληλουχίες (DNA, RNA, ή πρωτείνες), αποθηκευονται συνήθως σε ένα text format που ονομάζεται FASTA. Το χαρακτηριστικό ειναι το εξής:

>ΟΝΟΜΑ_ΑΛΛΗΛΟΥΧΙΑΣ
ACCCATATA

όπου ΟΝΟΜΑ_ΑΛΛΗΛΟΥΧΙΑΣ, παριστάνει το όνομα που έχει δοθεί στην αλληλουχία που ακολουθεί. Για παράδειγμα, αν αναφερόμαστε στην αλληλουχία του γονιδίου BDNF στον άνθρωπο, ένα πιθανό όνομα θα μπορούσε να ειναι BDNF_HUMAN. Το σημαντικό ειναι ο χαρακτήρας > που υπάρχει στην αρχή της γραμμής του ονόματος.

Τι είδους αλληλουχίες μπορείτε να βρείτε ως FASTA?

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

Από που μπορείτε να κατεβάσετε FASTA αλληλουχίες

ή για παράδειγμα:

http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/upstream1000.fa.gz

Το πρώτο αρχείο περιέχει στοιχισμένα σύνολα ομόλογων γονιδίων για 100 οργανισμών (και τον άνθρωπο). Το δεύτερο αρχείο περιέχει τα γονίδια του ανθρώπου.

Εμείς μπορούμε να κατεβάσουμε τα δεδομένα ρυθμιστικών περιοχών 1000 βάσεις upstream. Άρα το αρχείο:

http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/upstream1000.fa.gz

Τι είδους χαρακτήρες μπορεί να βρείτε σε ένα αρχείο FASTA?

Αν μιλάμε για DNA:

The DNA alphabet

The DNA alphabet

Πως σαρώνουμε αλληλουχίες για να ψάξουμε που υπάρχει τα μοτίβα;

Η ιδέα ειναι οτι σε κάθε διαδοχική βάση του FASTA περνάμε το PWM και υπολογίζουμε ένα σκορ. Πώς μπορούμε να σκοράρουμε;

Όμως πως αποφασίζουμε αν κάποια σκορς ειναι υψηλά αρκετά;

Η ιδέα είναι η εξής. Αν ειχαμε τυχαίες ακολουθίες που δεν περιμένουμε να έχουν το μοτίβο, τοτε τι σκορς θα βρίσκαμε;Από εκεί μπορούμε να αποφασίσουμε αν το σκορ στην πραγματική αλληλουχία ειναι αρκετά υψηλά.

Πως θα διαβάσουμε FASTA files

Τα FASTA αρχεία ειναι αρκετά απλά στην δομή. Στην πιο απλή τους μορφή, αποτελούνται από εγγραφες ονομάτων αλληλουχιών, και αλληλουχιών, καταλαμβάνοντας διαδοχικά δύο γραμμές. Αυτό μας δίνει την δυνατότητα να το διαβάσουμε ως read.table, παίρνοντας στην συνέχεια μόνο τις ζυγές γραμμές που περιέχουν τις αλληλουχίες ή χρησιμοποιώντας την συνάρτηση grep, που θα μας δώσει τις γραμμές εκείνες που δεν περιέχουν ονόματα

a = read.table("5utr.fa")
linesWithSeq = grep(pattern="^>", x = a[,1], invert = TRUE)
a[linesWithSeq,]

Έχουμε όμως δυνητικά το εξής πρόβλημα. Αν το αρχείο ειναι πολύ μεγάλο, που συχνά είναι κάτι που συμβαίνει, τότε δεν θα χωράει στην μνήμη του υπολογιστή. Δηλαδή η μεταβλητή a θα καταλαμβάνει πολύ μεγάλο χώρο στην κύρια μνήμη του υπολογιστή, και τότε η εκτέλεση πράξεων ή θα καθυστερήσει πολύ ή αναγκαστικα θα τερματιστεί. Επομένως, θα πρέπει να σκεφτούμε διαφορετικά. Η βασική λύση συνοψίζεται στο εξής: Διαβάζουμε και αναλύουμε το αρχείο γραμμή γραμμή. Έτσι στην μνήμη θα υπάρχει μόνο η τρέχουσα γραμμή. Ο παρακάτω κώδικας ειναι βοηθητικός γι αυτό.

processFile = function(filepath) {
  con = file(filepath, "r")
  while ( TRUE ) {
    line = readLines(con, n = 1)
    if ( length(line) == 0 ) {
      break
    }
    print(line)
  }

  close(con)
}

Με την συνάρτηση αυτή, ορίζουμε ως filepath τη διευθυνση του αρχείου που θέλουμε να διαβάζουμε. Στην συνέχεια ανοίγουμε το αρχείο με την συνάρτηση file, η οποία ειναι αντίστοιχη της open, σε άλλες γλώσσες προγραμματισμού. Στην συνέχεια, χρησιμοποιώντας μια while loop (μαζί με break), διαβάζουμε το αρχείο εφόσον η γραμμή που διαβάζουμε έχει κάποιο μήκος. Στο παραπάνω παράδειγμα απλά τυπώνουμε την γραμμή, όμως μπορούμε να την διαχειριστούμε όπως θέλουμε. Τώρα σκεφτόμαστε ως εξής: Θα πρέπει να φτιάξω μια συνάρτηση που αφενός χωρίζει την FASTA γραμμή και φτιάχνει ένα vector, και στην συνέχεια να σκανάρει την αλληλουχία σε τμήματα διαδοχικά όπως φαίνεται στην παρακάτω εικόνα

Scanning a fasta sequence for a motif

processFile = function(filepath, pwm) {
  con = file(filepath, "r")
  while ( TRUE ) {
    line = readLines(con, n = 1)
    if ( length(line) == 0 ) {
      break
    }
   isSequence = grep(pattern="^>", x=line, invert = TRUE)
   
   if(length(isSequence) > 0){
      ## we split the line and we create a vector 
      seq.vec = unlist(strsplit(line, split=""))
      ## code for analyzing sequence
      ### 
      ## It is necessary to go up to the "length_vec - length motif + 1"
      for(i in 1:(length(seq.vec)-ncol(pwm)+1)){
        # This function evaluates the score for a substring of the fasta
        ## seq.vec
        ## at position i
        ## for the position weight matrix pwm
        ## note that pwm consists of some columns. For this number of columns the score will be evaluated.
        getScore(seq.vec, i, pwm)
      }
   }
  }

  close(con)
}



processFile("5utr.fa", pwm=pwm)

Βέβαια, πρέπει να συντάξουμε πρώτα την συνάρτηση getScore πυ θα παίρνει σαν όρισμα το vector με το fasta, την θέση εκτίμησης i, και το pwm.

alphabet=c("A", "C", "G", "T")

getPwm = function(motifFile){
  a = read.table(motifFile, colClasses="character")
  for(i in 1:nrow(a)){
    v = unlist(strsplit(a[i,], ""))
    ## initialize the pwm at first iteration only
    if(i == 1){
      pwm = matrix(0, nrow=length(alphabet), ncol=length(v))
      rownames(pwm) = alphabet
    }
    for(j in 1:ncol(pwm)){
      letter = v[j] 
      pwm[letter, j] = pwm[letter, j] + 1
    }
  }
  pwm = t(t(pwm)/colSums(pwm))
  pwm = log2((pwm+1e-6)/0.25)
  return(pwm)
}

pwm = getPwm("pwm.txt")

getScore  = function(v, i, pwm){
  score = 0
  maxScore = 0
  minScore = 0
  for(j in 1:ncol(pwm)){
    MAX = max(pwm[,j])
    MIN = min(pwm[,j])
    maxScore = maxScore + MAX
    minScore = minScore + MIN
  }
  for(j in 1:ncol(pwm)){
    letter = v[i+j-1]
    score = score + pwm[letter, j] 
  }
  score = (score - minScore)/(maxScore - minScore)
  return(score)
}


processFile = function(filepath, pwm) {
  seq = 1
   seq.scores = list()
  con = file(filepath, "r")
  while ( TRUE ) {
    line = readLines(con, n = 1)
    if ( length(line) == 0 ) {
      break
    }
   isSequence = grep(pattern="^>", x=line, invert = TRUE)
   
  
   if(length(isSequence) > 0){
     
     ## initialize Scores for each line
     lineScores = c()
      ## we split the line and we create a vector 
      seq.vec = unlist(strsplit(line, split=""))
      ## code for analyzing sequence
      ### 
      ## It is necessary to go up to the "length_vec - length motif + 1"
      for(i in 1:(length(seq.vec)-ncol(pwm)+1)){
        # This function evaluates the score for a substring of the fasta
        ## seq.vec
        ## at position i
        ## for the position weight matrix pwm
        ## note that pwm consists of some columns. For this number of columns the score will be evaluated.
        score = getScore(seq.vec, i, pwm)
        lineScores = c(lineScores, score)
      }
     ## print(lineScores)
      max.score = max(lineScores)
     ## print(max.score)
      seq.scores[[seq]] = max.score
      ## increment the counter
      seq = seq + 1
   }
  }
  close(con)
  return(seq.scores)
}

seqScores = processFile("5utr.fa", pwm=pwm)
randomScores = processFile("randomseq.fa", pwm=pwm)

#plot(1:length(seqScores), unlist(seqScores))
#plot(1:length(randomScores), unlist(randomScores))


{plot(1:length(seqScores), unlist(seqScores), ylim=c(0.8, 1))
points(1:length(randomScores), unlist(randomScores), col="red")
abline(h=mean(unlist(seqScores)), col="black")
abline(h=mean(unlist(randomScores)), col="red")
legend("topleft", legend=c("sequences", "random sequences"), col=c("black", "red"), pch=19)}