forked from chrchang/plink-ng
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathplink_family.c
4791 lines (4713 loc) · 179 KB
/
plink_family.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#include "plink_common.h"
#include "plink_assoc.h"
#include "plink_cluster.h"
#include "plink_family.h"
#include "plink_stats.h"
void family_init(Family_info* fam_ip) {
fam_ip->mendel_max_trio_error = 1.0;
fam_ip->mendel_max_var_error = 1.0;
fam_ip->mendel_exclude_one_ratio = 0.0;
fam_ip->mendel_modifier = 0;
fam_ip->tdt_modifier = 0;
fam_ip->tdt_mperm_val = 0;
fam_ip->dfam_modifier = 0;
fam_ip->dfam_mperm_val = 0;
fam_ip->qfam_modifier = 0;
fam_ip->qfam_mperm_val = 0;
}
// bottom 2 bits of index = child genotype
// middle 2 bits of index = paternal genotype
// top 2 bits of index = maternal genotype
// bits 0-7 = child increment (always 1)
// bits 8-15 = father increment
// bits 16-23 = mother increment
// bits 24-31 = error code
const uint32_t mendel_error_table[] =
{0, 0, 0x1010101, 0x8000001,
0, 0, 0, 0x7010001,
0, 0, 0, 0x7010001,
0x3000101, 0, 0, 0x7010001,
0, 0, 0, 0x6000101,
0, 0, 0, 0,
0, 0, 0, 0,
0x3000101, 0, 0, 0,
0, 0, 0, 0x6000101,
0, 0, 0, 0,
0, 0, 0, 0,
0x3000101, 0, 0, 0,
0x4010001, 0, 0, 0x6000101,
0x4010001, 0, 0, 0,
0x4010001, 0, 0, 0,
0x5000001, 0, 0x2010101, 0};
// necessary to check child gender when dealing with error 9/10
const uint32_t mendel_error_table_x[] =
{0, 0, 0x1010101, 0x8000001,
0, 0, 0, 0x9010001,
0, 0, 0, 0x7010001,
0x3000101, 0, 0, 0x7010001,
0, 0, 0, 0x6000101,
0, 0, 0, 0,
0, 0, 0, 0,
0x3000101, 0, 0, 0,
0, 0, 0, 0x6000101,
0, 0, 0, 0,
0, 0, 0, 0,
0x3000101, 0, 0, 0,
0x4010001, 0, 0, 0x6000101,
0xa010001, 0, 0, 0,
0x4010001, 0, 0, 0,
0x5000001, 0, 0x2010101, 0};
int32_t get_trios_and_families(uintptr_t unfiltered_sample_ct, uintptr_t* sample_exclude, uintptr_t sample_ct, uintptr_t* founder_info, uintptr_t* sex_nm, uintptr_t* sex_male, char* sample_ids, uintptr_t max_sample_id_len, char* paternal_ids, uintptr_t max_paternal_id_len, char* maternal_ids, uintptr_t max_maternal_id_len, char** fids_ptr, uintptr_t* max_fid_len_ptr, char** iids_ptr, uintptr_t* max_iid_len_ptr, uint64_t** family_list_ptr, uint32_t* family_ct_ptr, uint64_t** trio_list_ptr, uintptr_t* trio_ct_ptr, uint32_t** trio_lookup_ptr, uint32_t include_duos, uint32_t toposort) {
// This mirrors linkRelateds() in genedrop.cpp, and parseTrios() in trio.cpp,
// in PLINK 1.07.
//
// family_list has paternal indices in low 32 bits, maternal indices in high
// 32, sorted in child ID order.
// trio_list has child IDs in low 32 bits, family_list indices in high 32
// bits, in PLINK 1.07 reporting order.
// If include_duos is set, missing parental IDs are coded as
// unfiltered_sample_ct. include_duos must be 0 or 1.
// If toposort is set, trio_lookup has child IDs in [4n], paternal IDs in
// [4n+1], maternal IDs in [4n+2], and reporting sequence index in [4n+3].
// Otherwise, it's [3n]/[3n+1]/[3n+2]. toposort must be 0 or 1.
//
// fids is a list of null-terminated FIDs using trio_list indices, and iids
// is a list of IIDs using regular unfiltered indices. If include_duos is
// set, iids has a trailing entry set to '0'. (fids_ptr, iids_ptr, and the
// corresponding lengths can be NULL.)
//
// PLINK 1.07 enforces <= 1 father and <= 1 mother per sample (and ambiguous
// sex parents are not permitted), but the IDs CAN be reversed in the .fam
// with no adverse consequences. For backward compatibility, we replicate
// this. (Possible todo: report a warning exactly once when this happens.)
// It won't be replicated in PLINK 2.0.
unsigned char* wkspace_mark = wkspace_base;
uint64_t* edge_list = NULL;
uint32_t* toposort_queue = NULL;
char* fids = NULL;
char* iids = NULL;
uintptr_t unfiltered_sample_ctl = (unfiltered_sample_ct + (BITCT - 1)) / BITCT;
uintptr_t unfiltered_sample_ctp1l = 1 + (unfiltered_sample_ct / BITCT);
uintptr_t sample_uidx = next_unset_unsafe(sample_exclude, 0);
// does *not* use populate_id_htable
uintptr_t htable_size = geqprime(2 * unfiltered_sample_ct + 1);
uintptr_t topsize = 0;
uintptr_t max_fid_len = 2;
uintptr_t max_iid_len = 2;
uint64_t family_code = 0;
uint32_t family_ct = 0;
uint32_t edge_ct = 0;
uint32_t remaining_edge_ct = 0;
uint32_t tqueue_start = 0;
uint32_t tqueue_end = 0;
int32_t retval = 0;
uintptr_t* founder_info2; // now decoupled from --make-founders
uint64_t* family_list;
uint64_t* trio_list_tmp;
// Initial "hash value" is twice the lower of the two parental uidxs. We use
// an open addressing scheme with prime table size >2n and a form of
// quadratic probing to ensure performance stays reasonable if someone has
// waaay too many concubines, etc. Specifically, given an initial hash
// collision at position 2j, and larger parent ID k (set to n if second
// parent is actually missing), next probe is (k+j) positions later, 3rd
// probe after that is 3(k+j) positions after the second, etc. (all numbers
// modulo table size).
uint64_t* family_htable;
uint64_t* trio_write;
uint64_t* edge_write;
uint32_t* family_idxs;
uint32_t* sample_id_map;
uint32_t* trio_lookup;
uint32_t* uiptr;
char* sorted_sample_ids;
char* idbuf;
char* idptr;
char* iidptr;
char* pidptr1;
char* pidptr2;
uint64_t trio_code;
uint64_t edge_code;
uint64_t ullii;
uintptr_t topsize_bak;
uintptr_t topsize_bak2;
uintptr_t family_idx;
uintptr_t trio_ct;
uintptr_t trio_idx;
uintptr_t fidlen;
uintptr_t sample_idx;
uintptr_t slen;
uintptr_t uidx1;
uintptr_t uidx2;
uintptr_t next_probe_incr;
uintptr_t probe_incr_incr;
uint32_t first_sex;
uint32_t uii;
int32_t sorted_idx;
founder_info2 = (uintptr_t*)top_alloc(&topsize, unfiltered_sample_ctp1l * sizeof(intptr_t));
if (!founder_info2) {
goto get_trios_and_families_ret_NOMEM;
}
memcpy(founder_info2, founder_info, unfiltered_sample_ctl * sizeof(intptr_t));
if (unfiltered_sample_ct & (BITCT - 1)) {
SET_BIT(founder_info2, unfiltered_sample_ct);
} else {
founder_info2[unfiltered_sample_ctl] = 1;
}
topsize_bak = topsize;
trio_list_tmp = (uint64_t*)top_alloc(&topsize, sample_ct * sizeof(int64_t));
if (!trio_list_tmp) {
goto get_trios_and_families_ret_NOMEM;
}
topsize_bak2 = topsize;
sorted_sample_ids = (char*)top_alloc(&topsize, sample_ct * max_sample_id_len);
if (!sorted_sample_ids) {
goto get_trios_and_families_ret_NOMEM;
}
sample_id_map = (uint32_t*)top_alloc(&topsize, sample_ct * sizeof(int32_t));
if (!sample_id_map) {
goto get_trios_and_families_ret_NOMEM;
}
idbuf = (char*)top_alloc(&topsize, max_sample_id_len);
if (!idbuf) {
goto get_trios_and_families_ret_NOMEM;
}
wkspace_left -= topsize;
if (sort_item_ids_noalloc(sorted_sample_ids, sample_id_map, unfiltered_sample_ct, sample_exclude, sample_ct, sample_ids, max_sample_id_len, 0, 0, strcmp_deref)) {
wkspace_left += topsize;
goto get_trios_and_families_ret_1;
}
// over-allocate here, we shrink family_list later when we know how many
// families there are
if (wkspace_alloc_ull_checked(&family_list, sample_ct * sizeof(int64_t)) ||
wkspace_alloc_ull_checked(&family_htable, htable_size * sizeof(int64_t)) ||
wkspace_alloc_ui_checked(&family_idxs, htable_size * sizeof(int32_t))) {
goto get_trios_and_families_ret_NOMEM2;
}
fill_ull_zero(family_htable, htable_size);
wkspace_left += topsize;
// 1. populate family_list (while using family_htable to track duplicates),
// determine max_iid_len, count qualifying trios
trio_write = trio_list_tmp;
for (sample_idx = 0; sample_idx < sample_ct; sample_uidx++, sample_idx++) {
next_unset_ul_unsafe_ck(sample_exclude, &sample_uidx);
idptr = &(sample_ids[sample_uidx * max_sample_id_len]);
iidptr = (char*)memchr(idptr, '\t', max_sample_id_len);
fidlen = (uintptr_t)((++iidptr) - idptr);
if (fidlen > max_fid_len) {
max_fid_len = fidlen;
}
slen = strlen(iidptr);
if (slen >= max_iid_len) {
max_iid_len = slen + 1;
}
if (IS_SET(founder_info, sample_uidx)) {
continue;
}
memcpy(idbuf, idptr, fidlen);
pidptr1 = &(paternal_ids[sample_uidx * max_paternal_id_len]);
pidptr2 = &(maternal_ids[sample_uidx * max_maternal_id_len]);
slen = strlen(pidptr1);
if (fidlen + slen < max_sample_id_len) {
memcpy(&(idbuf[fidlen]), pidptr1, slen + 1);
sorted_idx = bsearch_str(idbuf, fidlen + slen, sorted_sample_ids, max_sample_id_len, sample_ct);
} else {
sorted_idx = -1;
}
first_sex = 0;
if (sorted_idx == -1) {
if (!include_duos) {
SET_BIT(founder_info2, sample_uidx);
continue;
}
uidx1 = unfiltered_sample_ct;
} else {
uidx1 = sample_id_map[(uint32_t)sorted_idx];
if (uidx1 == sample_uidx) {
idbuf[fidlen - 1] = ' ';
LOGPREPRINTFWW("Error: '%s' is his/her own parent.\n", idbuf);
goto get_trios_and_families_ret_INVALID_FORMAT_2;
} else if (!IS_SET(sex_nm, uidx1)) {
idbuf[fidlen - 1] = ' ';
LOGPREPRINTFWW("Error: Parent '%s' has unspecified sex.\n", idbuf);
goto get_trios_and_families_ret_INVALID_FORMAT_2;
}
first_sex = 2 - IS_SET(sex_male, uidx1);
}
slen = strlen(pidptr2);
if (fidlen + slen < max_sample_id_len) {
memcpy(&(idbuf[fidlen]), pidptr2, slen);
sorted_idx = bsearch_str(idbuf, fidlen + slen, sorted_sample_ids, max_sample_id_len, sample_ct);
} else {
sorted_idx = -1;
}
if (sorted_idx == -1) {
if ((!include_duos) || (uidx1 == unfiltered_sample_ct)) {
SET_BIT(founder_info2, sample_uidx);
continue;
}
if (first_sex == 1) {
family_code = (((uint64_t)unfiltered_sample_ct) << 32) | ((uint64_t)uidx1);
} else {
family_code = (((uint64_t)uidx1) << 32) | ((uint64_t)unfiltered_sample_ct);
}
next_probe_incr = unfiltered_sample_ct + uidx1;
} else {
uidx2 = sample_id_map[(uint32_t)sorted_idx];
if (uidx2 == sample_uidx) {
idbuf[fidlen - 1] = ' ';
LOGPREPRINTFWW("Error: '%s' is their own parent.\n", idbuf);
goto get_trios_and_families_ret_INVALID_FORMAT_2;
} else if (!IS_SET(sex_nm, uidx2)) {
idbuf[fidlen - 1] = ' ';
LOGPREPRINTFWW("Error: Parent '%s' has unspecified sex.\n", idbuf);
goto get_trios_and_families_ret_INVALID_FORMAT_2;
}
uii = IS_SET(sex_male, uidx2);
if (2 - uii == first_sex) {
idptr[fidlen - 1] = ' ';
LOGPREPRINTFWW("Error: '%s' has two %sies.\n", idptr, (first_sex == 1)? "dadd" : "momm");
goto get_trios_and_families_ret_INVALID_FORMAT_2;
}
if (uii) {
family_code = (((uint64_t)uidx1) << 32) | ((uint64_t)uidx2);
} else {
family_code = (((uint64_t)uidx2) << 32) | ((uint64_t)uidx1);
}
next_probe_incr = uidx1 + uidx2;
if (uidx1 > uidx2) {
uidx1 = uidx2;
}
}
probe_incr_incr = next_probe_incr * 2;
uidx1 *= 2; // now current probe position
while (1) {
ullii = family_htable[uidx1];
if (!ullii) {
family_idx = family_ct++;
family_list[family_idx] = family_code;
family_htable[uidx1] = family_code;
family_idxs[uidx1] = family_idx;
break;
} else if (ullii == family_code) {
family_idx = family_idxs[uidx1];
break;
}
uidx1 += next_probe_incr;
if (uidx1 >= htable_size) {
uidx1 -= htable_size;
}
// quadratic probing implementation that avoids 64-bit modulus operation
// on 32-bit systems
next_probe_incr += probe_incr_incr;
if (next_probe_incr >= htable_size) {
next_probe_incr -= htable_size;
}
}
*trio_write++ = (((uint64_t)family_idx) << 32) | ((uint64_t)sample_uidx);
}
trio_ct = (uintptr_t)(trio_write - trio_list_tmp);
wkspace_reset(wkspace_mark);
wkspace_alloc(family_ct * sizeof(int64_t)); // family_list
topsize = topsize_bak2;
wkspace_left -= topsize;
if (wkspace_alloc_ull_checked(&trio_write, trio_ct * sizeof(int64_t))) {
goto get_trios_and_families_ret_NOMEM2;
}
wkspace_left += topsize;
topsize = topsize_bak;
memcpy(trio_write, trio_list_tmp, trio_ct * sizeof(int64_t));
#ifdef __cplusplus
std::sort((int64_t*)trio_write, (int64_t*)(&(trio_write[trio_ct])));
#else
qsort(trio_write, trio_ct, sizeof(int64_t), llcmp);
#endif
wkspace_left -= topsize;
if (wkspace_alloc_ui_checked(&trio_lookup, trio_ct * (3 + toposort) * sizeof(int32_t))) {
goto get_trios_and_families_ret_NOMEM2;
}
if (fids_ptr) {
if (wkspace_alloc_c_checked(&fids, trio_ct * max_fid_len) ||
wkspace_alloc_c_checked(&iids, (unfiltered_sample_ct + include_duos) * max_iid_len)) {
goto get_trios_and_families_ret_NOMEM2;
}
}
if (toposort) {
if (wkspace_alloc_ull_checked(&edge_list, trio_ct * 2 * sizeof(int64_t))) {
goto get_trios_and_families_ret_NOMEM2;
}
// Edge list excludes founder parents; edge codes have parental uidx in
// high 32 bits, trio idx in low 32.
edge_write = edge_list;
for (trio_idx = 0; trio_idx < trio_ct; trio_idx++) {
family_code = family_list[(uintptr_t)(trio_write[trio_idx] >> 32)];
if (!is_set(founder_info2, family_code >> 32)) {
*edge_write++ = (family_code & 0xffffffff00000000LLU) | ((uint64_t)trio_idx);
}
if (!is_set(founder_info2, (uint32_t)family_code)) {
*edge_write++ = ((family_code & 0xffffffffLLU) << 32) | ((uint64_t)trio_idx);
}
}
edge_ct = (uintptr_t)(edge_write - edge_list);
if (edge_ct) {
#ifdef __cplusplus
std::sort((int64_t*)edge_list, (int64_t*)(&(edge_list[edge_ct])));
#else
qsort(edge_list, edge_ct, sizeof(int64_t), llcmp);
#endif
}
wkspace_shrink_top(edge_list, edge_ct * sizeof(int64_t));
if (wkspace_alloc_ui_checked(&toposort_queue, trio_ct * sizeof(int32_t))) {
goto get_trios_and_families_ret_NOMEM2;
}
remaining_edge_ct = edge_ct;
}
wkspace_left += topsize;
*family_list_ptr = family_list;
*family_ct_ptr = family_ct;
*trio_list_ptr = trio_write;
*trio_ct_ptr = trio_ct;
*trio_lookup_ptr = trio_lookup;
if (max_fid_len_ptr) {
*max_fid_len_ptr = max_fid_len;
}
if (fids_ptr) {
*fids_ptr = fids;
*iids_ptr = iids;
*max_iid_len_ptr = max_iid_len;
sample_uidx = next_unset_unsafe(sample_exclude, 0);
for (sample_idx = 0; sample_idx < sample_ct; sample_uidx++, sample_idx++) {
next_unset_ul_unsafe_ck(sample_exclude, &sample_uidx);
idptr = &(sample_ids[sample_uidx * max_sample_id_len]);
iidptr = (char*)memchr(idptr, '\t', max_fid_len);
strcpy(&(iids[sample_uidx * max_iid_len]), &(iidptr[1]));
}
if (include_duos) {
memcpy(&(iids[unfiltered_sample_ct * max_iid_len]), "0", 2);
}
}
uiptr = trio_lookup;
for (trio_idx = 0; trio_idx < trio_ct; trio_idx++) {
trio_code = trio_write[trio_idx];
sample_uidx = (uint32_t)trio_code;
if (fids_ptr) {
idptr = &(sample_ids[sample_uidx * max_sample_id_len]);
iidptr = (char*)memchr(idptr, '\t', max_fid_len);
fidlen = (uintptr_t)(iidptr - idptr);
memcpyx(&(fids[trio_idx * max_fid_len]), idptr, fidlen, '\0');
}
family_code = family_list[(uintptr_t)(trio_code >> 32)];
if (!toposort) {
*uiptr++ = sample_uidx;
*uiptr++ = (uint32_t)family_code;
*uiptr++ = (uint32_t)(family_code >> 32);
} else if (is_set(founder_info2, family_code >> 32) && is_set(founder_info2, (uint32_t)family_code)) {
// Just populate the "no incoming edges" queue here.
toposort_queue[tqueue_end++] = trio_idx;
}
}
if (toposort) {
// Kahn A. B. (1962) Topological sorting of large networks.
// Communications of the ACM, 5.
// This is a breadth-first implementation.
while (tqueue_start < tqueue_end) {
trio_idx = toposort_queue[tqueue_start++];
trio_code = trio_write[trio_idx];
sample_uidx = (uint32_t)trio_code;
family_code = family_list[(uintptr_t)(trio_code >> 32)];
*uiptr++ = sample_uidx;
*uiptr++ = (uint32_t)family_code;
*uiptr++ = (uint32_t)(family_code >> 32);
*uiptr++ = trio_idx;
if (remaining_edge_ct) {
SET_BIT(founder_info2, sample_uidx);
ullii = ((uint64_t)sample_uidx) << 32;
uii = uint64arr_greater_than(edge_list, edge_ct, ullii);
ullii |= 0xffffffffLLU;
while (uii < edge_ct) {
edge_code = edge_list[uii];
if (edge_code > ullii) {
break;
}
remaining_edge_ct--;
// look up child, see if other parent is now a founder
trio_code = trio_write[(uint32_t)edge_code];
family_code = family_list[(uintptr_t)(trio_code >> 32)];
if (is_set(founder_info2, family_code >> 32) && is_set(founder_info2, (uint32_t)family_code)) {
toposort_queue[tqueue_end++] = (uint32_t)edge_code;
}
uii++;
}
}
}
if (remaining_edge_ct) {
logprint("Error: Pedigree graph is cyclic. Check for evidence of time travel abuse in\nyour cohort.\n");
goto get_trios_and_families_ret_INVALID_FORMAT;
}
wkspace_reset(edge_list);
}
while (0) {
get_trios_and_families_ret_NOMEM2:
wkspace_left += topsize;
get_trios_and_families_ret_NOMEM:
retval = RET_NOMEM;
break;
get_trios_and_families_ret_INVALID_FORMAT_2:
logprintb();
get_trios_and_families_ret_INVALID_FORMAT:
retval = RET_INVALID_FORMAT;
break;
}
get_trios_and_families_ret_1:
if (retval) {
wkspace_reset(wkspace_mark);
}
return retval;
}
uint32_t erase_mendel_errors(uintptr_t unfiltered_sample_ct, uintptr_t* loadbuf, uintptr_t* workbuf, uint32_t* trio_lookup, uint32_t trio_ct, uint32_t multigen) {
uint32_t* uiptr = trio_lookup;
uint32_t cur_errors = 0;
uint32_t trio_idx;
uint32_t lookup_idx;
uintptr_t ulii;
uintptr_t uljj;
uint32_t uii;
uint32_t ujj;
uint32_t ukk;
uint32_t umm;
uint32_t unn;
uint32_t uoo;
uint32_t upp;
memcpy(workbuf, loadbuf, (unfiltered_sample_ct + 3) / 4);
SET_BIT_DBL(workbuf, unfiltered_sample_ct);
if (!multigen) {
for (trio_idx = 0; trio_idx < trio_ct; trio_idx++) {
uii = *uiptr++;
ujj = *uiptr++;
ukk = *uiptr++;
umm = mendel_error_table[((workbuf[uii / BITCT2] >> (2 * (uii % BITCT2))) & 3) | (((workbuf[ujj / BITCT2] >> (2 * (ujj % BITCT2))) & 3) << 2) | (((workbuf[ukk / BITCT2] >> (2 * (ukk % BITCT2))) & 3) << 4)];
if (umm) {
ulii = loadbuf[uii / BITCT2];
uljj = ONELU << (2 * (uii % BITCT2));
loadbuf[uii / BITCT2] = (ulii & (~(3 * uljj))) | uljj;
if (umm & 0x100) {
ulii = loadbuf[ujj / BITCT2];
uljj = ONELU << (2 * (ujj % BITCT2));
loadbuf[ujj / BITCT2] = (ulii & (~(3 * uljj))) | uljj;
}
if (umm & 0x10000) {
ulii = loadbuf[ukk / BITCT2];
uljj = ONELU << (2 * (ukk % BITCT2));
loadbuf[ukk / BITCT2] = (ulii & (~(3 * uljj))) | uljj;
}
cur_errors++;
}
}
} else {
for (lookup_idx = 0; lookup_idx < trio_ct; lookup_idx++) {
uii = *uiptr++;
ujj = *uiptr++;
ukk = *uiptr++;
trio_idx = *uiptr++;
unn = uii / BITCT2;
uoo = ujj / BITCT2;
upp = ukk / BITCT2;
uii = 2 * (uii % BITCT2);
ujj = 2 * (ujj % BITCT2);
ukk = 2 * (ukk % BITCT2);
ulii = (workbuf[unn] >> uii) & 3;
uljj = ((workbuf[uoo] >> ujj) & 3) | (((workbuf[upp] >> ukk) & 3) << 2);
if (ulii != 1) {
umm = mendel_error_table[ulii | (uljj << 2)];
if (umm) {
ulii = loadbuf[unn];
uljj = ONELU << uii;
loadbuf[unn] = (ulii & (~(3 * uljj))) | uljj;
if (umm & 0x100) {
ulii = loadbuf[uoo];
uljj = ONELU << ujj;
loadbuf[uoo] = (ulii & (~(3 * uljj))) | uljj;
}
if (umm & 0x10000) {
ulii = loadbuf[upp];
uljj = ONELU << ukk;
loadbuf[upp] = (ulii & (~(3 * uljj))) | uljj;
}
cur_errors++;
}
} else if (!uljj) {
workbuf[unn] &= ~(ONELU << uii);
} else if (uljj == 15) {
workbuf[unn] |= (2 * ONELU) << uii;
}
}
}
return cur_errors;
}
void fill_mendel_errstr(uint32_t error_code, char** allele_ptrs, uint32_t* alens, char* wbuf, uint32_t* len_ptr) {
char* wptr;
uint32_t len;
if (error_code < 10) {
wbuf[0] = ' ';
wbuf[1] = error_code + '0';
} else {
memcpy(wbuf, "10", 2);
}
if (!alens[0]) {
// lazy fill
alens[0] = strlen(allele_ptrs[0]);
alens[1] = strlen(allele_ptrs[1]);
}
// PLINK 1.07 may fail to put space here when there are long allele codes; we
// don't replicate that. (We actually force two spaces here since the error
// column itself contains spaces.)
wptr = memseta(&(wbuf[2]), 32, 2);
switch (error_code) {
case 1:
len = 5 * alens[0] + alens[1] + 10;
if (len < 20) {
wptr = memseta(wptr, 32, 20 - len);
}
wptr = memcpyax(wptr, allele_ptrs[0], alens[0], '/');
wptr = memcpya(wptr, allele_ptrs[0], alens[0]);
wptr = memcpyl3a(wptr, " x ");
wptr = memcpyax(wptr, allele_ptrs[0], alens[0], '/');
wptr = memcpya(wptr, allele_ptrs[0], alens[0]);
wptr = memcpya(wptr, " -> ", 4);
wptr = memcpyax(wptr, allele_ptrs[0], alens[0], '/');
wptr = memcpya(wptr, allele_ptrs[1], alens[1]);
break;
case 2:
len = alens[0] + 5 * alens[1] + 10;
if (len < 20) {
wptr = memseta(wptr, 32, 20 - len);
}
wptr = memcpyax(wptr, allele_ptrs[1], alens[1], '/');
wptr = memcpya(wptr, allele_ptrs[1], alens[1]);
wptr = memcpyl3a(wptr, " x ");
wptr = memcpyax(wptr, allele_ptrs[1], alens[1], '/');
wptr = memcpya(wptr, allele_ptrs[1], alens[1]);
wptr = memcpya(wptr, " -> ", 4);
wptr = memcpyax(wptr, allele_ptrs[0], alens[0], '/');
wptr = memcpya(wptr, allele_ptrs[1], alens[1]);
break;
case 3:
len = 2 * alens[0] + 2 * alens[1] + 12;
if (len < 20) {
wptr = memseta(wptr, 32, 20 - len);
}
wptr = memcpyax(wptr, allele_ptrs[1], alens[1], '/');
wptr = memcpya(wptr, allele_ptrs[1], alens[1]);
wptr = memcpya(wptr, " x */* -> ", 10);
wptr = memcpyax(wptr, allele_ptrs[0], alens[0], '/');
wptr = memcpya(wptr, allele_ptrs[0], alens[0]);
break;
case 4:
case 10:
len = 2 * alens[0] + 2 * alens[1] + 12;
if (len < 20) {
wptr = memseta(wptr, 32, 20 - len);
}
wptr = memcpya(wptr, "*/* x ", 6);
wptr = memcpyax(wptr, allele_ptrs[1], alens[1], '/');
wptr = memcpya(wptr, allele_ptrs[1], alens[1]);
wptr = memcpya(wptr, " -> ", 4);
wptr = memcpyax(wptr, allele_ptrs[0], alens[0], '/');
wptr = memcpya(wptr, allele_ptrs[0], alens[0]);
break;
case 5:
len = 2 * alens[0] + 4 * alens[1] + 10;
if (len < 20) {
wptr = memseta(wptr, 32, 20 - len);
}
wptr = memcpyax(wptr, allele_ptrs[1], alens[1], '/');
wptr = memcpya(wptr, allele_ptrs[1], alens[1]);
wptr = memcpyl3a(wptr, " x ");
wptr = memcpyax(wptr, allele_ptrs[1], alens[1], '/');
wptr = memcpya(wptr, allele_ptrs[1], alens[1]);
wptr = memcpya(wptr, " -> ", 4);
wptr = memcpyax(wptr, allele_ptrs[0], alens[0], '/');
wptr = memcpya(wptr, allele_ptrs[0], alens[0]);
break;
case 6:
len = 2 * alens[0] + 2 * alens[1] + 12;
if (len < 20) {
wptr = memseta(wptr, 32, 20 - len);
}
wptr = memcpyax(wptr, allele_ptrs[0], alens[0], '/');
wptr = memcpya(wptr, allele_ptrs[0], alens[0]);
wptr = memcpya(wptr, " x */* -> ", 10);
wptr = memcpyax(wptr, allele_ptrs[1], alens[1], '/');
wptr = memcpya(wptr, allele_ptrs[1], alens[1]);
break;
case 7:
case 9:
len = 2 * alens[0] + 2 * alens[1] + 12;
if (len < 20) {
wptr = memseta(wptr, 32, 20 - len);
}
wptr = memcpya(wptr, "*/* x ", 6);
wptr = memcpyax(wptr, allele_ptrs[0], alens[0], '/');
wptr = memcpya(wptr, allele_ptrs[0], alens[0]);
wptr = memcpya(wptr, " -> ", 4);
wptr = memcpyax(wptr, allele_ptrs[1], alens[1], '/');
wptr = memcpya(wptr, allele_ptrs[1], alens[1]);
break;
case 8:
len = 4 * alens[0] + 2 * alens[1] + 10;
if (len < 20) {
wptr = memseta(wptr, 32, 20 - len);
}
wptr = memcpyax(wptr, allele_ptrs[0], alens[0], '/');
wptr = memcpya(wptr, allele_ptrs[0], alens[0]);
wptr = memcpyl3a(wptr, " x ");
wptr = memcpyax(wptr, allele_ptrs[0], alens[0], '/');
wptr = memcpya(wptr, allele_ptrs[0], alens[0]);
wptr = memcpya(wptr, " -> ", 4);
wptr = memcpyax(wptr, allele_ptrs[1], alens[1], '/');
wptr = memcpya(wptr, allele_ptrs[1], alens[1]);
break;
}
*wptr++ = '\n';
*len_ptr = (uintptr_t)(wptr - wbuf);
}
int32_t mendel_error_scan(Family_info* fam_ip, FILE* bedfile, uintptr_t bed_offset, char* outname, char* outname_end, uint32_t plink_maxfid, uint32_t plink_maxiid, uint32_t plink_maxsnp, uintptr_t unfiltered_marker_ct, uintptr_t* marker_exclude, uintptr_t* marker_exclude_ct_ptr, uintptr_t* marker_reverse, char* marker_ids, uintptr_t max_marker_id_len, char** marker_allele_ptrs, uintptr_t max_marker_allele_len, uintptr_t unfiltered_sample_ct, uintptr_t* sample_exclude, uintptr_t* sample_exclude_ct_ptr, uintptr_t* founder_info, uintptr_t* sex_nm, uintptr_t* sex_male, char* sample_ids, uintptr_t max_sample_id_len, char* paternal_ids, uintptr_t max_paternal_id_len, char* maternal_ids, uintptr_t max_maternal_id_len, uint32_t hh_exists, Chrom_info* chrom_info_ptr, uint32_t calc_mendel) {
unsigned char* wkspace_mark = wkspace_base;
FILE* outfile = NULL;
FILE* outfile_l = NULL;
uintptr_t* sample_male_include2 = NULL;
uintptr_t* error_locs = NULL;
char* varptr = NULL;
char* chrom_name_ptr = NULL;
unsigned char* cur_errors = NULL;
uint64_t* family_error_cts = NULL;
uint32_t* child_cts = NULL;
uintptr_t marker_ct = unfiltered_marker_ct - *marker_exclude_ct_ptr;
uintptr_t unfiltered_sample_ct4 = (unfiltered_sample_ct + 3) / 4;
uintptr_t unfiltered_sample_ctp1l2 = 1 + (unfiltered_sample_ct / BITCT2);
uintptr_t final_mask = get_final_mask(unfiltered_sample_ct);
uintptr_t sample_ct = unfiltered_sample_ct - *sample_exclude_ct_ptr;
uintptr_t marker_uidx = ~ZEROLU;
uint64_t tot_error_ct = 0;
uint32_t unfiltered_sample_ctl2m1 = (unfiltered_sample_ct - 1) / BITCT2;
uint32_t include_duos = (fam_ip->mendel_modifier / MENDEL_DUOS) & 1;
uint32_t multigen = (fam_ip->mendel_modifier / MENDEL_MULTIGEN) & 1;
uint32_t var_first = fam_ip->mendel_modifier & MENDEL_FILTER_VAR_FIRST;
uint32_t varlen = 0;
uint32_t chrom_name_len = 0;
uint32_t new_marker_exclude_ct = 0;
uint32_t error_ct_fill = 0;
int32_t retval = 0;
char chrom_name_buf[5];
char* errstrs[10];
uint32_t errstr_lens[11];
uint32_t alens[2];
const uint32_t* error_table_ptr;
char* fids;
char* iids;
char* wptr;
char* cptr;
uint64_t* family_list;
uint64_t* trio_list;
uintptr_t* loadbuf;
uint32_t* trio_lookup;
uint32_t* error_cts;
uint32_t* error_cts_tmp;
uint32_t* error_cts_tmp2;
uint32_t* uiptr;
#ifdef __LP64__
__m128i* vptr;
__m128i* vptr2;
#endif
uintptr_t trio_ct4;
uintptr_t max_fid_len;
uintptr_t max_iid_len;
uintptr_t trio_ct;
uintptr_t trio_ctl;
uintptr_t trio_idx;
uintptr_t lookup_idx;
uintptr_t ulii;
uintptr_t uljj;
double exclude_one_ratio;
double dxx;
double dyy;
uint64_t trio_code;
uint64_t family_code;
uint32_t family_ct;
uint32_t var_error_max;
uint32_t cur_error_ct;
uint32_t chrom_fo_idx;
uint32_t chrom_idx;
uint32_t chrom_end;
uint32_t is_x;
uint32_t uii;
uint32_t ujj;
uint32_t ukk;
uint32_t umm;
marker_ct -= count_non_autosomal_markers(chrom_info_ptr, marker_exclude, 0, 1);
if ((!marker_ct) || is_set(chrom_info_ptr->haploid_mask, 0)) {
logprint("Warning: Skipping --me/--mendel since there is no autosomal or Xchr data.\n");
goto mendel_error_scan_ret_1;
}
retval = get_trios_and_families(unfiltered_sample_ct, sample_exclude, sample_ct, founder_info, sex_nm, sex_male, sample_ids, max_sample_id_len, paternal_ids, max_paternal_id_len, maternal_ids, max_maternal_id_len, &fids, &max_fid_len, &iids, &max_iid_len, &family_list, &family_ct, &trio_list, &trio_ct, &trio_lookup, include_duos, multigen);
if (retval) {
goto mendel_error_scan_ret_1;
}
if (!trio_ct) {
LOGPRINTF("Warning: Skipping --me/--mendel since there are no %strios.\n", include_duos? "duos or " : "");
goto mendel_error_scan_ret_1;
}
trio_ct4 = (trio_ct + 3) / 4;
trio_ctl = (trio_ct + (BITCT - 1)) / BITCT;
var_error_max = (int32_t)(fam_ip->mendel_max_var_error * (1 + SMALL_EPSILON) * ((intptr_t)trio_ct));
if (wkspace_alloc_ul_checked(&loadbuf, unfiltered_sample_ctp1l2 * sizeof(intptr_t)) ||
wkspace_alloc_ui_checked(&error_cts, trio_ct * 3 * sizeof(int32_t)) ||
wkspace_alloc_ui_checked(&error_cts_tmp, trio_ct4 * 4 * sizeof(int32_t))) {
goto mendel_error_scan_ret_NOMEM;
}
if (!var_first) {
error_cts_tmp2 = error_cts_tmp;
} else {
if (wkspace_alloc_ui_checked(&error_cts_tmp2, trio_ct4 * 4 * sizeof(int32_t))) {
goto mendel_error_scan_ret_NOMEM;
}
fill_uint_zero(error_cts_tmp2, trio_ct4 * 4);
}
loadbuf[unfiltered_sample_ctp1l2 - 1] = 0;
fill_uint_zero(error_cts, trio_ct * 3);
fill_uint_zero(error_cts_tmp, trio_ct4 * 4);
hh_exists &= XMHH_EXISTS;
if (alloc_raw_haploid_filters(unfiltered_sample_ct, hh_exists, 0, sample_exclude, sex_male, NULL, &sample_male_include2)) {
goto mendel_error_scan_ret_NOMEM;
}
alens[0] = 0;
alens[1] = 0;
if (calc_mendel) {
// ugh, this calculation was totally off before.
// * max_marker_allele_len includes trailing null (though forgetting this
// was harmless)
// * minimum buffer size was 25, not 21, due to inclusion of error code and
// following double-space in the string (forgetting this was NOT
// harmless)
ulii = max_marker_allele_len * 6 + 9;
if (ulii < 25) {
ulii = 25;
}
if (wkspace_alloc_ull_checked(&family_error_cts, family_ct * 3 * sizeof(int64_t)) ||
wkspace_alloc_ui_checked(&child_cts, family_ct * sizeof(int32_t)) ||
wkspace_alloc_c_checked(&(errstrs[0]), ulii * 10)) {
goto mendel_error_scan_ret_NOMEM;
}
for (uii = 1; uii < 10; uii++) {
errstrs[uii] = &(errstrs[0][uii * ulii]);
}
if (multigen) {
if (wkspace_alloc_ul_checked(&error_locs, trio_ctl * sizeof(intptr_t)) ||
wkspace_alloc_uc_checked(&cur_errors, trio_ct)) {
goto mendel_error_scan_ret_NOMEM;
}
fill_ulong_zero(error_locs, trio_ctl);
}
memcpy(outname_end, ".mendel", 8);
if (fopen_checked(&outfile, outname, "w")) {
goto mendel_error_scan_ret_OPEN_FAIL;
}
sprintf(tbuf, "%%%us %%%us CHR %%%us CODE ERROR\n", plink_maxfid, plink_maxiid, plink_maxsnp);
fprintf(outfile, tbuf, "FID", "KID", "SNP");
memcpy(outname_end, ".lmendel", 9);
if (fopen_checked(&outfile_l, outname, "w")) {
goto mendel_error_scan_ret_OPEN_FAIL;
}
// replicate harmless 'N' misalignment bug
sprintf(tbuf, " CHR %%%us N\n", plink_maxsnp);
fprintf(outfile_l, tbuf, "SNP");
} else {
// suppress warning
fill_ulong_zero((uintptr_t*)errstrs, 10);
}
for (chrom_fo_idx = 0; chrom_fo_idx < chrom_info_ptr->chrom_ct; chrom_fo_idx++) {
chrom_idx = chrom_info_ptr->chrom_file_order[chrom_fo_idx];
is_x = (((uint32_t)chrom_info_ptr->x_code) == chrom_idx);
if ((IS_SET(chrom_info_ptr->haploid_mask, chrom_idx) && (!is_x)) || (((uint32_t)chrom_info_ptr->mt_code) == chrom_idx)) {
continue;
}
chrom_end = chrom_info_ptr->chrom_file_order_marker_idx[chrom_fo_idx + 1];
uii = next_unset(marker_exclude, chrom_info_ptr->chrom_file_order_marker_idx[chrom_fo_idx], chrom_end);
if (uii == chrom_end) {
continue;
}
if (!is_x) {
error_table_ptr = mendel_error_table;
} else {
error_table_ptr = mendel_error_table_x;
}
if (calc_mendel) {
chrom_name_ptr = chrom_name_buf5w4write(chrom_name_buf, chrom_info_ptr, chrom_idx, &chrom_name_len);
}
if (uii != marker_uidx) {
marker_uidx = uii;
goto mendel_error_scan_seek;
}
while (1) {
if (load_raw2(bedfile, loadbuf, unfiltered_sample_ct4, unfiltered_sample_ctl2m1, final_mask)) {
goto mendel_error_scan_ret_READ_FAIL;
}
if (IS_SET(marker_reverse, marker_uidx)) {
reverse_loadbuf((unsigned char*)loadbuf, unfiltered_sample_ct);
}
if (hh_exists && is_x) {
hh_reset((unsigned char*)loadbuf, sample_male_include2, unfiltered_sample_ct);
}
// missing parents are treated as having uidx equal to
// unfiltered_sample_ct, and we set the corresponding genotype to always
// be missing. This lets us avoid special-casing duos.
SET_BIT_DBL(loadbuf, unfiltered_sample_ct);
uiptr = trio_lookup;
cur_error_ct = 0;
if (calc_mendel) {
varptr = &(marker_ids[marker_uidx * max_marker_id_len]);
varlen = strlen(varptr);
alens[0] = 0;
alens[1] = 0;
fill_uint_zero(errstr_lens, 11);
}
if (!multigen) {
for (trio_idx = 0; trio_idx < trio_ct; trio_idx++) {
uii = *uiptr++;
ujj = *uiptr++;
ukk = *uiptr++;
umm = error_table_ptr[((loadbuf[uii / BITCT2] >> (2 * (uii % BITCT2))) & 3) | (((loadbuf[ujj / BITCT2] >> (2 * (ujj % BITCT2))) & 3) << 2) | (((loadbuf[ukk / BITCT2] >> (2 * (ukk % BITCT2))) & 3) << 4)];
if (umm) {
error_cts_tmp2[trio_idx] += umm & 0xffffff;
cur_error_ct++;
if (calc_mendel) {
umm >>= 24;
if ((umm > 8) && (!is_set(sex_male, uii))) {
umm = 34 - 3 * umm; // 9 -> 7, 10 -> 4
}
wptr = fw_strcpy(plink_maxfid, &(fids[trio_idx * max_fid_len]), tbuf);
*wptr++ = ' ';
wptr = fw_strcpy(plink_maxiid, &(iids[uii * max_iid_len]), wptr);
*wptr++ = ' ';
wptr = memcpyax(wptr, chrom_name_ptr, chrom_name_len, ' ');
wptr = fw_strcpyn(plink_maxsnp, varlen, varptr, wptr);
wptr = memseta(wptr, 32, 5);
if (fwrite_checked(tbuf, wptr - tbuf, outfile)) {
goto mendel_error_scan_ret_WRITE_FAIL;
}
if (!errstr_lens[umm]) {
fill_mendel_errstr(umm, &(marker_allele_ptrs[2 * marker_uidx]), alens, errstrs[umm - 1], &(errstr_lens[umm]));
}
if (fwrite_checked(errstrs[umm - 1], errstr_lens[umm], outfile)) {
goto mendel_error_scan_ret_WRITE_FAIL;
}
}
}
}
} else {
for (lookup_idx = 0; lookup_idx < trio_ct; lookup_idx++) {
uii = *uiptr++;
ujj = *uiptr++;
ukk = *uiptr++;
trio_idx = *uiptr++;
uljj = ((loadbuf[ujj / BITCT2] >> (2 * (ujj % BITCT2))) & 3) | (((loadbuf[ukk / BITCT2] >> (2 * (ukk % BITCT2))) & 3) << 2);
umm = uii / BITCT2;
ujj = 2 * (uii % BITCT2);
ulii = (loadbuf[umm] >> ujj) & 3;
if (ulii != 1) {
umm = error_table_ptr[ulii | (uljj << 2)];
if (umm) {
error_cts_tmp2[trio_idx] += umm & 0xffffff;
cur_error_ct++;
if (calc_mendel) {
set_bit(error_locs, trio_idx);
umm >>= 24;
if ((umm > 8) && (!is_set(sex_male, uii))) {
umm = 34 - 3 * umm;
}
cur_errors[trio_idx] = (unsigned char)umm;
}
}
} else if (!uljj) {
loadbuf[umm] &= ~(ONELU << ujj);
} else if (uljj == 15) {
loadbuf[umm] |= (2 * ONELU) << ujj;
}
}
if (calc_mendel && cur_error_ct) {
trio_idx = 0;
for (uii = 0; uii < cur_error_ct; trio_idx++, uii++) {
next_set_ul_unsafe_ck(error_locs, &trio_idx);
wptr = fw_strcpy(plink_maxfid, &(fids[trio_idx * max_fid_len]), tbuf);
*wptr++ = ' ';
wptr = fw_strcpy(plink_maxiid, &(iids[((uint32_t)trio_list[trio_idx]) * max_iid_len]), wptr);
*wptr++ = ' ';
wptr = memcpyax(wptr, chrom_name_ptr, chrom_name_len, ' ');
wptr = fw_strcpyn(plink_maxsnp, varlen, varptr, wptr);
wptr = memseta(wptr, 32, 5);
if (fwrite_checked(tbuf, wptr - tbuf, outfile)) {
goto mendel_error_scan_ret_WRITE_FAIL;
}
umm = cur_errors[trio_idx];
if (!errstr_lens[umm]) {
fill_mendel_errstr(umm, &(marker_allele_ptrs[2 * marker_uidx]), alens, errstrs[umm - 1], &(errstr_lens[umm]));
}
if (fwrite_checked(errstrs[umm - 1], errstr_lens[umm], outfile)) {
goto mendel_error_scan_ret_WRITE_FAIL;
}
}
fill_ulong_zero(error_locs, trio_ctl);
}
}
if (calc_mendel) {
if (fwrite_checked(chrom_name_ptr, chrom_name_len, outfile_l)) {
goto mendel_error_scan_ret_WRITE_FAIL;
}
tbuf[0] = ' ';
wptr = fw_strcpyn(plink_maxsnp, varlen, varptr, &(tbuf[1]));
*wptr++ = ' ';
wptr = uint32_writew4x(wptr, cur_error_ct, '\n');
if (fwrite_checked(tbuf, wptr - tbuf, outfile_l)) {
goto mendel_error_scan_ret_WRITE_FAIL;
}
}
if (cur_error_ct) {
if (cur_error_ct > var_error_max) {
SET_BIT(marker_exclude, marker_uidx);
new_marker_exclude_ct++;
}
if ((cur_error_ct <= var_error_max) || (!var_first)) {
if (var_first) {
#ifdef __LP64__
vptr = (__m128i*)error_cts_tmp;
vptr2 = (__m128i*)error_cts_tmp2;
for (trio_idx = 0; trio_idx < trio_ct4; trio_idx++) {
*vptr = _mm_add_epi64(*vptr, *vptr2++);
vptr++;
}
#else
for (trio_idx = 0; trio_idx < trio_ct; trio_idx++) {
error_cts_tmp[trio_idx] += error_cts_tmp2[trio_idx];
}