Nie możesz wybrać więcej, niż 25 tematów Tematy muszą się zaczynać od litery lub cyfry, mogą zawierać myślniki ('-') i mogą mieć do 35 znaków.
 
 
 
 

150 wiersze
5.2 KiB

  1. /* call.h -- variant calling declarations.
  2. Copyright (C) 2013-2015, 2019-2020 Genome Research Ltd.
  3. Author: Petr Danecek <pd3@sanger.ac.uk>
  4. Permission is hereby granted, free of charge, to any person obtaining a copy
  5. of this software and associated documentation files (the "Software"), to deal
  6. in the Software without restriction, including without limitation the rights
  7. to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  8. copies of the Software, and to permit persons to whom the Software is
  9. furnished to do so, subject to the following conditions:
  10. The above copyright notice and this permission notice shall be included in
  11. all copies or substantial portions of the Software.
  12. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  13. IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  14. FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  15. AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  16. LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  17. OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  18. THE SOFTWARE. */
  19. #ifndef __CALL_H__
  20. #define __CALL_H__
  21. #include <htslib/vcf.h>
  22. #include <htslib/synced_bcf_reader.h>
  23. #include "vcmp.h"
  24. #define CALL_KEEPALT 1
  25. #define CALL_VARONLY (1<<1)
  26. #define CALL_CONSTR_TRIO (1<<2)
  27. #define CALL_CONSTR_ALLELES (1<<3)
  28. //
  29. #define CALL_FMT_PV4 (1<<5)
  30. #define CALL_FMT_GQ (1<<6)
  31. #define CALL_FMT_GP (1<<7)
  32. #define FATHER 0
  33. #define MOTHER 1
  34. #define CHILD 2
  35. typedef struct
  36. {
  37. char *name;
  38. int sample[3]; // father, mother, child
  39. int type; // see FTYPE_* definitions in mcall.c
  40. }
  41. family_t;
  42. // For the single-sample and grouped -G calling
  43. typedef struct
  44. {
  45. double ref_lk, max_lk, lk_sum;
  46. float *qsum; // QS(quality sum) values
  47. int nqsum;
  48. uint32_t *smpl, nsmpl;
  49. uint32_t nals, als;
  50. }
  51. smpl_grp_t;
  52. // For the `-C alleles -i` constrained calling
  53. typedef struct
  54. {
  55. uint32_t n:31, used:1;
  56. char **allele;
  57. }
  58. tgt_als_t;
  59. typedef struct _ccall_t ccall_t;
  60. typedef struct
  61. {
  62. // mcall only
  63. int npdg;
  64. int *als_map, nals_map; // mapping from full set of alleles to trimmed set of alleles (old -> new)
  65. int *pl_map, npl_map; // same as above for PLs, but reverse (new -> old)
  66. char **als; // array to hold the trimmed set of alleles to appear on output
  67. int nals; // size of the als array
  68. int als_new, nals_new; // bitmask with final alleles and their number
  69. family_t *fams; // list of families and samples for trio calling
  70. int nfams, mfams;
  71. int ntrio[5][5]; // possible trio genotype combinations and their counts; first idx:
  72. uint16_t *trio[5][5]; // family type, second index: allele count (2-4, first two are unused)
  73. double *GLs;
  74. float *GPs; // FORMAT/GP: posterior probabilities
  75. int32_t *GQs, *ADs; // FORMAT/GQ: genotype qualities; AD: allelic depth for -G
  76. int32_t *itmp; // temporary int array, used for new PLs with CALL_CONSTR_ALLELES
  77. int n_itmp, nGPs, nADs;
  78. vcmp_t *vcmp;
  79. double trio_Pm_SNPs, trio_Pm_del, trio_Pm_ins; // P(mendelian) for trio calling, see mcall_call_trio_genotypes()
  80. int32_t *ugts, *cgts; // unconstraind and constrained GTs
  81. uint32_t output_tags;
  82. char *prior_AN, *prior_AC; // reference panel AF tags (AF=AC/AN)
  83. tgt_als_t *tgt_als; // for CALL_CONSTR_ALLELES
  84. char *sample_groups; // for single-sample or grouped calling with -G
  85. char *sample_groups_tag; // for -G [AD|QS:]
  86. smpl_grp_t *smpl_grp;
  87. int nsmpl_grp;
  88. // ccall only
  89. double indel_frac, min_perm_p, min_lrt;
  90. double prior_type, pref;
  91. int ngrp1_samples, n_perm;
  92. char *prior_file;
  93. ccall_t *cdat;
  94. // shared
  95. bcf_srs_t *srs; // BCF synced readers holding target alleles for CALL_CONSTR_ALLELES
  96. bcf1_t *rec;
  97. bcf_hdr_t *hdr;
  98. uint32_t flag; // One or more of the CALL_* flags defined above
  99. uint8_t *ploidy, all_diploid, unseen;
  100. double pl2p[256]; // PL to 10^(-PL/10) table
  101. int32_t *PLs; // VCF PL likelihoods (rw)
  102. int nPLs, mPLs, nac;
  103. int32_t *gts, *ac; // GTs and AC (w)
  104. double *pdg; // PLs converted to P(D|G)
  105. float *anno16; int n16; // see anno[16] in bam2bcf.h
  106. double theta; // prior
  107. }
  108. call_t;
  109. void error(const char *format, ...);
  110. /*
  111. * call() - return -1 value on critical error; -2 to skip the site; or the number of non-reference
  112. * alleles on success.
  113. */
  114. int mcall(call_t *call, bcf1_t *rec); // multiallic and rare-variant calling model
  115. int ccall(call_t *call, bcf1_t *rec); // the default consensus calling model
  116. int qcall(call_t *call, bcf1_t *rec); // QCall output
  117. void mcall_init(call_t *call);
  118. void ccall_init(call_t *call);
  119. void qcall_init(call_t *call);
  120. void mcall_destroy(call_t *call);
  121. void ccall_destroy(call_t *call);
  122. void qcall_destroy(call_t *call);
  123. void call_init_pl2p(call_t *call);
  124. uint32_t *call_trio_prep(int is_x, int is_son);
  125. void init_allele_trimming_maps(call_t *call, int nals_ori, int als_out);
  126. void mcall_trim_and_update_numberR(call_t *call, bcf1_t *rec, int nals_ori, int nals_new);
  127. #endif