Browse code

added code to allow genotype calling for Illumina chips without the need for region specific package

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@126739 bc3139a8-67e5-0310-9ffc-ced21a209358

unknown authored on 14/02/2017 03:47:52
Showing1 changed files
... ...
@@ -422,12 +422,12 @@ SEXP krlmmComputeM(SEXP A, SEXP B){
422 422
 
423 423
     INTERNAL VARIABLES
424 424
     -------- ---------
425
-    rowsAB, colsAB: dimensions of A and B, which is number of SNP and number of sample
425
+    num_SNP, num_sample: dimensions of A and B, which is number of SNP and number of sample
426 426
   */
427 427
 
428
-  long rowsAB, colsAB;
429
-  rowsAB = INTEGER(getAttrib(A, R_DimSymbol))[0];
430
-  colsAB = INTEGER(getAttrib(A, R_DimSymbol))[1];
428
+  long num_SNP, num_sample;
429
+  num_SNP = INTEGER(getAttrib(A, R_DimSymbol))[0];
430
+  num_sample = INTEGER(getAttrib(A, R_DimSymbol))[1];
431 431
 
432 432
   A = coerceVector(A, REALSXP);
433 433
   B = coerceVector(B, REALSXP);
... ...
@@ -435,16 +435,16 @@ SEXP krlmmComputeM(SEXP A, SEXP B){
435 435
   long i, j;
436 436
 
437 437
   SEXP Rval;
438
-  PROTECT(Rval = allocMatrix(REALSXP, rowsAB, colsAB));
438
+  PROTECT(Rval = allocMatrix(REALSXP, num_SNP, num_sample));
439 439
 
440 440
   double *ptr2M;
441 441
   long elepos;
442 442
   ptr2M = NUMERIC_POINTER(Rval);
443 443
 
444
-  for (i = 1; i <= rowsAB; i++){
445
-    for (j = 1; j <= colsAB; j++){
444
+  for (i = 1; i <= num_SNP; i++){
445
+    for (j = 1; j <= num_sample; j++){
446 446
       // elepos is an index for A, B and M
447
-      elepos = CMatrixElementPosition(i, j, rowsAB);
447
+      elepos = CMatrixElementPosition(i, j, num_SNP);
448 448
       ptr2M[elepos] = (log2(REAL(A)[elepos])-log2(REAL(B)[elepos]));
449 449
     }
450 450
   }
... ...
@@ -459,34 +459,34 @@ SEXP krlmmComputeS(SEXP A, SEXP B){
459 459
   /*
460 460
     ARGUMENTS
461 461
     ---------
462
-    A: intensity matrix for allele A11
462
+    A: intensity matrix for allele A
463 463
     B: intensity matrix for allele B
464 464
     S: average log-intensity for a particular SNP (outgoing)
465 465
 
466 466
     INTERNAL VARIABLES
467 467
     -------- ---------
468
-    rowsAB, colsAB: dimensions of A and B, which is number of SNP and number of sample
468
+    num_SNP, num_sample: dimensions of A and B, which is number of SNP and number of sample
469 469
   */
470 470
 
471
-  long rowsAB, colsAB;
472
-  rowsAB = INTEGER(getAttrib(A, R_DimSymbol))[0];
473
-  colsAB = INTEGER(getAttrib(A, R_DimSymbol))[1];
471
+  long num_SNP, num_sample;
472
+  num_SNP = INTEGER(getAttrib(A, R_DimSymbol))[0];
473
+  num_sample = INTEGER(getAttrib(A, R_DimSymbol))[1];
474 474
 
475 475
   A = coerceVector(A, REALSXP);
476 476
   B = coerceVector(B, REALSXP); 
477 477
   long i, j;
478 478
 
479 479
   SEXP Rval;
480
-  PROTECT(Rval = allocMatrix(REALSXP, rowsAB, colsAB));
480
+  PROTECT(Rval = allocMatrix(REALSXP, num_SNP, num_sample));
481 481
 
482 482
   double *ptr2S;
483 483
   long elepos;
484 484
   ptr2S = NUMERIC_POINTER(Rval);
485 485
 
486
-  for (i = 1; i <= rowsAB; i++){
487
-    for (j = 1; j <= colsAB; j++){
486
+  for (i = 1; i <= num_SNP; i++){
487
+    for (j = 1; j <= num_sample; j++){
488 488
       // elepos is an index for A, B and S
489
-      elepos = CMatrixElementPosition(i, j, rowsAB);
489
+      elepos = CMatrixElementPosition(i, j, num_SNP);
490 490
       ptr2S[elepos] = (log2(REAL(A)[elepos])+log2(REAL(B)[elepos]))/2;
491 491
     }
492 492
   }
... ...
@@ -495,15 +495,16 @@ SEXP krlmmComputeS(SEXP A, SEXP B){
495 495
   return Rval;
496 496
 }
497 497
 
498
-void calculate_multiple_cluster_scores(int row, double *intensity, double mean_intensity, int *clustering, int num_SNP, int num_sample, int *ptr, double **dist, int *clustercount){
499
-    int p, q, o;
498
+void calculate_multiple_cluster_scores(long row, double *intensity, double mean_intensity, int *clustering, long num_SNP, long num_sample, int *ptr, double **dist, long *clustercount){
499
+    long p, q;
500
+    int o;
500 501
     long vectorPos;
501 502
     double a_dist;
502 503
     for(p=1; p <= num_sample; p++){
503 504
         dist[p][p] = 0;
504 505
     }
505
-    for(p=1; p<num_sample; p++){
506
-        for(q=p+1; q<=num_sample; q++){
506
+    for(p=1; p < num_sample; p++){
507
+        for(q = p+1; q <= num_sample; q++){
507 508
    	    a_dist = fabs(intensity[CMatrixElementPosition(row, p, num_SNP)] - intensity[CMatrixElementPosition(row, q, num_SNP)]);
508 509
             dist[p][q] = a_dist;
509 510
             dist[q][p] = a_dist;
... ...
@@ -521,7 +522,7 @@ void calculate_multiple_cluster_scores(int row, double *intensity, double mean_i
521 522
         sum[2] = 0;
522 523
         sum[3] = 0;
523 524
         for (q=1; q<=num_sample; q++){
524
-	  sum[clustering[CMatrixElementPosition(row, q, num_SNP)]] += dist[p][q];
525
+            sum[clustering[CMatrixElementPosition(row, q, num_SNP)]] += dist[p][q];
525 526
         }
526 527
        
527 528
         within_cluster = sum[clustering[vectorPos]] / (clustercount[clustering[vectorPos]]- 1);
... ...
@@ -546,9 +547,9 @@ void calculate_multiple_cluster_scores(int row, double *intensity, double mean_i
546 547
     }   
547 548
 }
548 549
 
549
-void calculate_unique_cluster_scores(int row, double *intensity, double mean_intensity, double intensity_range, int num_SNP, int num_sample, int *ptr)				   
550
+void calculate_unique_cluster_scores(long row, double *intensity, double mean_intensity, double intensity_range, long num_SNP, long num_sample, int *ptr)				   
550 551
 {
551
-    int p;
552
+    long p;
552 553
     long vectorPos;
553 554
 
554 555
     for(p = 1; p <= num_sample; p++){
... ...
@@ -558,24 +559,24 @@ void calculate_unique_cluster_scores(int row, double *intensity, double mean_int
558 559
 }
559 560
 
560 561
 
561
-double calculate_SNP_mean(int row, double *intensity, int num_SNP, int num_sample)
562
+double calculate_SNP_mean(long row, double *intensity, long num_SNP, long num_sample)
562 563
 {
563 564
   double sum_intensity;
564 565
   sum_intensity = 0;
565
-  int p;
566
+  long p;
566 567
   for(p = 1; p <= num_sample; p++){
567 568
     sum_intensity = sum_intensity + intensity[CMatrixElementPosition(row, p, num_SNP)];
568 569
   }
569 570
   return(sum_intensity / num_sample);
570 571
 }
571 572
 
572
-double calculate_SNP_range(int row, double *intensity, int num_SNP, int num_sample)
573
+double calculate_SNP_range(long row, double *intensity, long num_SNP, long num_sample)
573 574
 {
574 575
   double min_intensity;
575 576
   double max_intensity;
576 577
   max_intensity = intensity[CMatrixElementPosition(row, 1, num_SNP)];
577 578
   min_intensity = intensity[CMatrixElementPosition(row, 1, num_SNP)];
578
-  int p;  
579
+  long p;  
579 580
   for(p = 2; p <= num_sample; p++){
580 581
     if (intensity[CMatrixElementPosition(row, p, num_SNP)] < min_intensity){
581 582
       min_intensity = intensity[CMatrixElementPosition(row, p, num_SNP)];  
... ...
@@ -590,7 +591,7 @@ double calculate_SNP_range(int row, double *intensity, int num_SNP, int num_samp
590 591
 
591 592
 SEXP krlmmConfidenceScore(SEXP M, SEXP clustering)
592 593
 {
593
-    int num_SNP, num_sample;
594
+    long num_SNP, num_sample;
594 595
     num_SNP = INTEGER(getAttrib(M, R_DimSymbol))[0];
595 596
     num_sample = INTEGER(getAttrib(M, R_DimSymbol))[1];
596 597
 
... ...
@@ -598,7 +599,7 @@ SEXP krlmmConfidenceScore(SEXP M, SEXP clustering)
598 599
     ptr2cluster = INTEGER_POINTER(AS_INTEGER(clustering));
599 600
     M = coerceVector(M, REALSXP);
600 601
   
601
-    int i, j;
602
+    long i, j;
602 603
     int cluster;
603 604
     double mean_intensity;
604 605
     double intensity_range;
... ...
@@ -617,34 +618,31 @@ SEXP krlmmConfidenceScore(SEXP M, SEXP clustering)
617 618
       dist[i] = (double *)malloc((num_sample + 1) * sizeof(double)); 
618 619
     }
619 620
 
620
-    int cluster_count[4]; // extra element to cope with zero-base notation
621
-   // int clid[4];
621
+    long cluster_count[4]; // extra element to cope with zero-base notation
622 622
 
623 623
     for (i = 1; i <= num_SNP; i++) {
624 624
         cluster_count[1] = 0;
625 625
         cluster_count[2] = 0;
626 626
         cluster_count[3] = 0;
627 627
         for (j=1; j<= num_sample; j++){
628
-	    cluster = ptr2cluster[CMatrixElementPosition(i, j, num_SNP)];
628
+            cluster = ptr2cluster[CMatrixElementPosition(i, j, num_SNP)];
629 629
             cluster_count[cluster]++;
630 630
         }
631
-	mean_intensity = calculate_SNP_mean(i, REAL(M), num_SNP, num_sample);
632
-	intensity_range = calculate_SNP_range(i, REAL(M), num_SNP, num_sample);
631
+        mean_intensity = calculate_SNP_mean(i, REAL(M), num_SNP, num_sample);
632
+        intensity_range = calculate_SNP_range(i, REAL(M), num_SNP, num_sample);
633 633
 
634 634
         k = 0;
635 635
         for (j=1; j<=3; j++){
636 636
             if (cluster_count[j] > 0){
637 637
                 k++;
638
-            //    clid[k] = j;
639 638
             }
640 639
         }
641 640
 
642
-	if (k==1) {
643
-	  calculate_unique_cluster_scores(i, REAL(M), mean_intensity, intensity_range, num_SNP, num_sample, ptr2score);
644
-	} else {
645
-	  calculate_multiple_cluster_scores(i, REAL(M), mean_intensity, ptr2cluster, num_SNP, num_sample, ptr2score, dist, cluster_count);
646
-	}       
647
-      
641
+        if (k==1) {
642
+            calculate_unique_cluster_scores(i, REAL(M), mean_intensity, intensity_range, num_SNP, num_sample, ptr2score);
643
+        } else {
644
+	        calculate_multiple_cluster_scores(i, REAL(M), mean_intensity, ptr2cluster, num_SNP, num_sample, ptr2score, dist, cluster_count);			
645
+        }       
648 646
     } 
649 647
 
650 648
     UNPROTECT(1);
Browse code

Added nopackage option for krlmm

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@117087 bc3139a8-67e5-0310-9ffc-ced21a209358

unknown authored on 03/05/2016 23:39:54
Showing1 changed files
... ...
@@ -425,28 +425,27 @@ SEXP krlmmComputeM(SEXP A, SEXP B){
425 425
     rowsAB, colsAB: dimensions of A and B, which is number of SNP and number of sample
426 426
   */
427 427
 
428
-  int rowsAB, colsAB;
428
+  long rowsAB, colsAB;
429 429
   rowsAB = INTEGER(getAttrib(A, R_DimSymbol))[0];
430 430
   colsAB = INTEGER(getAttrib(A, R_DimSymbol))[1];
431 431
 
432
-  int i, j;
433
-
434
-  int *ptr2A, *ptr2B;
435
-  ptr2A = INTEGER_POINTER(AS_INTEGER(A));
436
-  ptr2B = INTEGER_POINTER(AS_INTEGER(B));
432
+  A = coerceVector(A, REALSXP);
433
+  B = coerceVector(B, REALSXP);
434
+  
435
+  long i, j;
437 436
 
438 437
   SEXP Rval;
439 438
   PROTECT(Rval = allocMatrix(REALSXP, rowsAB, colsAB));
440 439
 
441 440
   double *ptr2M;
442
-  long ele;
441
+  long elepos;
443 442
   ptr2M = NUMERIC_POINTER(Rval);
444 443
 
445 444
   for (i = 1; i <= rowsAB; i++){
446 445
     for (j = 1; j <= colsAB; j++){
447
-      // elem is an index for A, B and M
448
-      ele = Cmatrix(i, j, rowsAB);
449
-      ptr2M[ele] = (log2(ptr2A[ele])-log2(ptr2B[ele]));
446
+      // elepos is an index for A, B and M
447
+      elepos = CMatrixElementPosition(i, j, rowsAB);
448
+      ptr2M[elepos] = (log2(REAL(A)[elepos])-log2(REAL(B)[elepos]));
450 449
     }
451 450
   }
452 451
   
... ...
@@ -460,7 +459,7 @@ SEXP krlmmComputeS(SEXP A, SEXP B){
460 459
   /*
461 460
     ARGUMENTS
462 461
     ---------
463
-    A: intensity matrix for allele A
462
+    A: intensity matrix for allele A11
464 463
     B: intensity matrix for allele B
465 464
     S: average log-intensity for a particular SNP (outgoing)
466 465
 
... ...
@@ -469,28 +468,26 @@ SEXP krlmmComputeS(SEXP A, SEXP B){
469 468
     rowsAB, colsAB: dimensions of A and B, which is number of SNP and number of sample
470 469
   */
471 470
 
472
-  int rowsAB, colsAB;
471
+  long rowsAB, colsAB;
473 472
   rowsAB = INTEGER(getAttrib(A, R_DimSymbol))[0];
474 473
   colsAB = INTEGER(getAttrib(A, R_DimSymbol))[1];
475 474
 
476
-  int i, j;
477
-
478
-  int *ptr2A, *ptr2B;
479
-  ptr2A = INTEGER_POINTER(AS_INTEGER(A));
480
-  ptr2B = INTEGER_POINTER(AS_INTEGER(B));
475
+  A = coerceVector(A, REALSXP);
476
+  B = coerceVector(B, REALSXP); 
477
+  long i, j;
481 478
 
482 479
   SEXP Rval;
483 480
   PROTECT(Rval = allocMatrix(REALSXP, rowsAB, colsAB));
484 481
 
485 482
   double *ptr2S;
486
-  long ele;
483
+  long elepos;
487 484
   ptr2S = NUMERIC_POINTER(Rval);
488 485
 
489 486
   for (i = 1; i <= rowsAB; i++){
490 487
     for (j = 1; j <= colsAB; j++){
491
-      // elem is an index for A, B and S
492
-      ele = Cmatrix(i, j, rowsAB);
493
-      ptr2S[ele] = (log2(ptr2A[ele])+log2(ptr2B[ele]))/2;
488
+      // elepos is an index for A, B and S
489
+      elepos = CMatrixElementPosition(i, j, rowsAB);
490
+      ptr2S[elepos] = (log2(REAL(A)[elepos])+log2(REAL(B)[elepos]))/2;
494 491
     }
495 492
   }
496 493
   
... ...
@@ -507,7 +504,7 @@ void calculate_multiple_cluster_scores(int row, double *intensity, double mean_i
507 504
     }
508 505
     for(p=1; p<num_sample; p++){
509 506
         for(q=p+1; q<=num_sample; q++){
510
-   	    a_dist = fabs(intensity[Cmatrix(row, p, num_SNP)] - intensity[Cmatrix(row, q, num_SNP)]);
507
+   	    a_dist = fabs(intensity[CMatrixElementPosition(row, p, num_SNP)] - intensity[CMatrixElementPosition(row, q, num_SNP)]);
511 508
             dist[p][q] = a_dist;
512 509
             dist[q][p] = a_dist;
513 510
         }
... ...
@@ -519,12 +516,12 @@ void calculate_multiple_cluster_scores(int row, double *intensity, double mean_i
519 516
     double temp_between_cluster;
520 517
     double max;
521 518
     for (p=1; p <= num_sample; p++){
522
-        vectorPos = Cmatrix(row, p, num_SNP);
519
+        vectorPos = CMatrixElementPosition(row, p, num_SNP);
523 520
         sum[1] = 0;
524 521
         sum[2] = 0;
525 522
         sum[3] = 0;
526 523
         for (q=1; q<=num_sample; q++){
527
-	  sum[clustering[Cmatrix(row, q, num_SNP)]] += dist[p][q];
524
+	  sum[clustering[CMatrixElementPosition(row, q, num_SNP)]] += dist[p][q];
528 525
         }
529 526
        
530 527
         within_cluster = sum[clustering[vectorPos]] / (clustercount[clustering[vectorPos]]- 1);
... ...
@@ -555,7 +552,7 @@ void calculate_unique_cluster_scores(int row, double *intensity, double mean_int
555 552
     long vectorPos;
556 553
 
557 554
     for(p = 1; p <= num_sample; p++){
558
-        vectorPos = Cmatrix(row, p, num_SNP); 
555
+        vectorPos = CMatrixElementPosition(row, p, num_SNP); 
559 556
         ptr[vectorPos] = genotypeConfidence2(1 -  fabs(fabs(intensity[vectorPos] - mean_intensity) / intensity_range));
560 557
     }
561 558
 }
... ...
@@ -567,7 +564,7 @@ double calculate_SNP_mean(int row, double *intensity, int num_SNP, int num_sampl
567 564
   sum_intensity = 0;
568 565
   int p;
569 566
   for(p = 1; p <= num_sample; p++){
570
-    sum_intensity = sum_intensity + intensity[Cmatrix(row, p, num_SNP)];
567
+    sum_intensity = sum_intensity + intensity[CMatrixElementPosition(row, p, num_SNP)];
571 568
   }
572 569
   return(sum_intensity / num_sample);
573 570
 }
... ...
@@ -576,15 +573,15 @@ double calculate_SNP_range(int row, double *intensity, int num_SNP, int num_samp
576 573
 {
577 574
   double min_intensity;
578 575
   double max_intensity;
579
-  max_intensity = intensity[Cmatrix(row, 1, num_SNP)];
580
-  min_intensity = intensity[Cmatrix(row, 1, num_SNP)];
576
+  max_intensity = intensity[CMatrixElementPosition(row, 1, num_SNP)];
577
+  min_intensity = intensity[CMatrixElementPosition(row, 1, num_SNP)];
581 578
   int p;  
582 579
   for(p = 2; p <= num_sample; p++){
583
-    if (intensity[Cmatrix(row, p, num_SNP)] < min_intensity){
584
-      min_intensity = intensity[Cmatrix(row, p, num_SNP)];  
580
+    if (intensity[CMatrixElementPosition(row, p, num_SNP)] < min_intensity){
581
+      min_intensity = intensity[CMatrixElementPosition(row, p, num_SNP)];  
585 582
     }	
586
-    if (intensity[Cmatrix(row, p, num_SNP)] > max_intensity){
587
-      max_intensity = intensity[Cmatrix(row, p, num_SNP)];
583
+    if (intensity[CMatrixElementPosition(row, p, num_SNP)] > max_intensity){
584
+      max_intensity = intensity[CMatrixElementPosition(row, p, num_SNP)];
588 585
     }
589 586
   }
590 587
   return(max_intensity - min_intensity);
... ...
@@ -597,11 +594,10 @@ SEXP krlmmConfidenceScore(SEXP M, SEXP clustering)
597 594
     num_SNP = INTEGER(getAttrib(M, R_DimSymbol))[0];
598 595
     num_sample = INTEGER(getAttrib(M, R_DimSymbol))[1];
599 596
 
600
-    double *ptr2M;
601 597
     int *ptr2cluster;
602
-    ptr2M = NUMERIC_POINTER(AS_NUMERIC(M));
603 598
     ptr2cluster = INTEGER_POINTER(AS_INTEGER(clustering));
604
-
599
+    M = coerceVector(M, REALSXP);
600
+  
605 601
     int i, j;
606 602
     int cluster;
607 603
     double mean_intensity;
... ...
@@ -629,11 +625,11 @@ SEXP krlmmConfidenceScore(SEXP M, SEXP clustering)
629 625
         cluster_count[2] = 0;
630 626
         cluster_count[3] = 0;
631 627
         for (j=1; j<= num_sample; j++){
632
-	    cluster = ptr2cluster[Cmatrix(i, j, num_SNP)];
628
+	    cluster = ptr2cluster[CMatrixElementPosition(i, j, num_SNP)];
633 629
             cluster_count[cluster]++;
634 630
         }
635
-	mean_intensity = calculate_SNP_mean(i, ptr2M, num_SNP, num_sample);
636
-	intensity_range = calculate_SNP_range(i, ptr2M, num_SNP, num_sample);
631
+	mean_intensity = calculate_SNP_mean(i, REAL(M), num_SNP, num_sample);
632
+	intensity_range = calculate_SNP_range(i, REAL(M), num_SNP, num_sample);
637 633
 
638 634
         k = 0;
639 635
         for (j=1; j<=3; j++){
... ...
@@ -644,9 +640,9 @@ SEXP krlmmConfidenceScore(SEXP M, SEXP clustering)
644 640
         }
645 641
 
646 642
 	if (k==1) {
647
-	  calculate_unique_cluster_scores(i, ptr2M, mean_intensity, intensity_range, num_SNP, num_sample, ptr2score);
643
+	  calculate_unique_cluster_scores(i, REAL(M), mean_intensity, intensity_range, num_SNP, num_sample, ptr2score);
648 644
 	} else {
649
-	  calculate_multiple_cluster_scores(i, ptr2M, mean_intensity, ptr2cluster, num_SNP, num_sample, ptr2score, dist, cluster_count);
645
+	  calculate_multiple_cluster_scores(i, REAL(M), mean_intensity, ptr2cluster, num_SNP, num_sample, ptr2score, dist, cluster_count);
650 646
 	}       
651 647
       
652 648
     } 
Browse code

fix conflict in description

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@80947 bc3139a8-67e5-0310-9ffc-ced21a209358

Rob Scharp authored on 30/09/2013 14:07:17
Showing1 changed files
1 1
old mode 100644
2 2
new mode 100755
... ...
@@ -88,7 +88,7 @@ SEXP gtypeCallerPart1(SEXP A, SEXP B, SEXP fIndex, SEXP mIndex,
88 88
   double buffer;
89 89
 
90 90
   // All pointers appear here
91
-  int *ptr2nIndexes, *ptr2A, *ptr2B, *ptr2Ns, *ptr2df, *ptr2mIndex, *ptr2fIndex, *ptr2ncIndexes;
91
+  int *ptr2nIndexes, *ptr2A, *ptr2B, *ptr2Ns, *ptr2df, *ptr2mIndex, *ptr2fIndex; //*ptr2ncIndexes;
92 92
   double *ptr2Smedian, *ptr2knots, *ptr2params, *ptr2Centers, *ptr2Scales, *ptr2probs, *ptr2trim;
93 93
   ptr2nIndexes = INTEGER_POINTER(AS_INTEGER(nIndexes));
94 94
   ptr2A = INTEGER_POINTER(AS_INTEGER(A));
... ...
@@ -97,7 +97,7 @@ SEXP gtypeCallerPart1(SEXP A, SEXP B, SEXP fIndex, SEXP mIndex,
97 97
   ptr2df = INTEGER_POINTER(AS_INTEGER(df));
98 98
   ptr2mIndex = INTEGER_POINTER(AS_INTEGER(mIndex));
99 99
   ptr2fIndex = INTEGER_POINTER(AS_INTEGER(fIndex));
100
-  ptr2ncIndexes = INTEGER_POINTER(AS_INTEGER(ncIndexes));
100
+ // ptr2ncIndexes = INTEGER_POINTER(AS_INTEGER(ncIndexes));
101 101
   ptr2Smedian = NUMERIC_POINTER(AS_NUMERIC(SMEDIAN));
102 102
   ptr2knots = NUMERIC_POINTER(AS_NUMERIC(knots));
103 103
   ptr2params = NUMERIC_POINTER(AS_NUMERIC(mixtureParams));
... ...
@@ -311,8 +311,9 @@ SEXP gtypeCallerPart2(SEXP A, SEXP B, SEXP fIndex, SEXP mIndex,
311 311
   ib3 = GET_LENGTH(noInfo);
312 312
 
313 313
   // All pointers appear here
314
-  int *ptr2nIndexes, *ptr2A, *ptr2B, *ptr2Ns, *ptr2df, *ptr2mIndex, *ptr2fIndex, *ptr2ncIndexes, *ptr2noTraining, *ptr2noInfo;
315
-  double *ptr2Smedian, *ptr2knots, *ptr2params, *ptr2Centers, *ptr2Scales, *ptr2probs, *ptr2trim;
314
+  int *ptr2nIndexes, *ptr2A, *ptr2B, *ptr2Ns, *ptr2df, *ptr2mIndex, *ptr2fIndex; // *ptr2ncIndexes, 
315
+  int *ptr2noTraining, *ptr2noInfo;
316
+  double *ptr2Smedian, *ptr2knots, *ptr2params, *ptr2Centers, *ptr2Scales, *ptr2probs; // *ptr2trim;
316 317
   ptr2nIndexes = INTEGER_POINTER(AS_INTEGER(nIndexes));
317 318
   ptr2A = INTEGER_POINTER(AS_INTEGER(A));
318 319
   ptr2B = INTEGER_POINTER(AS_INTEGER(B));
... ...
@@ -320,7 +321,7 @@ SEXP gtypeCallerPart2(SEXP A, SEXP B, SEXP fIndex, SEXP mIndex,
320 321
   ptr2df = INTEGER_POINTER(AS_INTEGER(df));
321 322
   ptr2mIndex = INTEGER_POINTER(AS_INTEGER(mIndex));
322 323
   ptr2fIndex = INTEGER_POINTER(AS_INTEGER(fIndex));
323
-  ptr2ncIndexes = INTEGER_POINTER(AS_INTEGER(ncIndexes));
324
+ // ptr2ncIndexes = INTEGER_POINTER(AS_INTEGER(ncIndexes));
324 325
   ptr2noTraining = INTEGER_POINTER(AS_INTEGER(noTraining));
325 326
   ptr2noInfo = INTEGER_POINTER(AS_INTEGER(noInfo));
326 327
 
... ...
@@ -330,7 +331,7 @@ SEXP gtypeCallerPart2(SEXP A, SEXP B, SEXP fIndex, SEXP mIndex,
330 331
   ptr2Centers = NUMERIC_POINTER(AS_NUMERIC(theCenters));
331 332
   ptr2Scales = NUMERIC_POINTER(AS_NUMERIC(theScales));
332 333
   ptr2probs = NUMERIC_POINTER(AS_NUMERIC(probs));
333
-  ptr2trim = NUMERIC_POINTER(AS_NUMERIC(trim));
334
+  // ptr2trim = NUMERIC_POINTER(AS_NUMERIC(trim));
334 335
 
335 336
   // End pointers
336 337
 
... ...
@@ -408,3 +409,249 @@ SEXP gtypeCallerPart2(SEXP A, SEXP B, SEXP fIndex, SEXP mIndex,
408 409
   }
409 410
   return(R_NilValue);
410 411
 }
412
+
413
+
414
+SEXP krlmmComputeM(SEXP A, SEXP B){
415
+
416
+  /*
417
+    ARGUMENTS
418
+    ---------
419
+    A: intensity matrix for allele A
420
+    B: intensity matrix for allele B
421
+    M: log-ratio for a particular SNP (outgoing)
422
+
423
+    INTERNAL VARIABLES
424
+    -------- ---------
425
+    rowsAB, colsAB: dimensions of A and B, which is number of SNP and number of sample
426
+  */
427
+
428
+  int rowsAB, colsAB;
429
+  rowsAB = INTEGER(getAttrib(A, R_DimSymbol))[0];
430
+  colsAB = INTEGER(getAttrib(A, R_DimSymbol))[1];
431
+
432
+  int i, j;
433
+
434
+  int *ptr2A, *ptr2B;
435
+  ptr2A = INTEGER_POINTER(AS_INTEGER(A));
436
+  ptr2B = INTEGER_POINTER(AS_INTEGER(B));
437
+
438
+  SEXP Rval;
439
+  PROTECT(Rval = allocMatrix(REALSXP, rowsAB, colsAB));
440
+
441
+  double *ptr2M;
442
+  long ele;
443
+  ptr2M = NUMERIC_POINTER(Rval);
444
+
445
+  for (i = 1; i <= rowsAB; i++){
446
+    for (j = 1; j <= colsAB; j++){
447
+      // elem is an index for A, B and M
448
+      ele = Cmatrix(i, j, rowsAB);
449
+      ptr2M[ele] = (log2(ptr2A[ele])-log2(ptr2B[ele]));
450
+    }
451
+  }
452
+  
453
+  UNPROTECT(1);
454
+  return Rval;
455
+}
456
+
457
+
458
+SEXP krlmmComputeS(SEXP A, SEXP B){
459
+
460
+  /*
461
+    ARGUMENTS
462
+    ---------
463
+    A: intensity matrix for allele A
464
+    B: intensity matrix for allele B
465
+    S: average log-intensity for a particular SNP (outgoing)
466
+
467
+    INTERNAL VARIABLES
468
+    -------- ---------
469
+    rowsAB, colsAB: dimensions of A and B, which is number of SNP and number of sample
470
+  */
471
+
472
+  int rowsAB, colsAB;
473
+  rowsAB = INTEGER(getAttrib(A, R_DimSymbol))[0];
474
+  colsAB = INTEGER(getAttrib(A, R_DimSymbol))[1];
475
+
476
+  int i, j;
477
+
478
+  int *ptr2A, *ptr2B;
479
+  ptr2A = INTEGER_POINTER(AS_INTEGER(A));
480
+  ptr2B = INTEGER_POINTER(AS_INTEGER(B));
481
+
482
+  SEXP Rval;
483
+  PROTECT(Rval = allocMatrix(REALSXP, rowsAB, colsAB));
484
+
485
+  double *ptr2S;
486
+  long ele;
487
+  ptr2S = NUMERIC_POINTER(Rval);
488
+
489
+  for (i = 1; i <= rowsAB; i++){
490
+    for (j = 1; j <= colsAB; j++){
491
+      // elem is an index for A, B and S
492
+      ele = Cmatrix(i, j, rowsAB);
493
+      ptr2S[ele] = (log2(ptr2A[ele])+log2(ptr2B[ele]))/2;
494
+    }
495
+  }
496
+  
497
+  UNPROTECT(1);
498
+  return Rval;
499
+}
500
+
501
+void calculate_multiple_cluster_scores(int row, double *intensity, double mean_intensity, int *clustering, int num_SNP, int num_sample, int *ptr, double **dist, int *clustercount){
502
+    int p, q, o;
503
+    long vectorPos;
504
+    double a_dist;
505
+    for(p=1; p <= num_sample; p++){
506
+        dist[p][p] = 0;
507
+    }
508
+    for(p=1; p<num_sample; p++){
509
+        for(q=p+1; q<=num_sample; q++){
510
+   	    a_dist = fabs(intensity[Cmatrix(row, p, num_SNP)] - intensity[Cmatrix(row, q, num_SNP)]);
511
+            dist[p][q] = a_dist;
512
+            dist[q][p] = a_dist;
513
+        }
514
+    }
515
+
516
+    double sum[4];
517
+    double within_cluster;
518
+    double between_cluster;
519
+    double temp_between_cluster;
520
+    double max;
521
+    for (p=1; p <= num_sample; p++){
522
+        vectorPos = Cmatrix(row, p, num_SNP);
523
+        sum[1] = 0;
524
+        sum[2] = 0;
525
+        sum[3] = 0;
526
+        for (q=1; q<=num_sample; q++){
527
+	  sum[clustering[Cmatrix(row, q, num_SNP)]] += dist[p][q];
528
+        }
529
+       
530
+        within_cluster = sum[clustering[vectorPos]] / (clustercount[clustering[vectorPos]]- 1);
531
+        between_cluster = -1.0;
532
+        for (o=1; o<=3; o++){
533
+            if ((o != clustering[vectorPos]) && (clustercount[o] > 0)){
534
+                temp_between_cluster = sum[o] / clustercount[o];
535
+                if ((between_cluster < 0) || (temp_between_cluster < between_cluster)) {
536
+                    between_cluster = temp_between_cluster;
537
+                }                                           
538
+            }                                               
539
+        }
540
+       
541
+        if (clustercount[clustering[vectorPos]] > 1) {
542
+            if (between_cluster > within_cluster) {
543
+                max = between_cluster;
544
+            } else {
545
+                max = within_cluster;
546
+            }
547
+            ptr[vectorPos] = genotypeConfidence2((between_cluster - within_cluster) / max);                                                   
548
+        }
549
+    }   
550
+}
551
+
552
+void calculate_unique_cluster_scores(int row, double *intensity, double mean_intensity, double intensity_range, int num_SNP, int num_sample, int *ptr)				   
553
+{
554
+    int p;
555
+    long vectorPos;
556
+
557
+    for(p = 1; p <= num_sample; p++){
558
+        vectorPos = Cmatrix(row, p, num_SNP); 
559
+        ptr[vectorPos] = genotypeConfidence2(1 -  fabs(fabs(intensity[vectorPos] - mean_intensity) / intensity_range));
560
+    }
561
+}
562
+
563
+
564
+double calculate_SNP_mean(int row, double *intensity, int num_SNP, int num_sample)
565
+{
566
+  double sum_intensity;
567
+  sum_intensity = 0;
568
+  int p;
569
+  for(p = 1; p <= num_sample; p++){
570
+    sum_intensity = sum_intensity + intensity[Cmatrix(row, p, num_SNP)];
571
+  }
572
+  return(sum_intensity / num_sample);
573
+}
574
+
575
+double calculate_SNP_range(int row, double *intensity, int num_SNP, int num_sample)
576
+{
577
+  double min_intensity;
578
+  double max_intensity;
579
+  max_intensity = intensity[Cmatrix(row, 1, num_SNP)];
580
+  min_intensity = intensity[Cmatrix(row, 1, num_SNP)];
581
+  int p;  
582
+  for(p = 2; p <= num_sample; p++){
583
+    if (intensity[Cmatrix(row, p, num_SNP)] < min_intensity){
584
+      min_intensity = intensity[Cmatrix(row, p, num_SNP)];  
585
+    }	
586
+    if (intensity[Cmatrix(row, p, num_SNP)] > max_intensity){
587
+      max_intensity = intensity[Cmatrix(row, p, num_SNP)];
588
+    }
589
+  }
590
+  return(max_intensity - min_intensity);
591
+}
592
+
593
+
594
+SEXP krlmmConfidenceScore(SEXP M, SEXP clustering)
595
+{
596
+    int num_SNP, num_sample;
597
+    num_SNP = INTEGER(getAttrib(M, R_DimSymbol))[0];
598
+    num_sample = INTEGER(getAttrib(M, R_DimSymbol))[1];
599
+
600
+    double *ptr2M;
601
+    int *ptr2cluster;
602
+    ptr2M = NUMERIC_POINTER(AS_NUMERIC(M));
603
+    ptr2cluster = INTEGER_POINTER(AS_INTEGER(clustering));
604
+
605
+    int i, j;
606
+    int cluster;
607
+    double mean_intensity;
608
+    double intensity_range;
609
+    int k;
610
+    
611
+    SEXP Rval;
612
+    PROTECT(Rval = allocMatrix(INTSXP, num_SNP, num_sample));
613
+
614
+    int *ptr2score;
615
+    ptr2score = INTEGER_POINTER(Rval);
616
+    
617
+    double **dist;
618
+    // allocate memory to dist
619
+    dist = (double **)malloc((num_sample + 1)*sizeof(double *));
620
+    for (i = 1; i <= num_sample; i++){
621
+      dist[i] = (double *)malloc((num_sample + 1) * sizeof(double)); 
622
+    }
623
+
624
+    int cluster_count[4]; // extra element to cope with zero-base notation
625
+   // int clid[4];
626
+
627
+    for (i = 1; i <= num_SNP; i++) {
628
+        cluster_count[1] = 0;
629
+        cluster_count[2] = 0;
630
+        cluster_count[3] = 0;
631
+        for (j=1; j<= num_sample; j++){
632
+	    cluster = ptr2cluster[Cmatrix(i, j, num_SNP)];
633
+            cluster_count[cluster]++;
634
+        }
635
+	mean_intensity = calculate_SNP_mean(i, ptr2M, num_SNP, num_sample);
636
+	intensity_range = calculate_SNP_range(i, ptr2M, num_SNP, num_sample);
637
+
638
+        k = 0;
639
+        for (j=1; j<=3; j++){
640
+            if (cluster_count[j] > 0){
641
+                k++;
642
+            //    clid[k] = j;
643
+            }
644
+        }
645
+
646
+	if (k==1) {
647
+	  calculate_unique_cluster_scores(i, ptr2M, mean_intensity, intensity_range, num_SNP, num_sample, ptr2score);
648
+	} else {
649
+	  calculate_multiple_cluster_scores(i, ptr2M, mean_intensity, ptr2cluster, num_SNP, num_sample, ptr2score, dist, cluster_count);
650
+	}       
651
+      
652
+    } 
653
+
654
+    UNPROTECT(1);
655
+    return Rval;
656
+}
657
+
Browse code

Fixing several warnings and notes

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@40930 bc3139a8-67e5-0310-9ffc-ced21a209358

Benilton Carvalho authored on 31/07/2009 04:51:13
Showing1 changed files
... ...
@@ -84,7 +84,7 @@ SEXP gtypeCallerPart1(SEXP A, SEXP B, SEXP fIndex, SEXP mIndex,
84 84
   //const int lenLists=3;
85 85
 
86 86
   // Buffers
87
-  int intbuffer, ibv1[colsAB], ibv2[colsAB], ib2;
87
+  int intbuffer, ibv1[colsAB], ib2;
88 88
   double buffer;
89 89
 
90 90
   // All pointers appear here
... ...
@@ -303,11 +303,8 @@ SEXP gtypeCallerPart2(SEXP A, SEXP B, SEXP fIndex, SEXP mIndex,
303 303
   colsAB = INTEGER(getAttrib(A, R_DimSymbol))[1];
304 304
   double likelihood[colsAB*3], M[colsAB], S[colsAB], f[colsAB];
305 305
 
306
-  // Constants
307
-  const int lenLists=3;
308
-
309 306
   // Buffers
310
-  int intbuffer, ibv1[colsAB], ibv2[colsAB], ib2, ib3, ibSnpLevel1=0, ibSnpLevel2=0;
307
+  int intbuffer, ib2, ib3, ibSnpLevel1=0, ibSnpLevel2=0;
311 308
   double buffer;
312 309
 
313 310
   ib2 = GET_LENGTH(noTraining);
Browse code

Fixed minor bugs in verbosity

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@36675 bc3139a8-67e5-0310-9ffc-ced21a209358

Benilton Carvalho authored on 15/01/2009 04:04:31
Showing1 changed files
... ...
@@ -131,7 +131,7 @@ SEXP gtypeCallerPart1(SEXP A, SEXP B, SEXP fIndex, SEXP mIndex,
131 131
     int ibv[ib2];
132 132
     if (nSnpsPerClass[h] > 0)
133 133
       for (i=0; i < ptr2nIndexes[h]; i++){
134
-	if (i%100000 == 0) Rprintf("+");
134
+	/* if (i%100000 == 0) Rprintf("+");*/
135 135
 	// Substract 1, as coming from R it is 1-based and C is 0-based.
136 136
 	ii=INTEGER(VECTOR_ELT(Indexes, h))[i] - 1;
137 137
 	for (j=0; j < colsAB; j++){
... ...
@@ -348,7 +348,7 @@ SEXP gtypeCallerPart2(SEXP A, SEXP B, SEXP fIndex, SEXP mIndex,
348 348
   for (h=0; h < nSnpClasses; h++){
349 349
     if (nSnpsPerClass[h] > 0)
350 350
       for (i=0; i < ptr2nIndexes[h]; i++){
351
-	if (i%100000 == 0) Rprintf("+");
351
+	/* if (i%100000 == 0) Rprintf("+"); */
352 352
 	// Substract 1, as coming from R it is 1-based and C is 0-based.
353 353
 	ii=INTEGER(VECTOR_ELT(Indexes, h))[i] - 1;
354 354
 	intbuffer = ii+1;
Browse code

Adding crlmm package

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@36653 bc3139a8-67e5-0310-9ffc-ced21a209358

Benilton Carvalho authored on 14/01/2009 22:00:52
Showing1 changed files
1 1
new file mode 100644
... ...
@@ -0,0 +1,413 @@
1
+#include <math.h>
2
+#include <R.h>
3
+#include <Rdefines.h>
4
+#include <Rmath.h>
5
+#include <Rinternals.h>
6
+
7
+#include "utils.h"
8
+
9
+static double mydt(double x, int df){
10
+  return(pow(1.0+pow(x, 2.0)/ (double)df, -((double)df+1.0)/2.0));
11
+}
12
+
13
+SEXP gtypeCallerPart1(SEXP A, SEXP B, SEXP fIndex, SEXP mIndex,
14
+		      SEXP theCenters, SEXP theScales, SEXP theNs,
15
+		      SEXP Indexes, SEXP cIndexes, SEXP nIndexes,
16
+		      SEXP ncIndexes, SEXP SMEDIAN,
17
+		      SEXP knots, SEXP mixtureParams, SEXP df,
18
+		      SEXP probs, SEXP trim){
19
+  /*
20
+    ARGUMENTS
21
+    ---------
22
+    A: intensity matrix for allele A
23
+    B: intensity matrix for allele B
24
+    fIndex: indexes for females (columns in A/B for females)
25
+    mIndex: indexes for males (columns in A/B for males)
26
+    theCenters: matrix with SNP-specific centers (3 columns: AA/AB/BB)
27
+    theScales: matrix with SNP-specific scales (3 columns: AA/AB/BB)
28
+    theNs: matrix with SNP-specific counts (3 columns: AA/AB/BB)
29
+    Indexes: list with 3 elements (autosomeIndex, XIndex, YIndex) for SNPs
30
+    cIndexes: list with 3 elements (keepIndex, keepIndexFemale, keepIndexMale) for arrays
31
+    SMEDIAN: scalar (median S)
32
+    knots: knots for mixture
33
+    mixtureParams: mixture parameters
34
+    probs: genotype priors (1/3) for *ALL* SNPs. It's a vector of length 3
35
+    trim: drop rate to estimate means
36
+
37
+    ASSUMPTIONS
38
+    -----------
39
+    A and B have the same dimensions
40
+    fIndex and mIndex are in a valid range for A/B
41
+    Number of rows of theCenters, theScales, theNs match the number of rows in A/B
42
+    Length of Indexes and cIndexes is 3
43
+    3 knots
44
+    4xNSAMPLES parameters
45
+    priors has length 3 and is the same for *ALL* SNPs
46
+    trim in (0, .5)
47
+
48
+    INTERNAL VARIABLES
49
+    -------- ---------
50
+    likelihood: matrix (nsample rows x 3 columns)
51
+    rowsAB, colsAB: dimensions of A and B
52
+    lenLists = 3, length of Indexes and cIndexes
53
+    h, i: iteration
54
+    nIndex: length of Indexes[[h]]
55
+    ii: particular value of Indexes
56
+    M: log-ratio for a particular SNP
57
+    S: adjusted average log-intensity for a particular SNP
58
+    f: f for a particular SNP
59
+
60
+    TODO
61
+    ----
62
+    - Get length of a vector from a list within C (adding nIndexes and ncIndexes for the moment)
63
+  */
64
+
65
+  /*
66
+    ==========================================
67
+              VARIABLE DECLARATION
68
+    ==========================================
69
+  */
70
+  // Organizing variables
71
+  int nSnpClasses, k;
72
+  nSnpClasses = GET_LENGTH(Indexes);
73
+  int nSnpsPerClass[nSnpClasses];
74
+  for (k=0; k < nSnpClasses; k++)
75
+    nSnpsPerClass[k] = GET_LENGTH(VECTOR_ELT(Indexes, k));
76
+				  
77
+  // General Variables
78
+  int rowsAB, colsAB, h, i, ii, j, elem, nMales, nFemales;
79
+  rowsAB = INTEGER(getAttrib(A, R_DimSymbol))[0];
80
+  colsAB = INTEGER(getAttrib(A, R_DimSymbol))[1];
81
+  double likelihood[colsAB*3], M[colsAB], S[colsAB], f[colsAB];
82
+
83
+  // Constants
84
+  //const int lenLists=3;
85
+
86
+  // Buffers
87
+  int intbuffer, ibv1[colsAB], ibv2[colsAB], ib2;
88
+  double buffer;
89
+
90
+  // All pointers appear here
91
+  int *ptr2nIndexes, *ptr2A, *ptr2B, *ptr2Ns, *ptr2df, *ptr2mIndex, *ptr2fIndex, *ptr2ncIndexes;
92
+  double *ptr2Smedian, *ptr2knots, *ptr2params, *ptr2Centers, *ptr2Scales, *ptr2probs, *ptr2trim;
93
+  ptr2nIndexes = INTEGER_POINTER(AS_INTEGER(nIndexes));
94
+  ptr2A = INTEGER_POINTER(AS_INTEGER(A));
95
+  ptr2B = INTEGER_POINTER(AS_INTEGER(B));
96
+  ptr2Ns = INTEGER_POINTER(AS_INTEGER(theNs));
97
+  ptr2df = INTEGER_POINTER(AS_INTEGER(df));
98
+  ptr2mIndex = INTEGER_POINTER(AS_INTEGER(mIndex));
99
+  ptr2fIndex = INTEGER_POINTER(AS_INTEGER(fIndex));
100
+  ptr2ncIndexes = INTEGER_POINTER(AS_INTEGER(ncIndexes));
101
+  ptr2Smedian = NUMERIC_POINTER(AS_NUMERIC(SMEDIAN));
102
+  ptr2knots = NUMERIC_POINTER(AS_NUMERIC(knots));
103
+  ptr2params = NUMERIC_POINTER(AS_NUMERIC(mixtureParams));
104
+  ptr2Centers = NUMERIC_POINTER(AS_NUMERIC(theCenters));
105
+  ptr2Scales = NUMERIC_POINTER(AS_NUMERIC(theScales));
106
+  ptr2probs = NUMERIC_POINTER(AS_NUMERIC(probs));
107
+  ptr2trim = NUMERIC_POINTER(AS_NUMERIC(trim));
108
+  // End pointers
109
+
110
+  // These will be returned to R
111
+  double *ptr2e1, *ptr2e2, *ptr2e3;
112
+  SEXP estimates1, estimates2, estimates3, output;
113
+  PROTECT(estimates1 = allocMatrix(REALSXP, rowsAB, 3));
114
+  PROTECT(estimates2 = allocMatrix(REALSXP, rowsAB, 3));
115
+  PROTECT(estimates3 = allocMatrix(REALSXP, rowsAB, 3));
116
+  
117
+  ptr2e1 = NUMERIC_POINTER(estimates1);
118
+  ptr2e2 = NUMERIC_POINTER(estimates2);
119
+  ptr2e3 = NUMERIC_POINTER(estimates3);
120
+  /*
121
+    ==========================================
122
+            END VARIABLE DECLARATION
123
+    ==========================================
124
+  */
125
+  nMales = GET_LENGTH(mIndex);
126
+  nFemales = GET_LENGTH(fIndex);
127
+
128
+  for (h=0; h < nSnpClasses; h++){
129
+    ib2 = GET_LENGTH(VECTOR_ELT(cIndexes, h));
130
+    double dbv[ib2];
131
+    int ibv[ib2];
132
+    if (nSnpsPerClass[h] > 0)
133
+      for (i=0; i < ptr2nIndexes[h]; i++){
134
+	if (i%100000 == 0) Rprintf("+");
135
+	// Substract 1, as coming from R it is 1-based and C is 0-based.
136
+	ii=INTEGER(VECTOR_ELT(Indexes, h))[i] - 1;
137
+	for (j=0; j < colsAB; j++){
138
+	  // j is an index for vectors whose length is number of samples (SAMPLE)
139
+	  // elem is an index for A and B **only** (or objs with SNP rows and SAMPLE columns)
140
+	  //Rprintf("J %d I %d Rows %d\n", j, ii, rowsAB);
141
+
142
+	  elem = rowsAB * j + ii;
143
+	  //Rprintf("\nElemt %d ", elem);
144
+	  
145
+	  M[j] = (log2(ptr2A[elem])-log2(ptr2B[elem]));
146
+	  S[j] = (log2(ptr2A[elem])+log2(ptr2B[elem]))/2 - ptr2Smedian[0];
147
+	  buffer = fmax(fmin(S[j], ptr2knots[2]), ptr2knots[0]);
148
+	  f[j] = ptr2params[j*4+0]+ptr2params[j*4+1]*buffer+ptr2params[j*4+2]*pow(buffer, 2.0)+ptr2params[j*4+3]*pow(fmax(0, buffer-ptr2knots[1]), 2.0);
149
+	  //Rprintf("M %f S %f f %f ", M[j], S[j], f[j]);
150
+
151
+	  // buffer here is sigma
152
+	  // likelihood for AA
153
+	  // All likelihoods already multiplied by prior to save time
154
+	  buffer = ptr2Scales[ii] * sdCorrection(&ptr2Ns[ii]);
155
+	  likelihood[j] = mydt( ((M[j]-f[j])-ptr2Centers[ii])/buffer, ptr2df[0])*ptr2probs[0];
156
+
157
+	  //Rprintf("L1 %2.4f ", likelihood[j]);
158
+
159
+	  // likelihood for AB
160
+	  buffer = ptr2Scales[ii+rowsAB] * sdCorrection(&ptr2Ns[ii+rowsAB]);
161
+	  likelihood[j+colsAB] = mydt( (M[j]-ptr2Centers[ii+rowsAB])/buffer, ptr2df[0])*ptr2probs[1];
162
+	  
163
+	  // intbuffer (here) is the subject ID as in R (ie. 1-based)
164
+	  intbuffer = j+1;
165
+	  if (nMales > 0)
166
+	    if (h > 0)
167
+	      if (intInSet(&intbuffer, ptr2mIndex, &nMales) > 0)
168
+		likelihood[j+colsAB] = 0;
169
+
170
+	  //Rprintf("L2 %2.4f ", likelihood[j+colsAB]);
171
+
172
+
173
+	  // likelihood for BB
174
+	  buffer = ptr2Scales[ii+2*rowsAB] * sdCorrection(&ptr2Ns[ii+2*rowsAB]);
175
+	  likelihood[j+2*colsAB] = mydt( ((M[j]+f[j])-ptr2Centers[ii+2*rowsAB])/buffer, ptr2df[0])*ptr2probs[2];
176
+	  
177
+	  // Females on Y: 1 to avoid NAs. Later made 0 (RI)
178
+	  // To save some time: 1*priors = priors
179
+	  if (nFemales > 0)
180
+	    if (h == 2)
181
+	      if (intInSet(&intbuffer, ptr2fIndex, &nFemales) >0){
182
+		likelihood[j] = ptr2probs[2];
183
+		likelihood[j+colsAB] = ptr2probs[2];
184
+		likelihood[j+2*colsAB] = ptr2probs[2];
185
+	      }
186
+
187
+	  //Rprintf("L3 %2.4f ", likelihood[j+2*colsAB]);
188
+
189
+
190
+	  // Compute simple posterior
191
+	  buffer = likelihood[j]+likelihood[j+colsAB]+likelihood[j+2*colsAB];
192
+	  likelihood[j]/=buffer;
193
+	  likelihood[j+colsAB]/=buffer;
194
+	  likelihood[j+2*colsAB]/=buffer;
195
+
196
+	  if (nFemales > 0)
197
+	    if (h == 2)
198
+	      if (intInSet(&intbuffer, ptr2fIndex, &nFemales) >0){
199
+		likelihood[j] = 0;
200
+		likelihood[j+colsAB] = 0;
201
+		likelihood[j+2*colsAB] = 0;
202
+	      }
203
+
204
+	  ibv1[j] = genotypeCall(&likelihood[j], &likelihood[j+colsAB], &likelihood[j+2*colsAB]);
205
+	}
206
+      
207
+	for (j=0; j < ib2; j++){
208
+	  intbuffer = INTEGER(VECTOR_ELT(cIndexes, h))[j]-1;
209
+	  dbv[j]=M[intbuffer]-f[intbuffer];
210
+	  ibv[j]=ibv1[intbuffer];
211
+	}
212
+	trimmed_mean(dbv, ibv, 1, ptr2trim[0], GET_LENGTH(VECTOR_ELT(cIndexes, h)), rowsAB, ptr2e1, ptr2e2, ptr2e3, ii);
213
+	for (j=0; j < ib2; j++){
214
+	  intbuffer = INTEGER(VECTOR_ELT(cIndexes, h))[j]-1;
215
+	  dbv[j]=M[intbuffer];
216
+	}
217
+	trimmed_mean(dbv, ibv, 2, ptr2trim[0], GET_LENGTH(VECTOR_ELT(cIndexes, h)), rowsAB, ptr2e1, ptr2e2, ptr2e3, ii);
218
+	for (j=0; j < ib2; j++){
219
+	  intbuffer = INTEGER(VECTOR_ELT(cIndexes, h))[j]-1;
220
+	  dbv[j] = M[intbuffer]+f[intbuffer];
221
+	}
222
+	trimmed_mean(dbv, ibv, 3, ptr2trim[0], GET_LENGTH(VECTOR_ELT(cIndexes, h)), rowsAB, ptr2e1, ptr2e2, ptr2e3, ii);
223
+      } /* for Snp */
224
+  } /* for SnpClass */
225
+  PROTECT(output = allocVector(VECSXP,3));
226
+  SET_VECTOR_ELT(output, 0, estimates1);
227
+  SET_VECTOR_ELT(output, 1, estimates2);
228
+  SET_VECTOR_ELT(output, 2, estimates3);
229
+  UNPROTECT(4);
230
+  return(output);
231
+}
232
+
233
+SEXP gtypeCallerPart2(SEXP A, SEXP B, SEXP fIndex, SEXP mIndex,
234
+		      SEXP theCenters, SEXP theScales, SEXP theNs,
235
+		      SEXP Indexes, SEXP cIndexes, SEXP nIndexes,
236
+		      SEXP ncIndexes, SEXP SMEDIAN,
237
+		      SEXP knots, SEXP mixtureParams, SEXP df,
238
+		      SEXP probs, SEXP trim, SEXP noTraining, SEXP noInfo){
239
+  /*
240
+    WARNING!!! REMEMBER TO MODIFY MY TWIN TOO!
241
+
242
+    ARGUMENTS
243
+    ---------
244
+    A: intensity matrix for allele A
245
+    B: intensity matrix for allele B
246
+    fIndex: indexes for females (columns in A/B for females)
247
+    mIndex: indexes for males (columns in A/B for males)
248
+    theCenters: matrix with SNP-specific centers (3 columns: AA/AB/BB)
249
+    theScales: matrix with SNP-specific scales (3 columns: AA/AB/BB)
250
+    theNs: matrix with SNP-specific counts (3 columns: AA/AB/BB)
251
+    Indexes: list with 3 elements (autosomeIndex, XIndex, YIndex) for SNPs
252
+    cIndexes: list with 3 elements (keepIndex, keepIndexFemale, keepIndexMale) for arrays
253
+    SMEDIAN: scalar (median S)
254
+    knots: knots for mixture
255
+    mixtureParams: mixture parameters
256
+    probs: genotype priors (1/3) for *ALL* SNPs. It's a vector of length 3
257
+    trim: drop rate to estimate means
258
+
259
+    ASSUMPTIONS
260
+    -----------
261
+    A and B have the same dimensions
262
+    fIndex and mIndex are in a valid range for A/B
263
+    Number of rows of theCenters, theScales, theNs match the number of rows in A/B
264
+    Length of Indexes and cIndexes is 3
265
+    3 knots
266
+    4xNSAMPLES parameters
267
+    priors has length 3 and is the same for *ALL* SNPs
268
+    trim in (0, .5)
269
+    Indexes and cIndexes have the same length.
270
+
271
+    INTERNAL VARIABLES
272
+    -------- ---------
273
+    likelihood: matrix (nsample rows x 3 columns)
274
+    rowsAB, colsAB: dimensions of A and B
275
+    lenLists = 3, length of Indexes and cIndexes
276
+    h, i: iteration
277
+    nIndex: length of Indexes[[h]]
278
+    ii: particular value of Indexes
279
+    M: log-ratio for a particular SNP
280
+    S: adjusted average log-intensity for a particular SNP
281
+    f: f for a particular SNP
282
+
283
+    TODO
284
+    ----
285
+    - Get length of a vector from a list within C (adding nIndexes and ncIndexes for the moment)
286
+  */
287
+
288
+  /*
289
+    ==========================================
290
+              VARIABLE DECLARATION
291
+    ==========================================
292
+  */
293
+  // Organizing variables
294
+  int nSnpClasses, k;
295
+  nSnpClasses = GET_LENGTH(Indexes);
296
+  int nSnpsPerClass[nSnpClasses];
297
+  for (k=0; k < nSnpClasses; k++)
298
+    nSnpsPerClass[k] = GET_LENGTH(VECTOR_ELT(Indexes, k));
299
+
300
+  // General Variables
301
+  int rowsAB, colsAB, h, i, ii, j, elem, nMales, nFemales;
302
+  rowsAB = INTEGER(getAttrib(A, R_DimSymbol))[0];
303
+  colsAB = INTEGER(getAttrib(A, R_DimSymbol))[1];
304
+  double likelihood[colsAB*3], M[colsAB], S[colsAB], f[colsAB];
305
+
306
+  // Constants
307
+  const int lenLists=3;
308
+
309
+  // Buffers
310
+  int intbuffer, ibv1[colsAB], ibv2[colsAB], ib2, ib3, ibSnpLevel1=0, ibSnpLevel2=0;
311
+  double buffer;
312
+
313
+  ib2 = GET_LENGTH(noTraining);
314
+  ib3 = GET_LENGTH(noInfo);
315
+
316
+  // All pointers appear here
317
+  int *ptr2nIndexes, *ptr2A, *ptr2B, *ptr2Ns, *ptr2df, *ptr2mIndex, *ptr2fIndex, *ptr2ncIndexes, *ptr2noTraining, *ptr2noInfo;
318
+  double *ptr2Smedian, *ptr2knots, *ptr2params, *ptr2Centers, *ptr2Scales, *ptr2probs, *ptr2trim;
319
+  ptr2nIndexes = INTEGER_POINTER(AS_INTEGER(nIndexes));
320
+  ptr2A = INTEGER_POINTER(AS_INTEGER(A));
321
+  ptr2B = INTEGER_POINTER(AS_INTEGER(B));
322
+  ptr2Ns = INTEGER_POINTER(AS_INTEGER(theNs));
323
+  ptr2df = INTEGER_POINTER(AS_INTEGER(df));
324
+  ptr2mIndex = INTEGER_POINTER(AS_INTEGER(mIndex));
325
+  ptr2fIndex = INTEGER_POINTER(AS_INTEGER(fIndex));
326
+  ptr2ncIndexes = INTEGER_POINTER(AS_INTEGER(ncIndexes));
327
+  ptr2noTraining = INTEGER_POINTER(AS_INTEGER(noTraining));
328
+  ptr2noInfo = INTEGER_POINTER(AS_INTEGER(noInfo));
329
+
330
+  ptr2Smedian = NUMERIC_POINTER(AS_NUMERIC(SMEDIAN));
331
+  ptr2knots = NUMERIC_POINTER(AS_NUMERIC(knots));
332
+  ptr2params = NUMERIC_POINTER(AS_NUMERIC(mixtureParams));
333
+  ptr2Centers = NUMERIC_POINTER(AS_NUMERIC(theCenters));
334
+  ptr2Scales = NUMERIC_POINTER(AS_NUMERIC(theScales));
335
+  ptr2probs = NUMERIC_POINTER(AS_NUMERIC(probs));
336
+  ptr2trim = NUMERIC_POINTER(AS_NUMERIC(trim));
337
+
338
+  // End pointers
339
+
340
+  /*
341
+    ==========================================
342
+            END VARIABLE DECLARATION
343
+    ==========================================
344
+  */
345
+  nMales = GET_LENGTH(mIndex);
346
+  nFemales = GET_LENGTH(fIndex);
347
+
348
+  for (h=0; h < nSnpClasses; h++){
349
+    if (nSnpsPerClass[h] > 0)
350
+      for (i=0; i < ptr2nIndexes[h]; i++){
351
+	if (i%100000 == 0) Rprintf("+");
352
+	// Substract 1, as coming from R it is 1-based and C is 0-based.
353
+	ii=INTEGER(VECTOR_ELT(Indexes, h))[i] - 1;
354
+	intbuffer = ii+1;
355
+	if (intInSet(&intbuffer, ptr2noTraining, &ib2) > 0) ibSnpLevel1 = 1;
356
+	if (intInSet(&intbuffer, ptr2noInfo, &ib3) > 0) ibSnpLevel2 = 1;
357
+	for (j=0; j < colsAB; j++){
358
+	  // j is an index for vectors whose length is number of samples (SAMPLE)
359
+	  // elem is an index for A and B **only** (or objs with SNP rows and SAMPLE columns)
360
+	  elem = rowsAB * j + ii;
361
+	  M[j] = (log2(ptr2A[elem])-log2(ptr2B[elem]));
362
+	  S[j] = (log2(ptr2A[elem])+log2(ptr2B[elem]))/2 - ptr2Smedian[0];
363
+	  buffer = fmax(fmin(S[j], ptr2knots[2]), ptr2knots[0]);
364
+	  f[j] = ptr2params[j*4+0]+ptr2params[j*4+1]*buffer+ptr2params[j*4+2]*pow(buffer, 2.0)+ptr2params[j*4+3]*pow(fmax(0, buffer-ptr2knots[1]), 2.0);
365
+	  
366
+	  // buffer here is sigma
367
+	  // likelihood for AA
368
+	  // All likelihoods already multiplied by prior to save time
369
+	  buffer = ptr2Scales[ii] * sdCorrection(&ptr2Ns[ii]);
370
+	  likelihood[j] = mydt( ((M[j]-f[j])-ptr2Centers[ii])/buffer, ptr2df[0])*ptr2probs[0];
371
+
372
+	  // likelihood for AB
373
+	  buffer = ptr2Scales[ii+rowsAB] * sdCorrection(&ptr2Ns[ii+rowsAB]);
374
+	  likelihood[j+colsAB] = mydt( (M[j]-ptr2Centers[ii+rowsAB])/buffer, ptr2df[0])*ptr2probs[1];
375
+	  
376
+	  // intbuffer (here) is the subject ID as in R (ie. 1-based)
377
+	  // h > 0 is chr X or Y
378
+	  intbuffer = j+1;
379
+	  if (nMales > 0)
380
+	    if (h > 0 && intInSet(&intbuffer, ptr2mIndex, &nMales) > 0) likelihood[j+colsAB] = 0;
381
+
382
+	  // likelihood for BB
383
+	  buffer = ptr2Scales[ii+2*rowsAB] * sdCorrection(&ptr2Ns[ii+2*rowsAB]);
384
+	  likelihood[j+2*colsAB] = mydt( ((M[j]+f[j])-ptr2Centers[ii+2*rowsAB])/buffer, ptr2df[0])*ptr2probs[2];
385
+
386
+	  // Females on Y: 1 to avoid NAs. Later made 0 (RI)
387
+	  // To save some time: 1*priors = priors
388
+	  if (nFemales > 0)
389
+	    if (h == 2 && intInSet(&intbuffer, ptr2fIndex, &nFemales) > 0)
390
+	      likelihood[j] = likelihood[j+colsAB] = likelihood[j+2*colsAB] = ptr2probs[2];
391
+	  
392
+	  // Compute simple posterior
393
+	  buffer = likelihood[j]+likelihood[j+colsAB]+likelihood[j+2*colsAB];
394
+	  likelihood[j]/=buffer;
395
+	  likelihood[j+colsAB]/=buffer;
396
+	  likelihood[j+2*colsAB]/=buffer;
397
+
398
+	  if (nFemales > 0)
399
+	    if (h == 2 && intInSet(&intbuffer, ptr2fIndex, &nFemales) > 0)
400
+	      likelihood[j] = likelihood[j+colsAB] = likelihood[j+2*colsAB] = 0;
401
+
402
+	  // IDENTICAL UNTIL HERE
403
+	  ptr2A[elem] = genotypeCall(&likelihood[j], &likelihood[j+colsAB], &likelihood[j+2*colsAB]);
404
+	  buffer = fmax(fmax(likelihood[j], likelihood[j+colsAB]), likelihood[j+2*colsAB]);
405
+	  if (ibSnpLevel1 == 1) buffer *= 0.995;
406
+	  if (ibSnpLevel2 == 1) buffer *= 0.000;
407
+	  ptr2B[elem] = genotypeConfidence(&buffer);
408
+	}
409
+	ibSnpLevel1 = ibSnpLevel2 = 0;
410
+      }
411
+  }
412
+  return(R_NilValue);
413
+}