diff --git a/noisyplotype.py b/noisyplotype.py index 055715f8b3a992f3519bdc6a4b51e64a6aa36348..dd6320d3f66ac097ae7c50fa85c75145ca46e0b9 100644 --- a/noisyplotype.py +++ b/noisyplotype.py @@ -209,7 +209,7 @@ def generate_individuals (nb_individuals, chromosome_size, marker_density, mean_ #markers_positions = sorted(random.sample(range(chromosome_size),size)) markers_positions = np.linspace(1, chromosome_size, markers_nb, dtype="int") - with open('Breakpoints.csv', 'w') as breakpointFile: + with open('Breakpoints_3x.csv', 'w') as breakpointFile: breakpointFile.write(",".join(["sample","average_bkp_position", "bkp_start_position", "bkp_stop_position", "transitionType"]) + "\n"); for i in range(nb_individuals): print("individual " + str(i+1)) @@ -241,8 +241,8 @@ chromosome_size = 44000000 cMsize = 180 ## Size of genetic map conversion_factor = chromosome_size/cMsize ## Corresponds to a bpPercM conversion ! Needs to be fixed... Does not produce the correct division! marker_density = 0.0055 -mean_depth = 3 -max_depth = 6 #TODO (better estimate of max?) +mean_depth = 1.5 +max_depth = 3 #TODO (better estimate of max?) errA = 0.005 errB = 0.005 @@ -276,7 +276,7 @@ header = [ ## Write the VCF files (with and without errors) -with open('test.vcf', 'w') as file: +with open('test_3x.vcf', 'w') as file: # Write the header of the file file.writelines(header) @@ -288,11 +288,11 @@ with open('test.vcf', 'w') as file: row = ["chr01", str(markers_positions[m]), ".", "A", "T", ".", ".", ".", "GT:DP:AD:RO:QR:AO:QA:GL"] for i in range(0, len(matrix)): row.append(matrix[i][m]) - row.append("0/0:"+ str(mean_depth) + ":" + str(mean_depth) + ",0:.:.:.:.") - row.append("1/1:"+ str(mean_depth) + ":0," + str(mean_depth) + ":.:.:.:.") + row.append("0/0:"+ str(round(mean_depth)) + ":" + str(round(mean_depth)) + ",0:.:.:.:.") + row.append("1/1:"+ str(round(mean_depth)) + ":0," + str(round(mean_depth)) + ":.:.:.:.") file.write("\t".join(row) + "\n"); -with open('test_error.vcf', 'w') as file: +with open('test_error_3x.vcf', 'w') as file: # Write the header of the file file.writelines(header) @@ -304,6 +304,6 @@ with open('test_error.vcf', 'w') as file: row = ["chr01", str(markers_positions[m]), ".", "A", "T", ".", ".", ".", "GT:DP:AD:RO:QR:AO:QA:GL"] for i in range(0, len(matrix_error)): row.append(matrix_error[i][m]) - row.append("0/0:"+ str(mean_depth) + ":" + str(mean_depth) + ",0:.:.:.:.") - row.append("1/1:"+ str(mean_depth) + ":0," + str(mean_depth) + ":.:.:.:.") + row.append("0/0:"+ str(round(mean_depth)) + ":" + str(round(mean_depth)) + ",0:.:.:.:.") + row.append("1/1:"+ str(round(mean_depth)) + ":0," + str(round(mean_depth)) + ":.:.:.:.") file.write("\t".join(row) + "\n"); \ No newline at end of file