| Previous changeset 0:2c498d40ecde (2012-04-09) Next changeset 2:41ef7e57c2fa (2012-04-09) |
|
Commit message:
add empty snifffers section to mitigate toolshed bug; new version of dpmix |
|
modified:
datatypes_conf.xml genome_diversity/src/dpmix.c |
| b |
| diff -r 2c498d40ecde -r 0df7a9b89f13 datatypes_conf.xml --- a/datatypes_conf.xml Mon Apr 09 12:03:06 2012 -0400 +++ b/datatypes_conf.xml Mon Apr 09 12:40:57 2012 -0400 |
| b |
| @@ -9,4 +9,5 @@ <datatype extension="wsf" type="galaxy.datatypes.wsf:SnpFile" display_in_upload="true"/> <datatype extension="wpf" type="galaxy.datatypes.wsf:SapFile" display_in_upload="true"/> </registration> + <sniffers/> </datatypes> |
| b |
| diff -r 2c498d40ecde -r 0df7a9b89f13 genome_diversity/src/dpmix.c --- a/genome_diversity/src/dpmix.c Mon Apr 09 12:03:06 2012 -0400 +++ b/genome_diversity/src/dpmix.c Mon Apr 09 12:40:57 2012 -0400 |
| [ |
| b'@@ -12,15 +12,15 @@\n * argv[6] = switch penalty (>= 0)\n * argv[7] = file giving heterochromatic intervals (\'-\' means that no file is\n *\t given)\n-* argv[8] = file name to write additional output\n-* argv[9], argv[10], ..., have the form "13:1:paul", "13:2:fred" or\n+* argv[8] = file name for additional output\n+* argv[9], argv[10], ..., have the form "13:1:Peter", "13:2:Paul" or\n *\t "13:0:Mary", meaning that the 13th and 14th columns (base 1)\n *\t give the allele counts for an individual that is in ancestral\n *\t population 1, ancestral population 2, or is a potentially admixed\n *\t individual, resp.\n \n What it does on Galaxy\n-The user specifies two "ancestral" populations (i.e., sources for chromosomes) and a set of potentially admixed individuals, and chooses between the sequence coverage or the estimated genotypes to measure the similarity of genomic intervals in admixed individuals to the two classes of ancestral chromosomes. The user also picks a "switch penalty", typically between 10 and 100. For each potentially admixed individual, the program divides the genome into three "genotypes": (0) homozygous for the second ancestral population (i.e., both chromosomes from that population), (1) heterozygous, or (2) homozygous for the second ancestral population. Parts of a chromosome that are labeled as "heterochromatic" via argv[7] are given the "non-genotype, 3. Smaller values of the switch penalty (corresponding to more ancient admixture events) generally lead to the reconstruction of more frequent changes between genotypes.\n+The user specifies two "ancestral" populations (i.e., sources for chromosomes) and a set of potentially admixed individuals, and chooses between the sequence coverage or the estimated genotypes to measure the similarity of genomic intervals in admixed individuals to the two classes of ancestral chromosomes. The user also picks a "switch penalty", typically between 10 and 100. For each potentially admixed individual, the program divides the genome into three "genotypes": (0) homozygous for the second ancestral population (i.e., both chromosomes from that population), (1) heterozygous, or (2) homozygous for the second ancestral population. Parts of a reference chromosome that are labeled as "heterochromatic" are given the non-genotype, 3. Smaller values of the switch penalty (corresponding to more ancient admixture events) generally lead to the reconstruction of more frequent changes between genotypes.\n */\n \n #include "lib.h"\n@@ -29,20 +29,23 @@\n // maximum length of a line from the table\n #define MOST 5000\n \n+// we create a linked list of "events" on a chromosome -- mostly SNPs, but\n+// also ends of hetorochomatic intervals\n struct snp {\n \tdouble F1, F2;\t// reference allele frequencies in the two populations\n-\tint pos, het_end, *g; // position and an array of admixed genotypes\n+\tint pos, *g,\t// position and an array of admixed genotypes\n+\t type;\t\t// 0 = SNP, 1 = start of het. interval, 2 = end\n \tstruct snp *prev;\t// we keep the list in order of decreasing pos\n } *last;\n \n+// array of potentially admixed individuals\n struct admixed {\n \tchar *name;\n \tint gcol, ge20, gt02;\n-\tlong long x[4];\n+\tlong long x[4];\t\t// number of reference bp in each state\n } A[MOST];\n \n-// information about the specified individuals, namely column and population\n-// for the ancestral populations and genotype columns for the admixed.\n+// information about "ancestral" individuals, namely column and population\n struct ances {\n \tint col, pop;\n \tchar *name;\n@@ -55,8 +58,7 @@\n } H[MOST];\n \n // global variables\n-double *S[4];\t// S[s][i] = best score ending in state s at SNP i\n-int *B[4],\t// backpointer to state at the previous SNP\n+int *B[4],\t// backpointer to state at the previous SNP (or event)\n *P;\t\t// chromosome position\n int nH, nI, nG, genotypes, nsnp, debug, chr_col, logs;\n char this_chr[100];\n@@ -94,7 +96,6 @@\n \t\t\tp = f1*(1.0-f2) + (1.0-f1)*f2;\n \t}\n \t\n-//fprintf(stderr, "%f %f %d %d =>'..b'ypes in the potential admixed individuals\n \t\tfor (i = A1 = B1 = A2 = B2 = 0; i < nI; ++i) {\n \t\t\tn = C[i].col;\n@@ -279,7 +313,7 @@\n \t\tnew->pos = X[chr_col+1];\n \t\tnew->F1 = F1 = (double)A1/(double)(A1+B1);\n \t\tnew->F2 = F2 = (double)A2/(double)(A2+B2);\n-\t\tnew->het_end = 0;\n+\t\tnew->type = 0;\n \t\tnew->g = ckalloc(nG*sizeof(int));\n \t\tfor (i = 0; i < nG; ++i) {\n \t\t\tg = new->g[i] = X[A[i].gcol];\n@@ -307,7 +341,7 @@\n /*\n printf("nsnp = %d\\n", nsnp);\n for (i = nsnp-1, new = last; i >= 0 && new != NULL; --i, new = new->prev) {\n-printf("%d %d ", new->pos, new->het_end);\n+printf("%d %d ", new->pos, new->type);\n printf("%g %g ", new->F1, new->F2);\n for (a = 0; a < nG; ++a)\n printf("%d", new->g[a]);\n@@ -318,14 +352,12 @@\n */\n \n \t// allocate arrays for the DP analysis\n-\tP = ckalloc(nsnp*sizeof(int));\n+\tP = ckalloc(nsnp*sizeof(int));\t// position of each event\n \tfor (i = nsnp-1, new = last; i >= 0 && new != NULL;\n \t --i, new = new->prev)\n \t\tP[i] = new->pos;\n \n-\tfor (i = 0; i < 4; ++i) {\n-\t\tS[i] = ckalloc((nsnp+1)*sizeof(double));\n-\t\tS[i][nsnp] = 0.0;\n+\tfor (i = 0; i < 4; ++i) {\t// space for back-pointers\n \t\tB[i] = ckalloc((nsnp+1)*sizeof(int));\n \t\tB[i][nsnp] = 0;\n \t}\n@@ -342,20 +374,18 @@\n \t\tfree(new);\n \t}\n \tfree(P);\n-\tfor (i = 0; i < 4; ++i) {\n-\t\tfree(S[i]);\n+\tfor (i = 0; i < 4; ++i)\n \t\tfree(B[i]);\n-\t}\n }\n \n int main(int argc, char **argv) {\n \tint n, i, j, k, saw[3];\n-\tlong long het_len;\n+\tlong long het_len, ref_len;\n \tfloat N;\n \tchar nam[100], *chr;\n \n-\tif (argc < 8)\n-\t\tfatal("args: table chr-col chr data-source switch heterochrom outfile n:1:name1 m:2:name2 ...");\n+\tif (argc < 9)\n+\t\tfatal("args: table chr-col chr data-source logs switch heterochrom outfile n:1:name1 m:2:name2 ...");\n \tif (same_string(argv[argc-1], "debug")) {\n \t\tdebug = 1;\n \t\t--argc;\n@@ -365,19 +395,18 @@\n \tchr_col = atoi(argv[2]);\n \tchr = argv[3];\n \tgenotypes = atoi(argv[4]);\n+\n \tlogs = atoi(argv[5]);\n \tif (logs)\n \t\tfatal("logarithms of probabilities -- under development");\n+\t//if (logs) switch_penalty = log(switch_penalty);\n+\n \tswitch_penalty = atof(argv[6]);\n-/*\n-\tif (logs)\n-\t\tswitch_penalty = log(switch_penalty);\n-*/\n \tif (switch_penalty < 0.0)\n \t\tfatal("negative switch penalty");\n \tout = ckopen(argv[8], "w");\n \n-\thet_len = 0;\n+\thet_len = ref_len = 0;\n \tif (!same_string(argv[7], "-")) {\n \t\tfp = ckopen(argv[7], "r");\n \t\twhile (fgets(buf, MOST, fp)) {\n@@ -388,19 +417,18 @@\n \t\t\tH[nH].chr = copy_string(nam);\n \t\t\tH[nH].b = i;\n \t\t\tH[nH].e = j;\n+\t\t\t// assumes last event per chrom. is a het. interval\n+\t\t\tif (nH > 0 && !same_string(nam, H[nH-1].chr))\n+\t\t\t\tref_len += j;\n \t\t\thet_len += (j - i);\n \t\t\t++nH;\n \t\t}\n \t\tfclose(fp);\n \t}\n-/*\n-for (i = 0; i < nH; ++i)\n-printf("%s %d %d\\n", H[i].chr, H[i].b, H[i].e);\n-exit(0);\n-*/\n+\tref_len += H[nH-1].e;\n \n+\t// populations must be disjoint\n \tsaw[1] = saw[2] = 0;\n-\t// populations must be disjoint\n \tfor (i = 9; i < argc; ++i) {\n \t\tif (sscanf(argv[i], "%d:%d:%s", &j, &k, nam) != 3)\n \t\t\tfatalf("not like 13:2:fred : %s", argv[i]);\n@@ -413,12 +441,12 @@\n \t\t\t;\n \t\tif (n < nI)\n \t\t\tfatal("populations are not disjoint");\n-\t\tif (k == 0) {\t// admixed individuals\n+\t\tif (k == 0) {\t// admixed individual\n \t\t\tif (nG >= MOST)\n \t\t\t\tfatal("Too many admixed individuals");\n \t\t\tA[nG].name = copy_string(nam);\n \t\t\tA[nG++].gcol = j+2;\n-\t\t} else {\t// ancestral populations\n+\t\t} else {\t// in an ancestral population\n \t\t\tif (nI >= MOST)\n \t\t\t\tfatal("Too many ancestral individuals");\n \t\t\tC[nI].col = j;\n@@ -426,7 +454,6 @@\n \t\t\tC[nI++].name = copy_string(nam);\n \t\t}\n \t}\n-\t\t\n \tif (saw[0] == 0)\n \t\tfatal("no admixed individual is specified");\n \tif (saw[1] == 0)\n@@ -466,6 +493,12 @@\n \t\t A[i].name, A[i].gt02);\n \t}\n \t// write fractions in each state to the output text file\n+\n+\tif (ref_len)\n+\t\tfprintf(out,\n+\t\t "%lld of %lld reference bp (%1.1f%%) are heterochromatin\\n\\n",\n+\t\t het_len, ref_len, 100.0*(float)het_len/(float)ref_len);\n+\n \tfor (i = 0; i < nG; ++i) {\n \t\tN = (float)(A[i].x[0] + A[i].x[1] + A[i].x[2])/100.0;\n \t\tfprintf(out, "%s: 0 = %1.1f%%, 1 = %1.1f%%, 2 = %1.1f%%\\n",\n' |