-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathdocumentation_Jensen_response.doc
501 lines (321 loc) · 23.9 KB
/
documentation_Jensen_response.doc
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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
#####################################################################################
# Measure H12 in 10KB windows under neutrality for ALL the different demographic models
####################################################################################
qsub ~/Jensen_response/scripts/qsub_constNe10e6_theta0
qsub ~/Jensen_response/scripts/qsub_constNe2.7e6_theta0
qsub ~/Jensen_response/scripts/qsub_dadi1_theta0
qsub ~/Jensen_response/scripts/qsub_dadi2_theta0
qsub ~/Jensen_response/scripts/qsub_admixture_bottleneck_mode_theta0
qsub ~/Jensen_response/scripts/qsub_admixture_mode_theta0
qsub ~/Jensen_response/scripts/qsub_admixture_posterior_95CI_theta0
qsub ~/Jensen_response/scripts/qsub_admixture_posterior_95CI_Harris_theta0
qsub ~/Jensen_response/scripts/qsub_admixture_posterior_95CI_Harris_theta0_rho5e-7
#qsub ~/Jensen_response/scripts/qsub_admixture_posterior_95CI_Harris_Nac_theta0
qsub ~/Jensen_response/scripts/qsub_admixture_mode_Garud_2015_theta0
qsub ~/Jensen_response/scripts/qsub_admixture_bot_mode_Garud_2015_theta0
qsub ~/Jensen_response/scripts/qsub_Arguello
# plot S, Pi, TajD, H12, and H2/H1 for these simulations vs the real data
python ~/Jensen_response/scripts/plot_original_models.py
# plot long-range LD for these scenarios too
#qsub ~/Jensen_response/scripts/qsub_constNe10e6_theta0
#qsub ~/Jensen_response/scripts/qsub_constNe2.2e6_theta0
#qsub ~/Jensen_response/scripts/qsub_dadi2_theta0_recoded_MS
#qsub ~/Jensen_response/scripts/qsub_admixture_mode_theta0
#qsub ~/Jensen_response/scripts/qsub_admixture_mode_rho5e-7_theta0
#qsub ~/Jensen_response/scripts/qsub_admixture_bottleneck_mode_theta0 # corrected this one
#qsub_admixture_mode_corrected_theta0
#qsub ~/Jensen_response/scripts/qsub_admixture_mode_theta0_MS
#qsub_admixture_mode_hardcoded_theta0
#qsub ~/Jensen_response/scripts/qsub_admixture_posterior_95CI_theta0 #corrected
#qsub ~/Jensen_response/scripts/qsub_admixture_posterior_95CI_theta0_MS
qsub ~/Jensen_response/scripts/qsub_admixture_uniform_95CI_theta0
qsub ~/Jensen_response/scripts/qsub_admixture_posterior_outside_95CI_theta0
### Run for X chr to check against Duchen's posterior ###
qsub ~/Jensen_response/scripts/qsub_admixture_95CI_DuchenRhoTheta_exactLength_theta0_Xchr
qsub ~/Jensen_response/scripts/qsub_admixture_fullPosterior_DuchenRhoTheta_exactLength_theta0_Xchr
##########################################
# H12 in 10kb windows for hard and soft #
# Admixture uniform, Admixture posterior #
##########################################
qsub ~/Jensen_response/scripts/qsub_admixture_posterior_95CI_theta0.01 #corrected
qsub ~/Jensen_response/scripts/qsub_admixture_posterior_95CI_theta10 #corrected
qsub ~/Jensen_response/scripts/qsub_admixture_uniform_95CI_theta0.01
qsub ~/Jensen_response/scripts/qsub_admixture_uniform_95CI_theta10
##############################################################
# H12 in 400kb windows for neutrality, hard, and soft sweeps #
# Admixture uniform, Admixture posterior #
##############################################################
qsub ~/Jensen_response/scripts/qsub_admixture_posterior_outside_95CI_theta0_400SNPs
qsub ~/Jensen_response/scripts/qsub_admixture_posterior_95CI_theta0_400SNPs #corrected
qsub ~/Jensen_response/scripts/qsub_admixture_posterior_95CI_theta0.01_400SNPs
qsub ~/Jensen_response/scripts/qsub_admixture_posterior_95CI_theta10_400SNPs
qsub ~/Jensen_response/scripts/qsub_admixture_uniform_95CI_theta0_400SNPs
qsub ~/Jensen_response/scripts/qsub_admixture_uniform_95CI_theta0.01_400SNPs
qsub ~/Jensen_response/scripts/qsub_admixture_uniform_95CI_theta10_400SNPs
qsub ~/Jensen_response/scripts/qsub_admixture_bottleneck_mode_theta0_400SNPs # corrected this one
########################################################################
# Check if there are 2 entries for Pi output, otherwise discard line: #
########################################################################
for file in admixture_bottleneck_mode_neutrality.txt_Pi_S.txt \
Admixture_mode_rho_5e-9_theta_0_selection_False_Pi_S.txt \
Admixture_posteriorDistn_95CI_rho_5e-9_theta_0.01_selection_True_Pi_S_bps.txt \ Admixture_posteriorDistn_95CI_rho_5e-9_theta_0_selection_False_Pi_S_bps.txt \
Admixture_posteriorDistn_95CI_rho_5e-9_theta_10_selection_True_Pi_S_bps.txt \
Admixture_posteriorDistn_outside_95CI_rho_5e-9_theta_0_selection_False_Pi_S_bps.txt \
Admixture_uniformPrior_rho_5e-9_theta_0.01_selection_True_Pi_S_bps.txt \
Admixture_uniformPrior_rho_5e-9_theta_0_selection_False_Pi_S_bps.txt \
Admixture_uniformPrior_rho_5e-9_theta_10_selection_True_Pi_S_bps.txt \
constNe106_neutrality_Pi_S.txt \
constNe2.2e6_neutrality_Pi_S.txt \
dadi1_neutrality_Pi_S.txt \
dadi2_neutrality_Pi_S.txt \
Admixture_mode_rho_5e-9_theta_0_selection_False_Pi_S_noSingletons_DuchenRhoTheta.txt \
Admixture_mode_rho_5e-9_theta_0_selection_False_Pi_S_Singletons_DuchenRhoTheta.txt ; do
for file in Admixture_mode_rho_5e-9_theta_0_selection_False_Pi_S_noSingletons.txt\
admixture_bottleneck_mode_neutrality.txt_Pi_S_noSingletons.txt \
Admixture_mode_rho_5e-9_theta_0_selection_False_Pi_S_noSingletons_DuchenRhoTheta.txt \
Admixture_mode_rho_5e-9_theta_0_selection_False_Pi_S_noSingletons.txt \
Admixture_posteriorDistn_95CI_rho_5e-9_theta_0_selection_False_Pi_S_noSingletons.txt \
Admixture_posteriorDistn_outside_95CI_rho_5e-9_theta_0_selection_False_Pi_S_noSingletons.txt \
Admixture_uniformPrior_rho_5e-9_theta_0_selection_False_Pi_S_noSingletons.txt; do
for file in constNe2.2e6_neutrality_Pi_S.txt dadi1_neutrality_Pi_S.txt dadi2_neutrality_Pi_S.txt; do
for file in Admixture_mode_corrected_rho_5e-9_theta_0_selection_False_Pi_S_noSingletons.txt Admixture_mode_corrected_rho_5e-9_theta_0_selection_False_Pi_S.txt; do
for file in admixture_bottleneck_mode_corrected_neutrality.txt_Pi_S_noSingletons.txt \
admixture_bottleneck_mode_corrected_neutrality.txt_Pi_S.txt \
Admixture_posteriorDistn_95CI_corrected_rho_5e-9_theta_0_selection_False_Pi_S_bps.txt \
Admixture_posteriorDistn_95CI_corrected_rho_5e-9_theta_0_selection_False_Pi_S_noSingletons.txt; do
echo $file
python ~/Jensen_response/scripts/count_fields_Pi_S.py ${file} tmpOut.txt
mv tmpOut.txt ${file}
done
#################################
# check if there are 3 entries #
#################################
for file in Admixture_mode_rho_5e-9_theta_0_selection_False_DuchenRhoTheta_exactLength_Pi_S_TajD_noSingletons_DuchenRhoTheta.txt Admixture_mode_rho_5e-9_theta_0_selection_False_DuchenRhoTheta_exactLength_Pi_S_TajD_Singletons_DuchenRhoTheta.txt Admixture_95CI_rho_5e-9_theta_0_selection_False_DuchenRhoTheta_exactLength_Pi_S_TajD_noSingletons_DuchenRhoTheta.txt Admixture_95CI_rho_5e-9_theta_0_selection_False_DuchenRhoTheta_exactLength_Pi_S_TajD_Singletons_DuchenRhoTheta.txt ; do
for file in Admixture_mode_corrected_hardcoded_rho_5e-9_theta_0_selection_False_Pi_S_TajD.txt Admixture_mode_corrected_hardcoded_rho_5e-9_theta_0_selection_False_Pi_S_TajD_noSingletons.txt ; do
for file in Admixture_95CI_rho_5e-9_theta_0_selection_False_DuchenRhoTheta_exactLength_Xchr_Pi_S_TajD_noSingletons_DuchenRhoTheta.txt \
Admixture_95CI_rho_5e-9_theta_0_selection_False_DuchenRhoTheta_exactLength_Xchr_Pi_S_TajD_Singletons_DuchenRhoTheta.txt \
Admixture_fullPosterior_rho_5e-9_theta_0_selection_False_DuchenRhoTheta_exactLength_Xchr_Pi_S_TajD_noSingletons_DuchenRhoTheta.txt \
Admixture_fullPosterior_rho_5e-9_theta_0_selection_False_DuchenRhoTheta_exactLength_Xchr_Pi_S_TajD_Singletons_DuchenRhoTheta.txt; do
for file in Admixture_posteriorDistn_95CI_corrected_rho_5e-9_theta_0_selection_False_Pi_S_TajD_bps.txt Admixture_posteriorDistn_95CI_corrected_rho_5e-9_theta_0_selection_False_Pi_S_TajD_noSingletons.txt Admixture_mode_corrected_rho_5e-9_theta_0_selection_False_Pi_S_TajD_noSingletons.txt Admixture_mode_corrected_rho_5e-9_theta_0_selection_False_Pi_S_TajD.txt; do
for file in admixture_bottleneck_mode_corrected_neutrality.txt_Pi_S_TajD_noSingletons.txt admixture_bottleneck_mode_corrected_neutrality.txt_Pi_S_TajD.txt; do
python ~/Jensen_response/scripts/count_fields_Pi_S_TajD.py ${file} tmpOut.txt
mv tmpOut.txt ${file}
done
##########################################################################
# Check if there are 17 entries for H12 output, otherwise discard line: #
##########################################################################
for file in admixture_bottleneck_mode_neutrality.txt_bps.txt Admixture_mode_rho_5e-9_theta_0_selection_False_bps.txt Admixture_posteriorDistn_95CI_rho_5e-9_theta_0.01_selection_True_bps.txt Admixture_posteriorDistn_95CI_rho_5e-9_theta_0.01_selection_True_snps.txt Admixture_posteriorDistn_95CI_rho_5e-9_theta_0_selection_False_bps.txt Admixture_posteriorDistn_95CI_rho_5e-9_theta_0_selection_False_snps.txt Admixture_posteriorDistn_95CI_rho_5e-9_theta_10_selection_True_bps.txt Admixture_posteriorDistn_95CI_rho_5e-9_theta_10_selection_True_snps.txt Admixture_uniformPrior_rho_5e-9_theta_0.01_selection_True_bps.txt Admixture_uniformPrior_rho_5e-9_theta_0.01_selection_True_snps.txt Admixture_uniformPrior_rho_5e-9_theta_0_selection_False_bps.txt Admixture_uniformPrior_rho_5e-9_theta_0_selection_False_snps.txt Admixture_uniformPrior_rho_5e-9_theta_10_selection_True_bps.txt Admixture_uniformPrior_rho_5e-9_theta_10_selection_True_snps.txt constNe106_neutrality_bps.txt
for file in constNe2.2e6_neutrality_bps.txt dadi1_neutrality_bps.txt dadi2_neutrality_bps.txt ; do
for file in admixture_bottleneck_mode_corrected_neutrality.txt_bps.txt \
Admixture_mode_corrected_rho_5e-9_theta_0_selection_False_bps.txt \
Admixture_posteriorDistn_95CI_corrected_rho_5e-9_theta_0.01_selection_True_bps.txt \
Admixture_posteriorDistn_95CI_corrected_rho_5e-9_theta_0.01_selection_True_snps.txt \
Admixture_posteriorDistn_95CI_corrected_rho_5e-9_theta_0_selection_False_bps.txt \
Admixture_posteriorDistn_95CI_corrected_rho_5e-9_theta_0_selection_False_snps.txt \
Admixture_posteriorDistn_95CI_corrected_rho_5e-9_theta_10_selection_True_bps.txt \
Admixture_posteriorDistn_95CI_corrected_rho_5e-9_theta_10_selection_True_snps.txt; do
for file in Admixture_mode_corrected_hardcoded_rho_5e-9_theta_0_selection_False_bps.txt; do
for file in Admixture_posteriorDistn_95CI_corrected_rho_5e-9_theta_0_selection_False_bps.txt Admixture_posteriorDistn_95CI_corrected_rho_5e-9_theta_0_selection_False_snps.txt Admixture_mode_corrected_rho_5e-9_theta_0_selection_False_bps.txt admixture_bottleneck_mode_corrected_neutrality.txt_bps.txt admixture_bottleneck_mode_corrected_neutrality.txt_snps.txt ; do
python ~/Jensen_response/scripts/count_fields.py ${file} tmpOut.txt
mv tmpOut.txt ${file}
done
for file in Admixture_mode_fixedPopSize_NAm1110000_EUr700000_rho_5e-9_theta_0_selection_False_MS_snps.txt Admixture_mode_fixedPopSize_NAm15984500_EUr700000_rho_5e-9_theta_0_selection_False_MS_snps.txt; do
python ~/Jensen_response/scripts/count_fields.py ${file} tmpOut.txt
mv tmpOut.txt ${file}
done
for file in Arguello_neutrality_locusLen100000_snps.txt admixture_posterior_Harris_neutrality_locusLen100000_snps.txt admixture_bot_mode_Garud2015_neutrality_locusLen100000_snps.txt admixture_posterior_neutrality_locusLen100000_snps.txt admixture_bot_mode_neutrality_locusLen100000_snps.txt constNe10e6_neutrality_locusLen100000_snps.txt admixture_mode_Garud2015_neutrality_locusLen100000_snps.txt constNe2.7e6_neutrality_locusLen100000_snps.txt admixture_mode_neutrality_locusLen100000_snps.txt dadi1_neutrality_locusLen100000_snps.txt admixture_posterior_Harris_Nac_neutrality_locusLen100000_snps.txt dadi2_neutrality_locusLen100000_snps.txt; do
python ~/Jensen_response/scripts/count_fields.py ${file} tmpOut.txt
mv tmpOut.txt ${file}
done
#compute Bayes factors
num_sims=16000
for win_type in bps snps; do
for model in posteriorDistn_95CI uniformPrior; do
hard_in=~/Jensen_response/analysis/Admixture_${model}_rho_5e-9_theta_0.01_selection_True_${win_type}.txt
soft_in=~/Jensen_response/analysis/Admixture_${model}_rho_5e-9_theta_10_selection_True_${win_type}.txt
outFile=~/Jensen_response/analysis/Admixture_${model}_rho_5e-9_theta_0.01_theta_10_selection_True_${win_type}_BFs_numSim_${num_sims}.txt
python ~/Jensen_response/scripts/compute_Bayes_factors.py $hard_in $soft_in $outFile $num_sims
done
done
#longer chromosomes
#for win_type in snps bps; do
for win_type in bps; do
for num_sims in 100000 16000; do
hard_in=~/Jensen_response/analysis/Admixture_rho_5e-9_theta_0.01_selection_True_${win_type}.txt
soft_in=~/Jensen_response/analysis/Admixture_rho_5e-9_theta_10_selection_True_${win_type}.txt
outFile=~/Jensen_response/analysis/Admixture_rho_5e-9_theta_0.01_theta_10_selection_True_${win_type}_BFs_numSim_${num_sims}.txt
python ~/Jensen_response/scripts/compute_Bayes_factors.py $hard_in $soft_in $outFile $num_sims
done
done
for win_type in bps; do
for num_sims in 16000; do
hard_in=~/Jensen_response/analysis/Admixture_rho_5e-9_theta_0.01_selection_True_jointDistn_${win_type}.txt
soft_in=~/Jensen_response/analysis/Admixture_rho_5e-9_theta_10_selection_True_jointDistn_${win_type}.txt
outFile=~/Jensen_response/analysis/Admixture_rho_5e-9_theta_0.01_theta_10_selection_True_jointDistn_${win_type}_BFs_numSim_${num_sims}.txt
python ~/Jensen_response/scripts/compute_Bayes_factors.py $hard_in $soft_in $outFile $num_sims
done
done
for win_type in bps; do
for num_sims in 16000; do
hard_in=~/Jensen_response/analysis/Admixture_rho_5e-9_theta_0.01_selection_True_mode_${win_type}.txt
soft_in=~/Jensen_response/analysis/Admixture_rho_5e-9_theta_10_selection_True_mode_${win_type}.txt
outFile=~/Jensen_response/analysis/Admixture_rho_5e-9_theta_0.01_theta_10_selection_True_mode_${win_type}_BFs_numSim_${num_sims}.txt
python ~/Jensen_response/scripts/compute_Bayes_factors.py $hard_in $soft_in $outFile $num_sims
done
done
#shorter chromosomes
for win_type in bps; do
for num_sims in 100000 16000; do
hard_in=~/Jensen_response/analysis/Admixture_rho_5e-9_theta_0.01_selection_True_${win_type}_10000bps.txt
soft_in=~/Jensen_response/analysis/Admixture_rho_5e-9_theta_10_selection_True_${win_type}_10000bps.txt
outFile=~/Jensen_response/analysis/Admixture_rho_5e-9_theta_0.01_theta_10_selection_True_${win_type}_BFs_10000bps_numSim_${num_sims}.txt
python ~/Jensen_response/scripts/compute_Bayes_factors.py $hard_in $soft_in $outFile ${num_sims}
done
done
# compute the top 50 peaks H12 and H2/H1 values in fixed 10kb windows
while read line; do
chr=`echo $line | cut -f1 -d' ' | cut -f2 -d'r'`
coord=`echo $line | cut -f2 -d' '`
echo $chr
echo $coord
python ~/Jensen_response/scripts/H12_H2H1_fixedWinOption.py /pollard/home/ngarud/westway-home/scalegen/NanditaWork/DGRP_processed/Variants_Sparse_${chr}.sample_swap_fixed_noHeader_noHaps_Ns_noInvariants.txt 145 -o tmp.txt -w 10000 -s $coord -f True
cat tmp.txt >> ~/Jensen_response/analysis/top50peaks_fixed_bp.txt
done < ~/Jensen_response/analysis/peaks_coordinates.txt
Pi_MS.py
H12=`cat Admixture_rho_5e-9_theta_0_selection_False_posteriorDistn_H12_command.txt | cut -f1`
Nac=`cat Admixture_rho_5e-9_theta_0_selection_False_posteriorDistn_H12_command.txt | cut -f2 | cut -f5 -d' '`
Nec=`cat Admixture_rho_5e-9_theta_0_selection_False_posteriorDistn_H12_command.txt | cut -f2 | cut -f15 -d' '`
Nnc=`cat Admixture_rho_5e-9_theta_0_selection_False_posteriorDistn_H12_command.txt | cut -f2 | cut -f21 -d' '`
while read line; do
H12=`echo $line | cut -f1 -d' '`
Nac=`echo $line | cut -f6 -d' '`
Nec=`echo $line | cut -f16 -d' '`
Nnc=`echo $line | cut -f22 -d' '`
echo -e "$H12\t$Nac\t$Nec\t$Nnc" >> H12_popsize_95.txt
done < Admixture_rho_5e-9_theta_0_selection_False_posteriorDistn_95CI_H12_command.txt
#Admixture_rho_5e-9_theta_0_selection_False_posteriorDistn_H12_command.txt
##############################################################
# Run admixture with exact rho, mu, and for the x chr to generate summary statistics for Duchen et al.
###########################################
###############
# experiments #
###############
#############
# migration #
#############
qsub_admixture_mode_theta0_migration_400SNPs
plot_admixture_migration.py
for m_NA in 0 0.1 0.25 0.5 0.75; do
file=~/Jensen_response/analysis/Admixture_mode_migration_mNA${m_NA}_mNE${m_NA}_rho_5e-9_theta_0_selection_False_LD_0.05_0.95.txt
python ~/Jensen_response/scripts/count_fields_any.py $file tmp.txt 8
# average the simulation results for LD:
cat tmp.txt | head -10000000 | sort -k4,4 -g > ~/Jensen_response/analysis/Admixture_mode_migration_mNA${m_NA}_mNE${m_NA}_rho_5e-9_theta_0_selection_False_LD_0.05_0.95_sorted.txt
python ~/Jensen_response/scripts/averageLD_realData2.py ~/Jensen_response/analysis/Admixture_mode_migration_mNA${m_NA}_mNE${m_NA}_rho_5e-9_theta_0_selection_False_LD_0.05_0.95_sorted.txt ~/Jensen_response/analysis/Admixture_mode_migration_mNA${m_NA}_mNE${m_NA}_rho_5e-9_theta_0_selection_False_LD_0.05_0.95_averaged.txt
done
##############################
# vary admixture proportions #
##############################
qsub_admixture_mode_theta0_diffProps_400SNPs
plot_admixture_diffProps.py
for prop in 0 0.1 0.25 0.4 0.5 0.75 0.85 0.9; do
file=~/Jensen_response/analysis/Admixture_mode_migration_diffProps${prop}_rho_5e-9_theta_0_selection_False_MS_LD_0.05_0.95.txt
python ~/Jensen_response/scripts/count_fields_any.py $file tmp.txt 8
# average the simulation results for LD:
cat tmp.txt | head -10000000 | sort -k4,4 -g > ~/Jensen_response/analysis/Admixture_mode_migration_diffProps${prop}_rho_5e-9_theta_0_selection_False_MS_LD_0.05_0.95_sorted.txt
python ~/Jensen_response/scripts/averageLD_realData2.py ~/Jensen_response/analysis/Admixture_mode_migration_diffProps${prop}_rho_5e-9_theta_0_selection_False_MS_LD_0.05_0.95_sorted.txt ~/Jensen_response/analysis/Admixture_mode_migration_diffProps${prop}_rho_5e-9_theta_0_selection_False_MS_LD_0.05_0.95_averaged.txt
done
##########################
# fixed population sizes #
##########################
qsub ~/Jensen_response/scripts/qsub_admixture_mode_theta0_fixedPopSize_400SNPs_MS
plot_admixture_fixedPopSize.py
# also test hard vs soft sweeps for fixed sizes (make BF plots)
# scenarios we want:
NAm=1110000
EUr=700000
qsub ~/Jensen_response/scripts/qsub_admixture_mode_theta10_fixedPopSize_400SNPs_MS
qsub ~/Jensen_response/scripts/qsub_admixture_mode_theta0.01_fixedPopSize_400SNPs_MS
# what do the trajectories look like for adaptive mutations
bash run_admixture_mode_fixedPopSize_400SNPs_MS_trajectory.sh
#plot the trajectories
python plot_trajectories.py
#try conditoining on the ending frequency of a sweep
qsub ~/Jensen_response/scripts/qsub_admixture_mode_theta10_fixedPopSize_400SNPs_MS_conditionEndFreq
qsub ~/Jensen_response/scripts/qsub_admixture_mode_theta0.01_fixedPopSize_400SNPs_MS_conditionEndFreq
#num_sims=20000
#for NAm in 61659 1110000 15984500; do
# for EUr in 500000 2000000; do
num_sims=100000
for NAm in 1110000 15984500; do
for EUr in 700000; do
hard_in=~/Jensen_response/analysis/Admixture_mode_fixedPopSize_NAm${NAm}_EUr${EUr}_rho_5e-9_theta_0.01_selection_True_MS_snps.txt
soft_in=~/Jensen_response/analysis/Admixture_mode_fixedPopSize_NAm${NAm}_EUr${EUr}_rho_5e-9_theta_10_selection_True_MS_snps.txt
outFile=~/Jensen_response/analysis/Admixture_mode_fixedPopSize_NAm${NAm}_EUr${EUr}_rho_5e-9_BFs_numSim_${num_sims}.txt
python ~/Jensen_response/scripts/compute_Bayes_factors.py $hard_in $soft_in $outFile $num_sims
done
done
for NAm in 2500 61659 1110000 15984500 28000000; do
echo 'NAm'
echo $NAm
for EUr in 16982 67608 700000 2000000 9550000; do
echo 'EUr'
echo $EUr
file=~/Jensen_response/analysis/Admixture_mode_fixedPopSize_NAm${NAm}_EUr${EUr}_rho_5e-9_theta_0_selection_False_MS_LD_0.05_0.95.txt
python ~/Jensen_response/scripts/count_fields_any.py $file tmp.txt 8
# average the simulation results for LD:
cat tmp.txt | head -10000000 | sort -k4,4 -g > ~/Jensen_response/analysis/Admixture_mode_fixedPopSize_NAm${NAm}_EUr${EUr}_rho_5e-9_theta_0_selection_False_MS_LD_0.05_0.95_sorted.txt
python ~/Jensen_response/scripts/averageLD_realData2.py ~/Jensen_response/analysis/Admixture_mode_fixedPopSize_NAm${NAm}_EUr${EUr}_rho_5e-9_theta_0_selection_False_MS_LD_0.05_0.95_sorted.txt ~/Jensen_response/analysis/Admixture_mode_fixedPopSize_NAm${NAm}_EUr${EUr}_rho_5e-9_theta_0_selection_False_MS_LD_0.05_0.95_averaged.txt
done
done
#####################
# vary growth rates #
#####################
#qsub_admixture_mode_theta0_varyGrowth_400SNPs_MS
qsub_admixture_mode_theta0_varyGrowth_400SNPs_MS_for_supp
plot_admixture_varyGrowth.py
for NAm_anc in 2500 61659; do
for NAm_cur in 1110000 15984500 28000000; do
for EUr_anc in 16982 67608; do
for EUr_cur in 700000 2000000 9550000; do
file=~/Jensen_response/analysis/Admixture_mode_varyGrowth_NAm_anc${NAm_anc}_NAm_cur${NAm_cur}_EUr_anc${EUr_anc}_EUr_cur${EUr_cur}_rho_5e-9_theta_0_selection_False_MS_LD_0.05_0.95.txt
python ~/Jensen_response/scripts/count_fields_any.py $file tmp.txt 8
# average the simulation results for LD:
cat tmp.txt | head -10000000 | sort -k4,4 -g > ~/Jensen_response/analysis/Admixture_mode_varyGrowth_NAm_anc${NAm_anc}_NAm_cur${NAm_cur}_EUr_anc${EUr_anc}_EUr_cur${EUr_cur}_rho_5e-9_theta_0_selection_False_MS_LD_0.05_0.95_sorted.txt
python ~/Jensen_response/scripts/averageLD_realData2.py ~/Jensen_response/analysis/Admixture_mode_varyGrowth_NAm_anc${NAm_anc}_NAm_cur${NAm_cur}_EUr_anc${EUr_anc}_EUr_cur${EUr_cur}_rho_5e-9_theta_0_selection_False_MS_LD_0.05_0.95_sorted.txt ~/Jensen_response/analysis/Admixture_mode_varyGrowth_NAm_anc${NAm_anc}_NAm_cur${NAm_cur}_EUr_anc${EUr_anc}_EUr_cur${EUr_cur}_rho_5e-9_theta_0_selection_False_MS_LD_0.05_0.95_averaged.txt
done
done
done
done
#####################
# Short intron data #
#####################
# recreate Pi, S, TajD statistics
# Test out S and Pi in Admixture as compared to the short intron data
# convert the shortIntronSNPs_noX.txt file format
#sed '1d' ~/Jensen_response/data/short_introns/shortIntronSNPs_noX.txt > ~/Jensen_response/data/short_introns/tmp.txt
#python ~/Jensen_response/scripts/convertFile.py ~/Jensen_response/data/short_introns/tmp.txt ~/Jensen_response/data/short_introns/shortIntronSNPs_noX_reformatted.txt
#python ~/Jensen_response/scripts/Pi_MS_TajD_perBp.py ~/Jensen_response/data/short_introns/shortIntronSNPs_noX_reformatted.txt ~/Jensen_response/analysis/shortIntronSNPs_S_Pi_TajD_perBp.txt 162 738024
python ~/Jensen_response/scripts/DaDi_compute_S_Pi_TajD.py
#################################
# Get the long range LD working #
#################################
# Check if there are 8 entries
for model in admixture_bot_mode admixture_mode admixture_posterior constNe10e6 constNe2.7e6 dadi1 dadi2 admixture_posterior_Harris admixture_mode_Garud2015 admixture_bot_mode_Garud2015 Arguello; do
echo $model
python ~/Jensen_response/scripts/count_fields_any.py ~/Jensen_response/analysis/${model}_neutrality_locusLen11000_LD_0.05_0.95.txt tmp.txt 8
# average the simulation results for LD:
cat tmp.txt | head -10000000 | sort -k4,4 -g > ~/Jensen_response/analysis/${model}_neutrality_locusLen11000_LD_0.05_0.95_sorted.txt
python ~/Jensen_response/scripts/averageLD_realData2.py ~/Jensen_response/analysis/${model}_neutrality_locusLen11000_LD_0.05_0.95_sorted.txt ~/Jensen_response/analysis/${model}_neutrality_locusLen11000_LD_0.05_0.95_averaged.txt
done
###########################
# Long range LD DGRP data #
###########################
cat ~/Jensen_response/analysis/DGRP_LD/DGRP_LD_1000_50_2_2_0_chr2R.txt ~/Jensen_response/analysis/DGRP_LD/DGRP_LD_1000_50_2_2_0_chr2L.txt ~/Jensen_response/analysis/DGRP_LD/DGRP_LD_1000_50_2_2_0_chr3R.txt ~/Jensen_response/analysis/DGRP_LD/DGRP_LD_1000_50_2_2_0_chr3L.txt | sort -k2,2 -g > alldataSorted
python ~/Jensen_response/scripts/averageLD_realData.py alldataSorted alldata_averaged
cat alldata_averaged | sort -k1,1 -g > alldata_averaged_sorted
# to do:
# average long range LD
# plot long range LD
# plot S, Pi, TajD, H12, and H2/H1
# Plot hard vs soft sweeps for left over scenarios
# More plots on varying growth rates.