t/CXGN/Cluster/cluster.t

   1 #!/usr/bin/perl
   2 use strict;
   3 use warnings;
   4
   5 use FindBin;
   6 use File::Temp qw/tempfile/;
   7 use Data::Dumper;
   8
   9 use Test::More;
  10
  11 use Bio::Index::Fasta;
  12
  13 use CXGN::Tools::List qw/str_in/;
  14 use IPC::Cmd qw/can_run/;
  15
  16 BEGIN {
  17   use_ok( 'CXGN::Cluster' )
  18     or plan skip_all => 'could not include the module being tested';
  19
  20   unless( can_run('phrap.longreads') ){
  21     plan skip_all => 'phrap.longreads not available';
  22   }
  23 }
  24
  25 plan tests => 19;
  26
  27 #now start testing with them
  28 my $set = CXGN::Cluster::ClusterSet->new;
  29
  30 my @seqnames = sort qw/monkey bonobo homo orangutan chimpanzee gorilla marmoset lemur/;
  31
  32 $set->add_match(@seqnames[0,1]);
  33 is_deeply( cs_contents($set),
  34            [[ @seqnames[0,1] ]],
  35            'one cluster with 2 members',
  36          );
  37
  38
  39 $set->add_match(@seqnames[0,2]);
  40 is_deeply( cs_contents($set),
  41            [[ @seqnames[0,1,2] ]],
  42            'then three members',
  43          );
  44
  45 $set->add_match(@seqnames[3,3]);
  46 is_deeply( cs_contents($set),
  47            [[ @seqnames[0,1,2] ],
  48             [ $seqnames[3] ],
  49            ],
  50            'self-match adds an unrelated cluster',
  51          );
  52
  53 $set->add_match(@seqnames[0,0]);
  54 is_deeply( cs_contents($set),
  55            [[ @seqnames[0,1,2] ],
  56             [ $seqnames[3] ],
  57            ],
  58            'self-match of already present sequence does not change anything',
  59          );
  60
  61 $set->add_match(@seqnames[1,1]);
  62 is_deeply( cs_contents($set),
  63            [[ @seqnames[0,1,2] ],
  64             [ $seqnames[3] ],
  65            ],
  66            'self-match of already present sequence does not change anything',
  67          );
  68
  69
  70 $set->add_match(@seqnames[3,7]);
  71 $set->add_match(@seqnames[3,6]);
  72 $set->add_match(@seqnames[3,5]);
  73 is_deeply( cs_contents($set),
  74            [
  75             [ @seqnames[3,5,6,7] ],
  76             [ @seqnames[0,1,2] ],
  77            ],
  78            'and add more to the second cluster',
  79          );
  80
  81
  82
  83 #find our file of test sequences
  84 my $test_seqs_filename = File::Spec->catfile($FindBin::Bin,'data','test.seq');
  85 die "can't open $test_seqs_filename" unless -r $test_seqs_filename;
  86
  87 #make a Bio::Index of them so we can retrieve individual ones
  88 my $test_seqs_index = Bio::Index::Fasta->new( -filename => do { my (undef,$tf) = tempfile(UNLINK => 1); $tf },
  89                                               -write_flag => 1
  90                                             );
  91 $test_seqs_index->make_index($test_seqs_filename);
  92
  93 #make a new set
  94 $set = CXGN::Cluster::ClusterSet->new;
  95
  96 #now test the cluster calculation using a few things we know should cluster in the test set
  97 my @known_clusters =
  98   (
  99    [sort qw/C04HBa0114G11.1 C04HBa0050I18.1 C04HBa0036C23.1 C04HBa0008H22.1/],
 100    [sort qw/C04HBa0024G05.1 C04HBa0020F17.1 /],
 101    [sort qw/C10HBa0020A12.1 C10HBa0248A13.3 /],
 102   );
 103
 104 #make a cluster with just one member
 105 $set->add_match($known_clusters[0][0],$known_clusters[0][0]);
 106 my @c = $set->get_clusters;
 107 is_deeply(cs_contents($set),
 108           [[$known_clusters[0][0]]],
 109           'single-member cluster loaded OK',
 110          );
 111
 112 is_deeply( [ $c[0]->get_contig_coords($test_seqs_index) ],
 113            [[['C04HBa0008H22.1',1,81572,1]]],
 114            'singleton contig coordinates appear correctly',
 115          );
 116
 117 #put in the first cluster
 118 $set->add_match($known_clusters[0][0],$_) foreach @{$known_clusters[0]};;
 119 is_deeply(cs_contents($set),
 120           [$known_clusters[0]],
 121           'first known cluster members loaded ok'
 122          );
 123
 124 #now get its coordinates, it should come out as just one contig
 125 @c = $set->get_clusters;
 126 #print Dumper $c[0]->get_contig_coords($test_seqs_index);
 127 is_deeply( [ $c[0]->get_contig_coords($test_seqs_index) ],
 128            [[
 129             [
 130              'C04HBa0114G11.1',
 131              1,
 132              '28696',
 133              1
 134             ],
 135             [
 136              'C04HBa0050I18.1',
 137              26697,
 138              '140337',
 139              1
 140             ],
 141             [
 142              'C04HBa0036C23.1',
 143              138338,
 144              '252489',
 145              1
 146             ],
 147             [
 148              'C04HBa0008H22.1',
 149              250490,
 150              '332061',
 151              1
 152             ]
 153            ]],
 154            'first test contig assembled OK'
 155          );
 156
 157
 158 my $bs = [ $c[0]->get_consensus_base_segments($test_seqs_index, min_segment_size => 100 ) ];
 159 is_deeply( $bs,
 160            [
 161                [
 162                    [
 163                        '1',
 164                        '26696',
 165                        'C04HBa0114G11.1',
 166                        1,
 167                        26696,
 168                        0
 169                       ],
 170                    [
 171                        '26697',
 172                        '138337',
 173                        'C04HBa0050I18.1',
 174                        1,
 175                        111641,
 176                        0
 177                       ],
 178                    [
 179                        '138338',
 180                        '250489',
 181                        'C04HBa0036C23.1',
 182                        1,
 183                        112152,
 184                        0
 185                       ],
 186                    [
 187                        '250490',
 188                        '332061',
 189                        'C04HBa0008H22.1',
 190                        1,
 191                        81572,
 192                        0
 193                       ]
 194                   ]
 195               ],
 196            'get_consensus_base_segments',
 197           )
 198     or diag Dumper $bs;
 199
 200 #put in the second cluster
 201 $set->add_match($known_clusters[1][0],$_) foreach @{$known_clusters[1]};
 202 is_deeply(cs_contents($set),
 203           [@known_clusters[0,1]],
 204           'second known cluster members loaded ok'
 205          );
 206
 207 #assemble the second cluster
 208 @c = $set->get_clusters;
 209 is_deeply( [ $c[0]->get_contig_coords($test_seqs_index) ],
 210            [[
 211              [
 212               'C04HBa0024G05.1',
 213               1,
 214               '85242',
 215               1
 216              ],
 217              [
 218               'C04HBa0020F17.1',
 219               83243,
 220               '211703',
 221               1
 222              ]
 223             ]],
 224            'second test contig assembled OK'
 225          );
 226
 227 $bs = [ $c[0]->get_consensus_base_segments($test_seqs_index) ];
 228 is_deeply( $bs,
 229            [
 230                [
 231                    [
 232                        '1',
 233                        '83242',
 234                        'C04HBa0024G05.1',
 235                        1,
 236                        83242,
 237                        0
 238                       ],
 239                    [
 240                        '83243',
 241                        '211703',
 242                        'C04HBa0020F17.1',
 243                        1,
 244                        128461,
 245                        0
 246                       ]
 247                   ]
 248            ],
 249           )
 250     or diag Dumper $bs;
 251
 252
 253 #make an artificial linkage between the first and second cluster, this
 254 #should assemble into two real contigs
 255 $set->add_match($known_clusters[0][0],$known_clusters[1][1]);
 256 is_deeply(cs_contents($set),
 257           [[sort @{$known_clusters[0]},@{$known_clusters[1]}]],
 258           'artificial cluster linkage results in joining of clusters'
 259          );
 260
 261 #now check that they assemble into two actual contigs
 262 @c = $set->get_clusters;
 263 #bprint Dumper $c[0]->get_contig_coords($test_seqs_index);
 264 is_deeply( [ $c[0]->get_contig_coords($test_seqs_index) ],
 265            [
 266             [
 267              [
 268               'C04HBa0024G05.1',
 269               1,
 270               '85242',
 271               1
 272              ],
 273              [
 274               'C04HBa0020F17.1',
 275               83243,
 276               '211703',
 277               1
 278              ]
 279             ],
 280             [
 281              [
 282               'C04HBa0114G11.1',
 283               1,
 284               '28696',
 285               1
 286              ],
 287              [
 288               'C04HBa0050I18.1',
 289               26697,
 290               '140337',
 291               1
 292              ],
 293              [
 294               'C04HBa0036C23.1',
 295               138338,
 296               '252489',
 297               1
 298              ],
 299              [
 300               'C04HBa0008H22.1',
 301               250490,
 302               '332061',
 303               1
 304              ]
 305             ]
 306            ],
 307            'erroneous precluster assembles into two actual contigs'
 308          );
 309
 310
 311 #make a new set
 312 $set = CXGN::Cluster::ClusterSet->new;
 313
 314 # now add fake matches between and among the known clusters to put
 315 # them all in one precluster
 316 $set->add_match( $known_clusters[0][0], $known_clusters[2][0] );
 317 $set->add_match( $known_clusters[2][0], $_ ) for @{$known_clusters[2]};
 318
 319 @c = $set->get_clusters;
 320 is( scalar(@c), 1, 'got one cluster' );
 321 # and check the base segments that are made
 322 $bs = [ $c[0]->get_consensus_base_segments($test_seqs_index, min_segment_size => 100 ) ];
 323 is_deeply( $bs,
 324            [
 325                [
 326                    [
 327                        1,
 328                        134537,
 329                        'C10HBa0020A12.1',
 330                        1,
 331                        134537,
 332                        0
 333                       ],
 334                    [
 335                        134538,
 336                        134541,
 337                        'C10HBa0248A13.3',
 338                        134538,
 339                        134541,
 340                        0
 341                       ],
 342                    [
 343                        134542,
 344                        135038,
 345                        'C10HBa0020A12.1',
 346                        134542,
 347                        135038,
 348                        0
 349                       ]
 350                   ],
 351                [
 352                    [
 353                        1,
 354                        81572,
 355                        'C04HBa0008H22.1',
 356                        1,
 357                        81572,
 358                        0
 359                       ]
 360                   ]
 361               ],
 362            'segment merging worked'
 363           )
 364     or diag Dumper $bs;
 365
 366 #make a sorted list of sorted arrayrefs representing the clusters in the set
 367 #clusters sorted in descending size order,
 368 #members sorted in alphabetical order
 369 sub cs_contents {
 370   my $set = shift;
 371   my @c = $set->get_clusters;
 372   return [ sort {@$b <=> @$a}  map [sort $_->get_members], @c ];
 373 }