root/trunk/ProjectFortress/demos/BirdCount0.fss

(*******************************************************************************
Copyright 2009 Michael Zody and Sun Microsystems, Inc. 
All rights reserved.

Sun Microsystems, Inc. is the Copyright owner of the Fortress programming language software, and Michael Zody is the developer of the algorithm which this software implements and the Copyright owner of the software implementation of the algorithm, to which Sun Microsystems, Inc. has a perpetual, irrevocable, royalty free right and license to use and distribute. 

Use is subject to license terms accompanying the software.

 ******************************************************************************)

component BirdCount0
import File.{...}
import FileSupport.{...}
import FlatString.{...}
import List.{...}
import Map.{...}
import Set.{...}
import System.{getProperty}

export Executable

N:ZZ32 = 4

conversion:Map[\String,ZZ32\] = 
    {[\String,ZZ32\]
         "AA"   |->0  ,
         "CA"   |->1  ,
         "GA"   |->2  ,
         "TA"   |->3  ,
         "AC"   |->1  ,
         "CC"   |->0  ,
         "GC"   |->3  ,
         "TC"   |->2  ,
         "AG"   |->2  ,
         "CG"   |->3  ,
         "GG"   |->0  ,
         "TG"   |->1  ,
         "AT"   |->3  ,
         "CT"   |->2  ,
         "GT"   |->1  ,
         "TT"   |->0  ,
         "AN"  |->N  , 
         "CN"  |->N  , 
         "GN"  |->N  , 
         "TN"  |->N  , 
         "NA"  |->N  , 
         "NC"  |->N  , 
         "NG"  |->N  , 
         "NT"  |->N  , 
         "NN"  |->N  
      }

inverseConversion:Map[\String,Char\] =
    {[\String, Char\]
        "A0" |-> 'A',
        "A1" |-> 'C',
        "A2" |-> 'G',
        "A3" |-> 'T',
        "C0" |-> 'C',
        "C1" |-> 'A',
        "C2" |-> 'T',
        "C3" |-> 'G',
        "G0" |-> 'G',
        "G1" |-> 'T',
        "G2" |-> 'A',
        "G3" |-> 'C',
        "T0" |-> 'T',
        "T1" |-> 'G',
        "T2" |-> 'C',
        "T3" |-> 'A'
    }
       

validTransitions:Set[\String\] = {[\String\] 
    "0011","0022","0033","0110","0123","0132","0213","0220","0231","0312","0321","0330",
    "1001","1023","1032","1100","1122","1133","1203","1221","1230","1302","1320","1331",
    "2002","2013","2031","2103","2112","2130","2200","2211","2233","2301","2310","2332",
    "3003","3012","3021","3102","3113","3120","3201","3210","3223","3300","3311","3322"
}

referenceChickenFile:String = getProperty("fortress.automhome",".") || "/ProjectFortress/demos/chickenData/chr10.fa.head" 
sampleFilesDir:String = getProperty("fortress.automhome",".") || "/ProjectFortress/demos/chickenData/test_data"
sampleFilesFileName:String = "chr10.csfasta.ma.sorted.head"

sampleChickens:List[\String\] = <|[\String\]
    "10a_LA_Frag35_20080704_white_leghorn_A",
    "10b_LA_Frag35_20080704_white_leghorn_B"
(*    "17_LA_Frag35_20080909_high_grow_line",
    "18_LA_Frag35_20080909_low_grow_line",
    "17.1_LA_Frag35_20080925_high_grow_line",
    "18.1_LA_Frag35_20080925_low_grow_line",
    "16_LA_Frag35_20080829_broiler",
    "22_LA_Frag35_20081016_RJF",
    "16.1_LA_Frag35_20080925_broiler",
    "22.1_LA_Frag35_20081114_RJF" *)
|>

convert(c:(Char, Char)):String = do
    (c1:Char,c2:Char) = c
   conversion.member("" || c1 || c2).get()
end

object Pairs[\T\](g: Generator[\T\]) extends Generator[\(T,T)\]
    generate[\R\](red: Reduction[\R\], m:(T,T)->R): R =
    if (l,v,r) <- g.generate[\Maybe[\(T,R,T)\]\](PairReduction[\T,R\](red,m),
                             fn (t:T) => Just[\(T,R,T)\](t, red.empty(), t))
        then v else red.empty() end
end

object PairReduction[\T,R\](red: Reduction[\R\], m:(T,T)->R)
        extends Reduction[\Maybe[\(T,R,T)\]\]
    empty():Maybe[\(T,R,T)\] = Nothing[\(T,R,T)\]
    join(left:Maybe[\(T,R,T)\], right:Maybe[\(T,R,T)\]) : Maybe[\(T,R,T)\]=
    if (l, v_l, m_l) <- left then
        if (m_r, v_r, r) <- right then
            Just[\(T,R,T)\](l, red.join(red.join(v_l, m(m_l, m_r)), v_r), r)
        else left end
    else right end
end 


trait dna comprises {dnaSequence, EmptySequence}
       getter colorSequence(): String
       opr OPLUS(self, other: dna) : dna
end


(* left and right are ACGT encoded, middle is color encoded *)

object dnaSequence(left:Char, middle:String, right:Char) extends dna
   getter colorSequence():String = middle
   opr OPLUS(self, other:dnaSequence):dnaSequence =
       dnaSequence(left, middle || convert(right, other.left) || other.middle(), other.right)
   opr OPLUS(self, other:EmptySequence) = self
end

object EmptySequence() extends dna
    getter colorSequence():String = ""
    opr OPLUS(self, other:dna) = other
end

object dnaSequenceReduction extends { MonoidReduction[\dna\],
                               ReductionWithZeroes[\dna,dna\] }
    getter asString() = "dnaSequenceReduction"
    empty(): dna = EmptySequence()
    join(a: dna, b: dna): dna = a OPLUS b
end

opr BIG OPLUS[\T\](): BigReduction[\dna, dna\] =
   BigReduction[\dna,dna\](dnaSequenceReduction)

opr BIG OPLUS[\T\](g:Generator[\dna\]) = 
   __bigOperatorSugar[\dna,dna,dna,dna\](BIG OPLUS[\T\](), g)


processReferenceFileLine(input:String):dna = do
    if input.get(0) = '>' 
        then EmptySequence()
    else 
        size = |input|
        println input " has " size " characters "  
        dnaSequence(input.get(0), BIG || [p<- Pairs(input)](convert(p)), input.get(size-1)) 
    end
end

processReferenceFile(name:String):String=do
   var rs:FileReadStream = FileReadStream(name)
   var res:dna  = BIG OPLUS [l<-rs.lines()] processReferenceFileLine(l)
   colors = res.colorSequence()
   println "result has " |colors| " characters which should be 1 less than input ACGT colors."
   colors
end

readReferenceFile(name:String):String=do
   var rs:FileReadStream = FileReadStream(name)
   var res:String = BIG || [l<-rs.lines()] l
   res
end

object Snip(header:String, sequence:String, name:String, pos:ZZ32, length:ZZ32, seqend:ZZ32, 
            referenceChickenColorsSnip:String, sampleChickenColorsSnip:String,
            referenceChickenAGCTSnip:String) 
    getter header():String = header
    getter sequence():String = sequence
    getter name():String = name
    getter pos():ZZ32 = pos
    getter length():ZZ32 = length
    getter seqend():ZZ32 = seqend
    getter referenceChickenColorsSnip():String = referenceChickenColorsSnip
    getter sampleChickenColorsSnip():String = sampleChickenColorsSnip
    getter referenceChickenAGCTSnip():String = referenceChickenAGCTSnip    
    getter asString():String = "\\n{Snip:" || name || "\\n" || referenceChickenColorsSnip || "\\n" sampleChickenColorsSnip || "\\n" || referenceChickenAGCTSnip 

    sampleChickenAGCTSnip():String = do

        l:ZZ32 = |sampleChickenColorsSnip|
        var i:ZZ32 := 0
        var result:String := ""
        var current:Char := sampleChickenColorsSnip.get(0)
        var c:Char = referenceChickenAGCTSnip.get(0)

        while (i < l) do
          result := result || c
          current := sampleChickenColorsSnip.get(i)
          temp:String = c.asString() || current.asString()
          c:= inverseConversion.member(temp).get
          i:= i+1
        end
        result := result || c
        result
    end

(*
       BIG || [(x,y) <- referenceChickenAGCTSnip.zip[\Char\](referenceChickenColorsSnip.javaSubstr(0,1) || sampleChickenColorsSnip)] inverseConversion.member(x.asString() || y.asString()).get()

*)

    printSnip():() = do
        println("Snip:" name)
        println(" " referenceChickenColorsSnip)
        println(" " sampleChickenColorsSnip)
        println(referenceChickenAGCTSnip)
(*        println(sampleChickenAGCTSnip())  *)
    end


end

reverse(sequence:String):String = BIG || [c<-sequence.reverse()] c

negativeOrientation(header:String):Boolean =  do
   pattern:String = ">\\w*,\\d*_-.*"
   header.javaRegExpMatches(pattern)
end

getAdjustedLocation(negativeOri:Boolean, location:String, length:ZZ32):ZZ32 = do
   if (negativeOri) then
      strToInt(location.asFlatString().javaSubstr(1), 10) - length + 1        
   else 
      strToInt(location) - 1
   end
end

getReferenceChickenColorsSnip(adjustedLoc:ZZ32, seqEnd:ZZ32):FlatString = do
   referenceChickenColors[adjustedLoc:seqEnd+1].asFlatString()
end

getAdjustedSequence(negativeOri:Boolean, sequence:String, refChickenColorsSnip:FlatString):FlatString = do
    temp:FlatString = sequence.asFlatString().javaSubstr(1)
    tempLength:ZZ32 = |temp|
    if (negativeOri) then
        ((reverse(temp).asFlatString().javaSubstr(0, tempLength - 1 )).asFlatString() || refChickenColorsSnip.asFlatString().javaSubstr(tempLength - 1, tempLength).asFlatString()).asFlatString()
    else
        (refChickenColorsSnip.asFlatString().javaSubstr(0,1) || temp.asFlatString().javaSubstr(1).asFlatString()).asFlatString()
    end
end


readASnip(r:ReadStream):Snip = do
   header:String = r.uncheckedReadLine()
   var sequence:String = r.uncheckedReadLine()

   if sequence = "" then Snip(header, sequence,"", 0, 0, 0, "","","") else
      name:String = header.asFlatString().javaRegExpSplit(",",0)
      loc:String = header.asFlatString().javaRegExpSplit(",",1).asFlatString().javaRegExpSplit("_",1).asFlatString().javaRegExpSplit("\\.",0)
      length:ZZ32 = sequence.asFlatString().size() - 1
      negativeOri:Boolean = negativeOrientation(header)
      adjustedLocation:ZZ32 = getAdjustedLocation(negativeOri, loc, length)
      sequenceEnd:ZZ32 = adjustedLocation + length - 1
      referenceChickenColorsSnip:FlatString = getReferenceChickenColorsSnip(adjustedLocation, sequenceEnd)
      adjustedSequence:FlatString := getAdjustedSequence(negativeOri, sequence, referenceChickenColorsSnip)
      referenceChickenAGCTSnip:String = referenceChickenAGCT.asFlatString().javaSubstr(adjustedLocation+6, sequenceEnd + 6 + 2)
      Snip(header, adjustedSequence, name, adjustedLocation, length, sequenceEnd, referenceChickenColorsSnip, adjustedSequence, referenceChickenAGCTSnip)
   end
end


ChunkSize:ZZ32 = 10000


(* Given two samples, this creates a list of the positions where they differ *)

object Event(s:Snip, diffs:String, startPos:ZZ32, endPos:ZZ32)
    getter snip():Snip = s
    getter startPos():ZZ32 = startPos
    getter endPos():ZZ32 = endPos
    getter diffs():String = diffs
    getter asString():String = "{Event:" || s.header() || ":" startPos || "," || endPos || "}"
    printEvent():() = do
       println("Event:" )
       s.printSnip()
       println(" " diffs)
       var i:ZZ32 := 0
       var c:Char := s.referenceChickenAGCTSnip().get(startPos)
       var result:String :=  c.asString()
       while (i < endPos) do
          next:Char = inverseConversion.member(c.asString() || s.sampleChickenColorsSnip.get(i+startPos).asString()).get()
          result := result || next
          c := next
          i := i+1
       end
       println(s.referenceChickenAGCTSnip.javaSubstr(startPos, endPos+startPos+1))
       println(result)
    end
end

(* This code is a little obtuse.  The idea is that if the sample differs from the reference 
   chicken in more than one consecutive location it may be an interesting mutation.  We look
   for these potential interesting mutations and return a list of them.
 *)

makeTuple(loc:ZZ32, diffs:String):(ZZ32,ZZ32) = do
    temp:String = diffs.asFlatString().javaSubstr(loc).javaRegExpSplit("\\.",0)
    (loc, |temp|)
end

makeEvent(s:Snip, diffs:String, loc:ZZ32):Event = do
    temp:String = diffs.asFlatString().javaSubstr(loc).javaRegExpSplit("\\.",0)
    Event(s, diffs, loc, |temp|)
end   

isInterestingPosition(loc:ZZ32, diffs:String):Boolean = do
   if loc = 0 OR: diffs.get(loc) = '.'  OR: loc = (|diffs| - 1) OR: (diffs.get(loc - 1) = 'X') OR: (|diffs.asFlatString().javaSubstr(loc).javaRegExpSplit("\\.",0)| < 2) then
      false
   else true
   end
end

isValidPairTransition(r:String, s:String) = do
    (r || s) IN validTransitions
end

isValidTransition(reference:String, sample:String) = do
    println("Reference = " reference)
    println("sample = " sample)
    var i:ZZ32 := 0
    var result:Boolean := true
    while (i < |reference| - 1) do
        r = reference.javaSubstr(i, 2)
        s = sample.javaSubstr(i, 2)
        result := result AND isValidPairTransition(r,s)
        i := i + 1
    end
    result
end

isValidTransition(s:Snip, loc:ZZ32, diffs:String):Boolean = do
   len:ZZ32 = |diffs.asFlatString().javaSubstr(loc).javaRegExpSplit("\\.",0)|
   sample:String = s.sampleChickenColorsSnip().javaSubstr(loc, loc+len)
   reference:String = s.referenceChickenColorsSnip().javaSubstr(loc, loc+len)
   isValidTransition(reference, sample)
end

isInteresting(s:Snip, loc:ZZ32, diffs:String):Boolean = do
   isInterestingPosition(loc, diffs) AND: isValidTransition(s, loc,diffs)
end




sampleCompare(ref:String, sample:String):String = do
   BIG || [(x,y) <- ref.zip[\Char\](sample)] (if (x =/= y) then "X" else "." end)
end

(* Given two samples this creates events *)

eventGenerator(s:Snip):Map[\ZZ32, List[\(ZZ32,ZZ32)\]\] = do
    diffs = sampleCompare(s.referenceChickenColorsSnip(),s.sampleChickenColorsSnip())
    stop = |diffs| - 1
    var database:Map[\ZZ32, List[\Event\]\] := {[\ZZ32, List[\Event\]\]}
    for x <- 0:stop do
        if isInteresting(s,x,diffs) then
           database := database.add(x + s.pos(), <| makeEvent(s, diffs, x) |>)
        end
    end
    database
end

combine(k:ZZ32, val1:List[\Event\], val2:List[\Event\]):List[\Event\] = val1 || val2

processSampleChicken(name:String):Map[\ZZ32, List[\Event\]\] = do
    max:ZZ32 = referenceChickenColors.size()
    var rs:FileReadStream = 
        FileReadStream(sampleFilesDir || "/" || name || "/" || sampleFilesFileName)
    fg:FileGenerator[\Snip\]  = FileGenerator[\Snip\](rs,60,readASnip)
    var database:Map[\ZZ32, List[\Event\]\] := {[\ZZ32, List[\Event\]\]}

    for s <- fg do
       database := database.union(combine, eventGenerator(s))
    end
    database
end

var referenceChickenColors:String := ""
var referenceChickenAGCT:String := ""

run() = do
    referenceChickenAGCT   := readReferenceFile(referenceChickenFile)
    referenceChickenColors := processReferenceFile(referenceChickenFile)
    var database:Map[\ZZ32, List[\Event\]\] := {[\ZZ32, List[\Event\]\]} 
    
    for sampleChicken <- sampleChickens do
        database := database.union(combine,processSampleChicken(sampleChicken))
    end
    println("Summary1")
    for (key,entryList) <- seq(database) do
       println(key)
       for entry <- seq(entryList) do
          entry.printEvent()
       end
    end

end  
end