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

     1  package genomeGraph
     2  
     3  import (
     4  	"fmt"
     5  	"log"
     6  	"testing"
     7  
     8  	"github.com/vertgenlab/gonomics/cigar"
     9  	"github.com/vertgenlab/gonomics/dna"
    10  	"github.com/vertgenlab/gonomics/dna/dnaTwoBit"
    11  	"github.com/vertgenlab/gonomics/giraf"
    12  )
    13  
    14  // TestGraph Structure
    15  //             n2          e0 = 1
    16  //         e1/    \e2      e1 = 0.05
    17  //      e0  /  e3  \       e2 = 1
    18  //  n0 --- n1 ----- n4     e3 = 0.8
    19  //          \      /       e4 = 0.15
    20  //         e4\    /e5      e5 = 1
    21  //             n3
    22  //
    23  //               A
    24  //             /    \
    25  //            /      \
    26  //  ATG --- CG ----- TAA
    27  //            \      /
    28  //             \    /
    29  //               T
    30  
    31  // Sam Header:
    32  //@HD     VN:1.6  SO:coordinate
    33  //@SQ     SN:n0 LN:3
    34  //@SQ     SN:n1 LN:2
    35  //@SQ     SN:n2 LN:1
    36  //@SQ     SN:n3 LN:1
    37  //@SQ     SN:n4 LN:3
    38  
    39  // Test Functions.
    40  func MakeTestGraph() *GenomeGraph {
    41  	graph := EmptyGraph()
    42  
    43  	var n0, n1, n2, n3, n4 Node
    44  	var e0, e1, e2, e3, e4, e5 Edge
    45  
    46  	// Make Nodes
    47  	n0 = Node{
    48  		Id:  0,
    49  		Seq: dna.StringToBases("ATG")}
    50  
    51  	n1 = Node{
    52  		Id:  1,
    53  		Seq: dna.StringToBases("CG")}
    54  
    55  	n2 = Node{
    56  		Id:  2,
    57  		Seq: dna.StringToBases("A")}
    58  
    59  	n3 = Node{
    60  		Id:  3,
    61  		Seq: dna.StringToBases("T")}
    62  
    63  	n4 = Node{
    64  		Id:  4,
    65  		Seq: dna.StringToBases("TAA")}
    66  
    67  	// Make Edges
    68  	e0 = Edge{
    69  		Dest: &n1,
    70  		Prob: 1}
    71  
    72  	e1 = Edge{
    73  		Dest: &n2,
    74  		Prob: 0.05}
    75  
    76  	e2 = Edge{
    77  		Dest: &n4,
    78  		Prob: 1}
    79  
    80  	e3 = Edge{
    81  		Dest: &n4,
    82  		Prob: 0.8}
    83  
    84  	e4 = Edge{
    85  		Dest: &n3,
    86  		Prob: 0.15}
    87  
    88  	e5 = Edge{
    89  		Dest: &n4,
    90  		Prob: 1}
    91  
    92  	// Define Paths
    93  	n0.Next = append(n0.Next, e0)
    94  	n1.Next = append(n1.Next, e1, e3, e4)
    95  	n1.Prev = append(n1.Prev, e0)
    96  	n2.Next = append(n2.Next, e2)
    97  	n2.Prev = append(n2.Prev, e1)
    98  	n3.Next = append(n3.Next, e5)
    99  	n3.Prev = append(n3.Prev, e4)
   100  	n4.Prev = append(n4.Prev, e2, e3, e5)
   101  
   102  	graph.Nodes = append(graph.Nodes, n0, n1, n2, n3, n4)
   103  
   104  	for i := 0; i < len(graph.Nodes); i++ {
   105  		graph.Nodes[i].SeqTwoBit = dnaTwoBit.NewTwoBit(graph.Nodes[i].Seq)
   106  	}
   107  
   108  	return graph
   109  }
   110  
   111  // check struct generated with parameters reads := RandGiraf(MakeTestGraph(), 1, 4, seed).
   112  var check = giraf.Giraf{
   113  	QName:     "0_3_2_2_-",
   114  	QStart:    0,
   115  	QEnd:      4,
   116  	PosStrand: false,                                                    // rev strand, must reverse complement
   117  	Path:      giraf.Path{TStart: 2, Nodes: []uint32{0, 1, 2}, TEnd: 1}, // Nodes 0->1->2, start base 3, end base 1
   118  	Cigar:     []cigar.ByteCigar{{RunLen: 4, Op: 'M'}},
   119  	AlnScore:  16607,
   120  	MapQ:      30,
   121  	Seq:       []dna.Base{3, 1, 2, 1}, // TCGC
   122  	Qual:      []uint8{16, 38, 38, 36},
   123  	Notes:     nil}
   124  
   125  func TestRandGiraf(t *testing.T) {
   126  	var seed int64 = 777
   127  
   128  	// Uncomment following line for random reads
   129  	//seed = time.Now().UnixNano()
   130  
   131  	reads := RandGiraf(MakeTestGraph(), 10, 4, seed)
   132  
   133  	if reads[0].QName != check.QName {
   134  		log.Fatalln("Reads do not match")
   135  	}
   136  	if reads[0].PosStrand != check.PosStrand {
   137  		log.Fatalln("Reads do not match")
   138  	}
   139  	if reads[0].Path.TStart != check.Path.TStart ||
   140  		reads[0].Path.TEnd != check.Path.TEnd ||
   141  		reads[0].Path.Nodes[0] != check.Path.Nodes[0] ||
   142  		reads[0].Path.Nodes[len(reads[0].Path.Nodes)-1] != check.Path.Nodes[len(reads[0].Path.Nodes)-1] {
   143  		log.Fatalln("Reads do not match")
   144  	}
   145  	if dna.CompareSeqsIgnoreCase(reads[0].Seq, check.Seq) == 1 {
   146  		log.Fatalln("Reads do not match")
   147  	}
   148  
   149  	fmt.Println("Sequences Match")
   150  
   151  	RandSomaticMutations(MakeTestGraph(), reads, 1, 0.75, seed)
   152  
   153  	if reads[0].Seq[1] != 1 {
   154  		log.Fatalln("Problem with simulating somatic mutations")
   155  	}
   156  
   157  	fmt.Println("Somatic Mutations Generated Correctly")
   158  
   159  	/*
   160  		for i := 0; i < len(reads); i++ {
   161  			fmt.Println(reads[i])
   162  			fmt.Println(reads[i].Path, reads[i].Aln[0])
   163  		}
   164  	*/
   165  }