|
2 | 2 | import sys, re, json |
3 | 3 | from collections import OrderedDict |
4 | 4 |
|
| 5 | +jsonOut=OrderedDict() |
5 | 6 | data=OrderedDict() |
6 | 7 |
|
| 8 | +## Function to convert a set of elements into floating point numbers, when possible, else leave them be. |
| 9 | +def make_float(x): |
| 10 | + print (x) |
| 11 | + output=[None for i in range(len(x))] |
| 12 | + ## If value for an estimate/error is -nan, replace with "NA". JSON does not accept NaN as a valid field. |
| 13 | + for i in range(len(x)): |
| 14 | + if x[i] == "-nan": |
| 15 | + output[i]="N/A" |
| 16 | + continue |
| 17 | + try: |
| 18 | + output[i]=float(x[i]) |
| 19 | + except: |
| 20 | + output[i]=x[i] |
| 21 | + |
| 22 | + return(tuple(output)) |
| 23 | + |
| 24 | + |
7 | 25 | Input_files=sys.argv[1:] |
8 | 26 |
|
9 | 27 | output = open("nuclear_contamination.txt", 'w') |
|
16 | 34 | ml2, err_ml2="N/A","N/A" |
17 | 35 | with open(fn, 'r') as f: |
18 | 36 | Estimates={} |
19 | | - Ind=re.sub('\.X.contamination.out$', '', fn) |
| 37 | + Ind=re.sub('\.X.contamination.out$', '', fn).split("/")[-1] |
20 | 38 | for line in f: |
21 | 39 | fields=line.strip().split() |
22 | 40 | if line.strip()[0:19] == "We have nSNP sites:": |
|
26 | 44 | err_mom1=fields[4].split(":")[1] |
27 | 45 | ml1=fields[5].split(":")[1] |
28 | 46 | err_ml1=fields[6].split(":")[1] |
29 | | - ## Sometimes angsd fails to run method 2, and the error is printed directly after the SE for ML. When that happens, exclude the first word in the error from the output. (Method 2 data will be shown as NA) |
| 47 | + ## Sometimes angsd fails to run method 2, and the error is printed directly after the SE for ML. When that happens, exclude the first word in the error from the output. (Method 2 jsonOut will be shown as NA) |
30 | 48 | if err_ml1.endswith("contamination"): |
31 | 49 | err_ml1 = err_ml1[:-13] |
32 | 50 | elif line.strip()[0:7] == "Method2" and line.strip()[9:16] == 'new_llh': |
33 | 51 | mom2=fields[3].split(":")[1] |
34 | 52 | err_mom2=fields[4].split(":")[1] |
35 | 53 | ml2=fields[5].split(":")[1] |
36 | 54 | err_ml2=fields[6].split(":")[1] |
37 | | - data[Ind]={ "Number_of_SNPs" : nSNPs, "Method1_MOM_estimate" : mom1, "Method1_MOM_SE" : err_mom1, "Method1_ML_estimate" : ml1, "Method1_ML_SE" : err_ml1, "Method2_MOM_estimate" : mom2, "Method2_MOM_SE" : err_mom2, "Method2_ML_estimate" : ml2, "Method2_ML_SE" : err_ml2 } |
| 55 | + ## Convert estimates and errors to floating point numbers |
| 56 | + (ml1, err_ml1, mom1, err_mom1, ml2, err_ml2, mom2, err_mom2) = make_float((ml1, err_ml1, mom1, err_mom1, ml2, err_ml2, mom2, err_mom2)) |
| 57 | + data[Ind]={ "Num_SNPs" : int(nSNPs), "Method1_MOM_estimate" : mom1, "Method1_MOM_SE" : err_mom1, "Method1_ML_estimate" : ml1, "Method1_ML_SE" : err_ml1, "Method2_MOM_estimate" : mom2, "Method2_MOM_SE" : err_mom2, "Method2_ML_estimate" : ml2, "Method2_ML_SE" : err_ml2 } |
38 | 58 | print (Ind, nSNPs, mom1, err_mom1, ml1, err_ml1, mom2, err_mom2, ml2, err_ml2, sep="\t", file=output) |
39 | 59 |
|
40 | | -with open('nuclear_contamination.json', 'w') as outfile: |
41 | | - json.dump(data, outfile) |
| 60 | + |
| 61 | +jsonOut = {"plot_type": "generalstats", "id": "nuclear_contamination", |
| 62 | + "pconfig": { |
| 63 | + "Num_SNPs" : {"title" : "Number of SNPs"}, |
| 64 | + "Method1_MOM_estimate" : {"title": "Contamination Estimate (Method1_MOM)"}, |
| 65 | + "Method1_MOM_SE" : {"title": "Estimate Error (Method1_MOM)"}, |
| 66 | + "Method1_ML_estimate" : {"title": "Contamination Estimate (Method1_ML)"}, |
| 67 | + "Method1_ML_SE" : {"title": "Estimate Error (Method1_ML)"}, |
| 68 | + "Method2_MOM_estimate" : {"title": "Contamination Estimate (Method2_MOM)"}, |
| 69 | + "Method2_MOM_SE" : {"title": "Estimate Error (Method2_MOM)"}, |
| 70 | + "Method2_ML_estimate" : {"title": "Contamination Estimate (Method2_ML)"}, |
| 71 | + "Method2_ML_SE" : {"title": "Estimate Error (Method2_ML)"} |
| 72 | + }, |
| 73 | + "data" : data |
| 74 | +} |
| 75 | +with open('nuclear_contamination_mqc.json', 'w') as outfile: |
| 76 | + json.dump(jsonOut, outfile) |
0 commit comments