adaptive-pulldown/exclude.c at master · DReichLab/adaptive-pulldown · GitHub

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
#include <exclude.h>
#include <math.h>
#include <mcio.h>

#define MAXRGN 1000

void
excluderegions (char *xregionname, SNP ** snps, int nsnps,
		char *deletesnpoutname)
{
  FILE *fp;

  int chr[MAXRGN];
  int lo[MAXRGN];
  int hi[MAXRGN];

  char line[MAXSTR];
  char *spt[MAXFF];
  int nsplit, nrgn, i, j;

  if ((fp = fopen (xregionname, "r")) == NULL) {
    printf ("excluderegions: can't open file %s\n", xregionname);
    return;
  }

  for (i = 0; i < MAXRGN; i++) {

    if (fgets (line, MAXSTR, fp) == NULL)
      break;

    nsplit = splitup (line, spt, MAXFF);
    if (nsplit != 3)
      continue;

    chr[i] = atoi (spt[0]);
    lo[i] = atoi (spt[1]);
    hi[i] = atoi (spt[2]);

  }
  fclose (fp);
  nrgn = i;

  for (i = 0; i < nsnps; i++) {
    SNP *cupt = snps[i];
    for (j = 0; j < nrgn; j++) {
      if (cupt->chrom == chr[j] && cupt->physpos >= lo[j]
	  && cupt->physpos <= hi[j]) {
	cupt->ignore = YES;
	if (deletesnpoutname != NULL) {
	  logdeletedsnp (cupt->ID, "xregion", deletesnpoutname);
	}
      }
    }
  }

  return;

}

void
hwfilter (SNP ** snps, int nsnps, int nindiv, double nhwfilter,
	  char *deletesnpoutname)
{

  int i, k;

  for (i = 0; i < nsnps; i++) {
    int num = 0, den = 0, het = 0, n0 = 0, n1 = 0, n2 = 0, nsamples;
    double p, Q, stdv;
    SNP *cupt = snps[i];

    for (k = 0; k < nindiv; k++) {
      int g = getgtypes (cupt, k);
      if (g >= 0) {
	num += g;
	den += 2;
      }
      if (g == 1) {
	het++;
	n1++;
      }
      else if (g == 0) {
	n0++;
      }
      else if (g == 2) {
	n2++;
      }
    }

    if ((nsamples = den / 2) == 0)
      continue;
    p = (double) num / den;
    Q = 2 * p * (1 - p);
    stdv = sqrt (Q * (1 - Q) / nsamples);
    if (fabs ((double) het / nsamples - Q) > nhwfilter * stdv) {
      printf ("SNP %s removed by Hardy-Weinberg filter\n", cupt->ID);
      cupt->ignore = YES;
      if (deletesnpoutname != NULL) {
	logdeletedsnp (cupt->ID, "hwfilt", deletesnpoutname);
      }
    }
  }


}