github.com/vertgenlab/gonomics@v1.0.0/genomeGraph/seed.go (about)

     1  package genomeGraph
     2  
     3  import (
     4  	"log"
     5  	"sort"
     6  
     7  	"github.com/vertgenlab/gonomics/fastq"
     8  )
     9  
    10  func extendCurrSeed(seed *SeedDev, gg *GenomeGraph, read fastq.Fastq, left bool, right bool) {
    11  	var newTStart, newQStart, newTEnd, newQEnd int32 = int32(seed.TargetStart) - 1, int32(seed.QueryStart) - 1, int32(seed.TargetStart + seed.Length), int32(seed.QueryStart + seed.Length)
    12  	//check to see if beginning is at index zero, if so do something like SeedDev.Prev
    13  	//if newStart < 0
    14  	if left {
    15  		for ; newTStart >= 0 && newQStart >= 0 && (gg.Nodes[seed.TargetId].Seq[newTStart] == read.Seq[newQStart]); newTStart, newQStart = newTStart-1, newQStart-1 {
    16  			seed.TargetStart = uint32(newTStart)
    17  			seed.QueryStart = uint32(newQStart)
    18  			seed.Length++
    19  		}
    20  	}
    21  	if right {
    22  		for ; int(newTEnd) < len(gg.Nodes[seed.TargetId].Seq) && int(newQEnd) < len(read.Seq) && (gg.Nodes[seed.TargetId].Seq[newTEnd] == read.Seq[newQEnd]); newTEnd, newQEnd = newTEnd+1, newQEnd+1 {
    23  			seed.Length++
    24  		}
    25  	}
    26  }
    27  
    28  func toTheRight(seed *SeedDev, gg *GenomeGraph, read fastq.Fastq) []*SeedDev {
    29  	//log.Printf("Depth of call is: %d for seed: %d", depth, seed.TargetId)
    30  	var answer []*SeedDev
    31  	extendCurrSeed(seed, gg, read, false, true)
    32  	var newTEnd, newQEnd int32 = int32(seed.TargetStart + seed.Length), int32(seed.QueryStart + seed.Length)
    33  	if (int(newTEnd) >= len(gg.Nodes[seed.TargetId].Seq) && int(newQEnd) < len(read.Seq)) && (len(gg.Nodes[seed.TargetId].Next) > 0) {
    34  		var newTStart int32 = 0
    35  		var newQStart int32 = newQEnd
    36  		var edgeSeeds []*SeedDev
    37  		var e int
    38  		for _, next := range gg.Nodes[seed.TargetId].Next {
    39  			//log.Printf("Number of nodes to check %d\n", len(gg.Nodes[seed.TargetId].Next))
    40  			if len(next.Dest.Seq) > 0 {
    41  				if next.Dest.Seq[newTStart] == read.Seq[newQStart] {
    42  					nextSeed := &SeedDev{TargetId: next.Dest.Id, TargetStart: uint32(newTStart), QueryStart: uint32(newQStart), Length: 1, PosStrand: seed.PosStrand, NextPart: nil}
    43  					edgeSeeds = toTheRight(nextSeed, gg, read)
    44  					for e = 0; e < len(edgeSeeds); e++ {
    45  						currSeed := &SeedDev{TargetId: seed.TargetId, TargetStart: seed.TargetStart, QueryStart: seed.QueryStart, Length: seed.Length, PosStrand: seed.PosStrand, NextPart: edgeSeeds[e]}
    46  						answer = append(answer, currSeed)
    47  					}
    48  				}
    49  			}
    50  		}
    51  	} else {
    52  		answer = append(answer, seed)
    53  	}
    54  	return answer
    55  }
    56  
    57  func toTheLeft(seed *SeedDev, gg *GenomeGraph, read fastq.Fastq) []*SeedDev {
    58  	var answer []*SeedDev
    59  	extendCurrSeed(seed, gg, read, true, false)
    60  	//var newTStart, newQStart int32 = int32(seed.TargetStart) - 1, int32(seed.QueryStart) - 1
    61  	if (seed.TargetStart <= 0 && seed.QueryStart > 0) && (len(gg.Nodes[seed.TargetId].Prev) > 0) {
    62  		var depthSeeds []*SeedDev
    63  		var prevLeft int
    64  		for _, prev := range gg.Nodes[seed.TargetId].Prev {
    65  			if len(prev.Dest.Seq) > 0 {
    66  				var newTStart int32 = int32(len(prev.Dest.Seq)) - 1
    67  				var newQStart int32 = int32(seed.QueryStart) - 1
    68  				if prev.Dest.Seq[newTStart] == read.Seq[newQStart] {
    69  					prevSeed := &SeedDev{TargetId: prev.Dest.Id, TargetStart: uint32(newTStart), QueryStart: uint32(newQStart), Length: 1, PosStrand: seed.PosStrand, NextPart: seed}
    70  					depthSeeds = toTheLeft(prevSeed, gg, read)
    71  					for prevLeft = 0; prevLeft < len(depthSeeds); prevLeft++ {
    72  						answer = append(answer, depthSeeds[prevLeft])
    73  					}
    74  				}
    75  			}
    76  		}
    77  	} else {
    78  		answer = append(answer, seed)
    79  	}
    80  	return answer
    81  }
    82  
    83  func extendSeedTogether(seed *SeedDev, gg *GenomeGraph, read fastq.Fastq) []*SeedDev {
    84  	var answer []*SeedDev
    85  	rightGraph := toTheRight(seed, gg, read)
    86  
    87  	for rSeeds := 0; rSeeds < len(rightGraph); rSeeds++ {
    88  		answer = append(answer, toTheLeft(rightGraph[rSeeds], gg, read)...)
    89  	}
    90  	return answer
    91  }
    92  
    93  func getLastPart(a *SeedDev) *SeedDev {
    94  	for ; a.NextPart != nil; a = a.NextPart {
    95  	}
    96  	return a
    97  }
    98  
    99  func CompareBlastScore(a *SeedDev, b *SeedDev, read fastq.Fastq, scoreMatrix [][]int64) int {
   100  	if BlastSeed(a, read, scoreMatrix) == BlastSeed(b, read, scoreMatrix) {
   101  		return 0
   102  	} else if BlastSeed(a, read, scoreMatrix) < BlastSeed(b, read, scoreMatrix) {
   103  		return -1
   104  	} else if BlastSeed(a, read, scoreMatrix) > BlastSeed(b, read, scoreMatrix) {
   105  		return 1
   106  	} else {
   107  		log.Fatalf("Error: SeedDev len compare failed on:%d %d %d, %d %d %d\n", a.TargetId, a.TargetStart, a.Length, b.TargetId, b.TargetStart, b.Length)
   108  		return 0
   109  	}
   110  }
   111  
   112  func SortBlastz(seeds []*SeedDev, read fastq.Fastq, scoreMatrix [][]int64) {
   113  	sort.Slice(seeds, func(i, j int) bool { return CompareBlastScore(seeds[i], seeds[j], read, scoreMatrix) == 1 })
   114  }
   115  
   116  func BlastSeed(seed *SeedDev, read fastq.Fastq, scoreMatrix [][]int64) int64 {
   117  	if seed.NextPart == nil {
   118  		return scoreSeed(seed, read, scoreMatrix)
   119  	} else {
   120  		return scoreSeed(seed, read, scoreMatrix) + scoreSeed(seed.NextPart, read, scoreMatrix)
   121  	}
   122  }