#!/bin/sh
#
#SBATCH --job-name=treemix # Job Name
#SBATCH --nodes=1 # 40 nodes
#SBATCH --ntasks-per-node=1 # 40 CPU allocation per Task
#SBATCH --partition=sixhour # Name of the Slurm partition used
#SBATCH --chdir=/home/d669d153/work/phil.dicaeum/treemix # Set working d$
#SBATCH --mem-per-cpu=5gb # memory requested
#SBATCH --time=100
#convert vcf into treemix file
#/home/d669d153/work/stacks-2.41/populations --in_vcf filtered.85.unlinked.vcf -O . --treemix -M treemix.popmap.txt
#remove stacks header
#echo "$(tail -n +2 filtered.85.unlinked.p.treemix)" > filtered.85.unlinked.p.treemix
#gzip file for input to treemix
#gzip filtered.85.unlinked.p.treemix
#run treemix with m0
/panfs/pfs.local/work/bi/bin/treemix-1.13/src/treemix -i filtered.85.unlinked.p.treemix.gz -root nigrilore -o treem0
#add 1 migration edge
/panfs/pfs.local/work/bi/bin/treemix-1.13/src/treemix -i filtered.85.unlinked.p.treemix.gz -m 1 -g treem0.vertices.gz treem0.edges.gz -o treem1
#add 2 migration edges
/panfs/pfs.local/work/bi/bin/treemix-1.13/src/treemix -i filtered.85.unlinked.p.treemix.gz -m 1 -g treem1.vertices.gz treem1.edges.gz -o treem2
#add 3 migration edges
/panfs/pfs.local/work/bi/bin/treemix-1.13/src/treemix -i filtered.85.unlinked.p.treemix.gz -m 1 -g treem2.vertices.gz treem2.edges.gz -o treem3Dicaeum TreeMix analysis
0.1 Step 1: move ulinked SNP dataset (as vcf) to the cluster, along with a popmap (text file) assigning all samples to tips.
popmap is a tab delimited text file with two columns, ‘sample ID’ and ‘population’.
0.2 Step 2: On the KU HPCC, use this job script to convert vcf to treemix input, then run TreeMix, and progressively add migration edges.
0.3 Step 3: copy the entire treemix outdirectory from the KU cluster to my local machine
scp -r d669d153@hpc.crc.ku.edu:/home/d669d153/work/phil.dicaeum/treemix /Users/devder/Desktop/phil.dicaeum/0.4 Step 4: move into the treemix output directory and plot trees
#source plotting functions that are distributed with treemix
source("~/Downloads/plotting_funcs.R")
#0 edge
plot_tree("~/Desktop/phil.dicaeum/treemix/treem0")Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem0"): NAs introduced
by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem0"): NAs introduced
by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem0"): NAs introduced
by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem0"): NAs introduced
by coercion V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 mind NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 16 4 1 1 1
3 3 zam NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 luz NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
5 15 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 2 2 3 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 15 1 16 3
V11
1 mind:0.000966099
2 (luz:0.00241644,mind:0.000966099):0.00171276
3 zam:0.0147955
4 luz:0.00241644
5 nigrilore:0.0738336
6 ((luz:0.00241644,mind:0.000966099):0.00171276,zam:0.0147955):0.0738336
7 (nigrilore:0.0738336,((luz:0.00241644,mind:0.000966099):0.00171276,zam:0.0147955):0.0738336);
x y ymin ymax
1 0.07651246 0.375 0.25 0.50
2 0.07554636 0.500 0.25 0.75
3 0.08862910 0.125 0.00 0.25
4 0.07796280 0.625 0.50 0.75
5 0.07383360 0.875 0.75 1.00
6 0.07383360 0.250 0.00 0.75
7 0.00000000 0.750 0.00 1.00Warning in max(e[e[, 5] == "MIG", 4]): no non-missing arguments to max;
returning -Inf
[1] 0.07651246 0.08862910 0.07796280 0.07383360
[1] 0.003
[1] "mse 0.0010874170625"$d
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 mind NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 16 4 1 1 1
3 3 zam NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 luz NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
5 15 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 2 2 3 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 15 1 16 3
V11
1 mind:0.000966099
2 (luz:0.00241644,mind:0.000966099):0.00171276
3 zam:0.0147955
4 luz:0.00241644
5 nigrilore:0.0738336
6 ((luz:0.00241644,mind:0.000966099):0.00171276,zam:0.0147955):0.0738336
7 (nigrilore:0.0738336,((luz:0.00241644,mind:0.000966099):0.00171276,zam:0.0147955):0.0738336);
x y ymin ymax
1 0.07651246 0.375 0.25 0.50
2 0.07554636 0.500 0.25 0.75
3 0.08862910 0.125 0.00 0.25
4 0.07796280 0.625 0.50 0.75
5 0.07383360 0.875 0.75 1.00
6 0.07383360 0.250 0.00 0.75
7 0.00000000 0.750 0.00 1.00
$e
V1 V2 V3 V4 V5 V6 V7
1 2 4 0.002416440 1 NOT_MIG 0 1
2 2 1 0.000966099 1 NOT_MIG 0 1
3 16 2 0.001712760 1 NOT_MIG 0 1
4 17 15 0.073833600 1 NOT_MIG 0 1
5 17 16 0.073833600 1 NOT_MIG 0 1
6 16 3 0.014795500 1 NOT_MIG 0 1#1 edge
plot_tree("~/Desktop/phil.dicaeum/treemix/treem1", plus = 0.02, arrow=.1, ybar = 0, scale=F, lwd=1.5)Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem1", plus = 0.02, :
NAs introduced by coercionWarning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem1", plus = 0.02, :
NAs introduced by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem1", plus = 0.02, :
NAs introduced by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem1", plus = 0.02, :
NAs introduced by coercion V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 mind NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 16 1 1 4 1
3 3 zam NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 luz NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
5 15 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 2 2 3 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 15 1 16 3
8 18 <NA> NOT_ROOT MIG NOT_TIP 2 4 NA NA NA
V11
1 mind:0.00096154
2 (mind:0.00096154,luz:0.002421):0.00190811
3 zam:0.0184821
4 luz:0.002421
5 nigrilore:0.0737359
6 ((mind:0.00096154,luz:0.002421):0.00190811,zam:0.0184821):0.0737359
7 (nigrilore:0.0737359,((mind:0.00096154,luz:0.002421):0.00190811,zam:0.0184821):0.0737359);
8 <NA>
x y ymin ymax
1 0.07660555 0.625 0.50 0.75
2 0.07564401 0.500 0.25 0.75
3 0.08870670 0.125 0.00 0.25
4 0.07806502 0.375 0.25 0.50
5 0.07373590 0.875 0.75 1.00
6 0.07373590 0.250 0.00 0.75
7 0.00000000 0.750 0.00 1.00
8 0.07573520 NA 0.25 0.50
[1] "0.037665 0.07564401 0.078065017"
[1] 0.07660555 0.08870670 0.07806502 0.07373590
[1] 0.003$d
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 mind NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 16 1 1 4 1
3 3 zam NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 luz NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
5 15 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 2 2 3 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 15 1 16 3
8 18 <NA> NOT_ROOT MIG NOT_TIP 2 4 NA NA NA
V11
1 mind:0.00096154
2 (mind:0.00096154,luz:0.002421):0.00190811
3 zam:0.0184821
4 luz:0.002421
5 nigrilore:0.0737359
6 ((mind:0.00096154,luz:0.002421):0.00190811,zam:0.0184821):0.0737359
7 (nigrilore:0.0737359,((mind:0.00096154,luz:0.002421):0.00190811,zam:0.0184821):0.0737359);
8 <NA>
x y ymin ymax
1 0.07660555 0.6250000 0.50 0.75
2 0.07564401 0.5000000 0.25 0.75
3 0.08870670 0.1250000 0.00 0.25
4 0.07806502 0.3750000 0.25 0.50
5 0.07373590 0.8750000 0.75 1.00
6 0.07373590 0.2500000 0.00 0.75
7 0.00000000 0.7500000 0.00 1.00
8 0.07573520 0.4952919 0.25 0.50
$e
V1 V2 V3 V4 V5 V6 V7
1 2 1 0.000961540 1.0000000 NOT_MIG 1.000000 1.000000
2 16 2 0.001908110 1.0000000 NOT_MIG 1.000000 1.000000
3 17 15 0.073735900 1.0000000 NOT_MIG 1.000000 1.000000
4 17 16 0.073735900 1.0000000 NOT_MIG 1.000000 1.000000
5 16 3 0.014970800 0.9000090 NOT_MIG 1.000000 0.810017
6 2 18 0.000091187 1.0000000 NOT_MIG 1.000000 1.000000
7 18 4 0.002329820 1.0000000 NOT_MIG 0.037665 1.000000
8 18 3 0.000000000 0.0999906 MIG 0.037665 0.000000#2 edges
plot_tree("~/Desktop/phil.dicaeum/treemix/treem2")Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem2"): NAs introduced
by coercionWarning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem2"): NAs introduced
by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem2"): NAs introduced
by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem2"): NAs introduced
by coercion V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 mind NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 16 1 1 4 1
3 3 zam NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 luz NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
5 15 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 2 2 3 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 15 1 16 3
8 18 <NA> NOT_ROOT MIG NOT_TIP 2 4 NA NA NA
V11
1 mind:0.00096154
2 (mind:0.00096154,luz:0.002421):0.0019081
3 zam:0.0184819
4 luz:0.002421
5 nigrilore:0.0737359
6 ((mind:0.00096154,luz:0.002421):0.0019081,zam:0.0184819):0.0737359
7 (nigrilore:0.0737359,((mind:0.00096154,luz:0.002421):0.0019081,zam:0.0184819):0.0737359);
8 <NA>
x y ymin ymax
1 0.07660554 0.625 0.50 0.75
2 0.07564400 0.500 0.25 0.75
3 0.08870670 0.125 0.00 0.25
4 0.07806500 0.375 0.25 0.50
5 0.07373590 0.875 0.75 1.00
6 0.07373590 0.250 0.00 0.75
7 0.00000000 0.750 0.00 1.00
8 0.07573519 NA 0.25 0.50
[1] "0.0376667 0.075644 0.0780650013"
[1] 0.07660554 0.08870670 0.07806500 0.07373590
[1] 0.003
[1] "mse 0.0010874170625"$d
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 mind NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 16 1 1 4 1
3 3 zam NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 luz NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
5 15 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 2 2 3 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 15 1 16 3
8 18 <NA> NOT_ROOT MIG NOT_TIP 2 4 NA NA NA
V11
1 mind:0.00096154
2 (mind:0.00096154,luz:0.002421):0.0019081
3 zam:0.0184819
4 luz:0.002421
5 nigrilore:0.0737359
6 ((mind:0.00096154,luz:0.002421):0.0019081,zam:0.0184819):0.0737359
7 (nigrilore:0.0737359,((mind:0.00096154,luz:0.002421):0.0019081,zam:0.0184819):0.0737359);
8 <NA>
x y ymin ymax
1 0.07660554 0.6250000 0.50 0.75
2 0.07564400 0.5000000 0.25 0.75
3 0.08870670 0.1250000 0.00 0.25
4 0.07806500 0.3750000 0.25 0.50
5 0.07373590 0.8750000 0.75 1.00
6 0.07373590 0.2500000 0.00 0.75
7 0.00000000 0.7500000 0.00 1.00
8 0.07573519 0.4952917 0.25 0.50
$e
V1 V2 V3 V4 V5 V6 V7
1 2 1 9.61540e-04 1.0000000 NOT_MIG 1.0000000 1.000000
2 16 2 1.90810e-03 1.0000000 NOT_MIG 1.0000000 1.000000
3 17 15 7.37359e-02 1.0000000 NOT_MIG 1.0000000 1.000000
4 17 16 7.37359e-02 1.0000000 NOT_MIG 1.0000000 1.000000
5 16 3 1.49708e-02 0.9000140 NOT_MIG 1.0000000 0.810025
6 2 18 9.11913e-05 1.0000000 NOT_MIG 1.0000000 1.000000
7 18 4 2.32981e-03 1.0000000 NOT_MIG 0.0376667 1.000000
8 18 3 0.00000e+00 0.0999859 MIG 0.0376667 0.000000#3 edges
plot_tree("~/Desktop/phil.dicaeum/treemix/treem3")Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem3"): NAs introduced
by coercionWarning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem3"): NAs introduced
by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem3"): NAs introduced
by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix/treem3"): NAs introduced
by coercion V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 mind NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 16 1 1 4 1
3 3 zam NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 luz NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
5 15 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 2 2 3 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 15 1 16 3
8 18 <NA> NOT_ROOT MIG NOT_TIP 2 4 NA NA NA
V11
1 mind:0.00096154
2 (mind:0.00096154,luz:0.002421):0.00190809
3 zam:0.0184817
4 luz:0.002421
5 nigrilore:0.0737359
6 ((mind:0.00096154,luz:0.002421):0.00190809,zam:0.0184817):0.0737359
7 (nigrilore:0.0737359,((mind:0.00096154,luz:0.002421):0.00190809,zam:0.0184817):0.0737359);
8 <NA>
x y ymin ymax
1 0.07660553 0.625 0.50 0.75
2 0.07564399 0.500 0.25 0.75
3 0.08870671 0.125 0.00 0.25
4 0.07806500 0.375 0.25 0.50
5 0.07373590 0.875 0.75 1.00
6 0.07373590 0.250 0.00 0.75
7 0.00000000 0.750 0.00 1.00
8 0.07573519 NA 0.25 0.50
[1] "0.0376685 0.07564399 0.0780649956"
[1] 0.07660553 0.08870671 0.07806500 0.07373590
[1] 0.003
[1] "mse 0.0010874170625"$d
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 mind NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 16 1 1 4 1
3 3 zam NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 luz NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
5 15 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 2 2 3 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 15 1 16 3
8 18 <NA> NOT_ROOT MIG NOT_TIP 2 4 NA NA NA
V11
1 mind:0.00096154
2 (mind:0.00096154,luz:0.002421):0.00190809
3 zam:0.0184817
4 luz:0.002421
5 nigrilore:0.0737359
6 ((mind:0.00096154,luz:0.002421):0.00190809,zam:0.0184817):0.0737359
7 (nigrilore:0.0737359,((mind:0.00096154,luz:0.002421):0.00190809,zam:0.0184817):0.0737359);
8 <NA>
x y ymin ymax
1 0.07660553 0.6250000 0.50 0.75
2 0.07564399 0.5000000 0.25 0.75
3 0.08870671 0.1250000 0.00 0.25
4 0.07806500 0.3750000 0.25 0.50
5 0.07373590 0.8750000 0.75 1.00
6 0.07373590 0.2500000 0.00 0.75
7 0.00000000 0.7500000 0.00 1.00
8 0.07573519 0.4952914 0.25 0.50
$e
V1 V2 V3 V4 V5 V6 V7
1 2 1 0.0009615400 1.0000000 NOT_MIG 1.0000000 1.000000
2 16 2 0.0019080900 1.0000000 NOT_MIG 1.0000000 1.000000
3 17 15 0.0737359000 1.0000000 NOT_MIG 1.0000000 1.000000
4 17 16 0.0737359000 1.0000000 NOT_MIG 1.0000000 1.000000
5 16 3 0.0149708091 0.9000190 NOT_MIG 1.0000000 0.810034
6 2 18 0.0000911956 1.0000000 NOT_MIG 1.0000000 1.000000
7 18 4 0.0023298100 1.0000000 NOT_MIG 0.0376685 1.000000
8 18 3 0.0000000000 0.0999812 MIG 0.0376685 0.0000000.5 Step 5: evaluate the optimal number of migration edges
#plot to see how much variance is explained by each edge
m=NULL
for(i in 0:3){
m[i+1] <- get_f(paste0("~/Desktop/phil.dicaeum/treemix/treem",i))
}
#print variance explained by each tree with 0,1,2,3 mig edges
m[1] 1 1 1 1#plot
plot(seq(0,3),m,pch="*",cex=2,col="blue", type="b",xlab="migration edge number", ylab="% explained variance")
0.6 Repeat procedure after removing the admixed samples from Mt. Busa
library(vcfR)
***** *** vcfR *** *****
This is vcfR 1.14.0
browseVignettes('vcfR') # Documentation
citation('vcfR') # Citation
***** ***** ***** *****library(SNPfiltR)This is SNPfiltR v.1.0.2
Detailed usage information is available at: devonderaad.github.io/SNPfiltR/
If you use SNPfiltR in your published work, please cite the following papers:
DeRaad, D.A. (2022), SNPfiltR: an R package for interactive and reproducible SNP filtering. Molecular Ecology Resources, 22, 2443-2453. http://doi.org/10.1111/1755-0998.13618
Knaus, Brian J., and Niklaus J. Grunwald. 2017. VCFR: a package to manipulate and visualize variant call format data in R. Molecular Ecology Resources, 17.1:44-53. http://doi.org/10.1111/1755-0998.12549#read in vcf
v<-read.vcfR("~/Desktop/phil.dicaeum/filtered.85.unlinked.vcf")Scanning file to determine attributes.
File attributes:
meta lines: 14
header_line: 15
variant count: 2590
column count: 69
Meta line 14 read in.
All meta lines processed.
gt matrix initialized.
Character matrix gt created.
Character matrix gt rows: 2590
Character matrix gt cols: 69
skip: 0
nrows: 2590
row_num: 0
Processed variant 1000
Processed variant 2000
Processed variant: 2590
All variants processedvcfR.trim<-v[,colnames(v@gt) != "D_hypoleucum_1956" & colnames(v@gt) != "D_hypoleucum_2208" & colnames(v@gt) != "D_hypoleucum_2067" & colnames(v@gt) != "D_hypoleucum_2229" & colnames(v@gt) != "D_hypoleucum_2253"]
vcfR.trim***** Object of Class vcfR *****
55 samples
2590 CHROMs
2,590 variants
Object size: 12 Mb
11.36 percent missing data
***** ***** *****vcfR.trim<-min_mac(vcfR.trim, min.mac = 1)4.36% of SNPs fell below a minor allele count of 1 and were removed from the VCF
#vcfR::write.vcf(vcfR.trim, file="~/Desktop/phil.dicaeum/unlinked.nomtbusa.vcf.gz")0.7 On the KU HPCC, use this job script to convert the downsampled vcf to treemix input, then run TreeMix, and progressively add migration edges.
#!/bin/sh
#
#SBATCH --job-name=treemix # Job Name
#SBATCH --nodes=1 # 40 nodes
#SBATCH --ntasks-per-node=1 # 40 CPU allocation per Task
#SBATCH --partition=sixhour # Name of the Slurm partition used
#SBATCH --chdir=/home/d669d153/work/phil.dicaeum/treemix.nomtbusa # Set working d$
#SBATCH --mem-per-cpu=5gb # memory requested
#SBATCH --time=100
#convert vcf into treemix file
/home/d669d153/work/stacks-2.41/populations --in_vcf unlinked.nomtbusa.vcf -O . --treemix -M treemix.popmap.txt
#remove stacks header
echo "$(tail -n +2 unlinked.nomtbusa.p.treemix)" > filtered.85.unlinked.p.treemix
#gzip file for input to treemix
gzip filtered.85.unlinked.p.treemix
#run treemix with m0
/panfs/pfs.local/work/bi/bin/treemix-1.13/src/treemix -i filtered.85.unlinked.p.treemix.gz -root nigrilore -o treem0
#add 1 migration edge
/panfs/pfs.local/work/bi/bin/treemix-1.13/src/treemix -i filtered.85.unlinked.p.treemix.gz -m 1 -g treem0.vertices.gz treem0.edges.gz -o treem1
#add 2 migration edges
/panfs/pfs.local/work/bi/bin/treemix-1.13/src/treemix -i filtered.85.unlinked.p.treemix.gz -m 1 -g treem1.vertices.gz treem1.edges.gz -o treem2
#add 3 migration edges
/panfs/pfs.local/work/bi/bin/treemix-1.13/src/treemix -i filtered.85.unlinked.p.treemix.gz -m 1 -g treem2.vertices.gz treem2.edges.gz -o treem30.8 Move into the treemix output directory and plot trees
#source plotting functions that are distributed with treemix
#source("~/Downloads/plotting_funcs.R")
#0 edge
plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem0")Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem0"): NAs
introduced by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem0"): NAs
introduced by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem0"): NAs
introduced by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem0"): NAs
introduced by coercion V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 zam NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 1 1 16 2
3 3 mind NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
5 15 luz NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 2 15 1 3 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 4 1 2 3
V11
1 zam:0.0155116
2 (zam:0.0155116,(luz:0.00238025,mind:0.00147542):0.00193747):0.0771807
3 mind:0.00147542
4 nigrilore:0.0771807
5 luz:0.00238025
6 (luz:0.00238025,mind:0.00147542):0.00193747
7 (nigrilore:0.0771807,(zam:0.0155116,(luz:0.00238025,mind:0.00147542):0.00193747):0.0771807);
x y ymin ymax
1 0.09269230 0.625 0.50 0.75
2 0.07718070 0.500 0.00 0.75
3 0.08059359 0.125 0.00 0.25
4 0.07718070 0.875 0.75 1.00
5 0.08149842 0.375 0.25 0.50
6 0.07911817 0.250 0.00 0.50
7 0.00000000 0.750 0.00 1.00Warning in max(e[e[, 5] == "MIG", 4]): no non-missing arguments to max;
returning -Inf
[1] 0.09269230 0.08059359 0.07718070 0.08149842
[1] 0.003
[1] "mse 0.001138198875"$d
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 zam NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 1 1 16 2
3 3 mind NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
5 15 luz NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 2 15 1 3 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 4 1 2 3
V11
1 zam:0.0155116
2 (zam:0.0155116,(luz:0.00238025,mind:0.00147542):0.00193747):0.0771807
3 mind:0.00147542
4 nigrilore:0.0771807
5 luz:0.00238025
6 (luz:0.00238025,mind:0.00147542):0.00193747
7 (nigrilore:0.0771807,(zam:0.0155116,(luz:0.00238025,mind:0.00147542):0.00193747):0.0771807);
x y ymin ymax
1 0.09269230 0.625 0.50 0.75
2 0.07718070 0.500 0.00 0.75
3 0.08059359 0.125 0.00 0.25
4 0.07718070 0.875 0.75 1.00
5 0.08149842 0.375 0.25 0.50
6 0.07911817 0.250 0.00 0.50
7 0.00000000 0.750 0.00 1.00
$e
V1 V2 V3 V4 V5 V6 V7
1 2 1 0.01551160 1 NOT_MIG 0 1
2 16 15 0.00238025 1 NOT_MIG 0 1
3 2 16 0.00193747 1 NOT_MIG 0 1
4 17 4 0.07718070 1 NOT_MIG 0 1
5 17 2 0.07718070 1 NOT_MIG 0 1
6 16 3 0.00147542 1 NOT_MIG 0 1#1 edge
plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem1", plus = 0.02, arrow=.1, ybar = 0, scale=F, lwd=1.5)Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem1", plus =
0.02, : NAs introduced by coercionWarning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem1", plus =
0.02, : NAs introduced by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem1", plus =
0.02, : NAs introduced by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem1", plus =
0.02, : NAs introduced by coercion V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 zam NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 1 1 16 2
3 3 mind NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
5 15 luz NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 2 3 1 15 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 4 1 2 3
8 18 <NA> NOT_ROOT MIG NOT_TIP 16 15 NA NA NA
V11
1 zam:0.0194395
2 (zam:0.0194395,(mind:0.00142987,luz:0.00242581):0.00220334):0.0770478
3 mind:0.00142987
4 nigrilore:0.0770478
5 luz:0.00242581
6 (mind:0.00142987,luz:0.00242581):0.00220334
7 (nigrilore:0.0770478,(zam:0.0194395,(mind:0.00142987,luz:0.00242581):0.00220334):0.0770478);
8 <NA>
x y ymin ymax
1 0.09279411 0.625 0.50 0.75
2 0.07704780 0.500 0.00 0.75
3 0.08068101 0.375 0.25 0.50
4 0.07704780 0.875 0.75 1.00
5 0.08167695 0.125 0.00 0.25
6 0.07925114 0.250 0.00 0.50
7 0.00000000 0.750 0.00 1.00
8 0.08016231 NA 0.00 0.25
[1] "0.375613 0.07925114 0.081676945"
[1] 0.09279411 0.08068101 0.07704780 0.08167695
[1] 0.003$d
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 zam NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 1 1 16 2
3 3 mind NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
5 15 luz NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 2 3 1 15 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 4 1 2 3
8 18 <NA> NOT_ROOT MIG NOT_TIP 16 15 NA NA NA
V11
1 zam:0.0194395
2 (zam:0.0194395,(mind:0.00142987,luz:0.00242581):0.00220334):0.0770478
3 mind:0.00142987
4 nigrilore:0.0770478
5 luz:0.00242581
6 (mind:0.00142987,luz:0.00242581):0.00220334
7 (nigrilore:0.0770478,(zam:0.0194395,(mind:0.00142987,luz:0.00242581):0.00220334):0.0770478);
8 <NA>
x y ymin ymax
1 0.09279411 0.6250000 0.50 0.75
2 0.07704780 0.5000000 0.00 0.75
3 0.08068101 0.3750000 0.25 0.50
4 0.07704780 0.8750000 0.75 1.00
5 0.08167695 0.1250000 0.00 0.25
6 0.07925114 0.2500000 0.00 0.50
7 0.00000000 0.7500000 0.00 1.00
8 0.08016230 0.2030484 0.00 0.25
$e
V1 V2 V3 V4 V5 V6 V7
1 2 1 0.015746310 0.9000090 NOT_MIG 1.000000 0.810017
2 2 16 0.002203340 1.0000000 NOT_MIG 1.000000 1.000000
3 17 4 0.077047800 1.0000000 NOT_MIG 1.000000 1.000000
4 17 2 0.077047800 1.0000000 NOT_MIG 1.000000 1.000000
5 16 3 0.001429870 1.0000000 NOT_MIG 1.000000 1.000000
6 16 18 0.000911165 1.0000000 NOT_MIG 1.000000 1.000000
7 18 15 0.001514640 1.0000000 NOT_MIG 0.375613 1.000000
8 18 1 0.000000000 0.0999906 MIG 0.375613 0.000000#2 edges
plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem2")Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem2"): NAs
introduced by coercionWarning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem2"): NAs
introduced by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem2"): NAs
introduced by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem2"): NAs
introduced by coercion V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 zam NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 1 1 16 2
3 3 mind NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
5 15 luz NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 2 3 1 15 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 4 1 2 3
8 18 <NA> NOT_ROOT MIG NOT_TIP 16 15 NA NA NA
V11
1 zam:0.0194393
2 (zam:0.0194393,(mind:0.00142987,luz:0.00242581):0.00220332):0.0770478
3 mind:0.00142987
4 nigrilore:0.0770478
5 luz:0.00242581
6 (mind:0.00142987,luz:0.00242581):0.00220332
7 (nigrilore:0.0770478,(zam:0.0194393,(mind:0.00142987,luz:0.00242581):0.00220332):0.0770478);
8 <NA>
x y ymin ymax
1 0.09279412 0.625 0.50 0.75
2 0.07704780 0.500 0.00 0.75
3 0.08068099 0.375 0.25 0.50
4 0.07704780 0.875 0.75 1.00
5 0.08167693 0.125 0.00 0.25
6 0.07925112 0.250 0.00 0.50
7 0.00000000 0.750 0.00 1.00
8 0.08016233 NA 0.00 0.25
[1] "0.375631 0.07925112 0.081676927"
[1] 0.09279412 0.08068099 0.07704780 0.08167693
[1] 0.003
[1] "mse 0.001138198875"$d
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 zam NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 1 1 16 2
3 3 mind NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
5 15 luz NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 2 3 1 15 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 4 1 2 3
8 18 <NA> NOT_ROOT MIG NOT_TIP 16 15 NA NA NA
V11
1 zam:0.0194393
2 (zam:0.0194393,(mind:0.00142987,luz:0.00242581):0.00220332):0.0770478
3 mind:0.00142987
4 nigrilore:0.0770478
5 luz:0.00242581
6 (mind:0.00142987,luz:0.00242581):0.00220332
7 (nigrilore:0.0770478,(zam:0.0194393,(mind:0.00142987,luz:0.00242581):0.00220332):0.0770478);
8 <NA>
x y ymin ymax
1 0.09279412 0.6250000 0.50 0.75
2 0.07704780 0.5000000 0.00 0.75
3 0.08068099 0.3750000 0.25 0.50
4 0.07704780 0.8750000 0.75 1.00
5 0.08167693 0.1250000 0.00 0.25
6 0.07925112 0.2500000 0.00 0.50
7 0.00000000 0.7500000 0.00 1.00
8 0.08016233 0.2030461 0.00 0.25
$e
V1 V2 V3 V4 V5 V6 V7
1 2 1 0.015746323 0.9000140 NOT_MIG 1.000000 0.810025
2 2 16 0.002203320 1.0000000 NOT_MIG 1.000000 1.000000
3 17 4 0.077047800 1.0000000 NOT_MIG 1.000000 1.000000
4 17 2 0.077047800 1.0000000 NOT_MIG 1.000000 1.000000
5 16 3 0.001429870 1.0000000 NOT_MIG 1.000000 1.000000
6 16 18 0.000911207 1.0000000 NOT_MIG 1.000000 1.000000
7 18 15 0.001514600 1.0000000 NOT_MIG 0.375631 1.000000
8 18 1 0.000000000 0.0999859 MIG 0.375631 0.000000#3 edges
plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem3")Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem3"): NAs
introduced by coercionWarning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem3"): NAs
introduced by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem3"): NAs
introduced by coercion
Warning in plot_tree("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem3"): NAs
introduced by coercion V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 zam NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 1 1 16 2
3 3 mind NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
5 15 luz NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 2 3 1 15 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 4 1 2 3
8 18 <NA> NOT_ROOT MIG NOT_TIP 16 15 NA NA NA
V11
1 zam:0.019439
2 (zam:0.019439,(mind:0.00142987,luz:0.00242581):0.00220331):0.0770478
3 mind:0.00142987
4 nigrilore:0.0770478
5 luz:0.00242581
6 (mind:0.00142987,luz:0.00242581):0.00220331
7 (nigrilore:0.0770478,(zam:0.019439,(mind:0.00142987,luz:0.00242581):0.00220331):0.0770478);
8 <NA>
x y ymin ymax
1 0.09279405 0.625 0.50 0.75
2 0.07704780 0.500 0.00 0.75
3 0.08068098 0.375 0.25 0.50
4 0.07704780 0.875 0.75 1.00
5 0.08167692 0.125 0.00 0.25
6 0.07925111 0.250 0.00 0.50
7 0.00000000 0.750 0.00 1.00
8 0.08016236 NA 0.00 0.25
[1] "0.375648 0.07925111 0.08167692"
[1] 0.09279405 0.08068098 0.07704780 0.08167692
[1] 0.003
[1] "mse 0.001138198875"$d
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1 1 zam NOT_ROOT NOT_MIG TIP 2 NA NA NA NA
2 2 <NA> NOT_ROOT NOT_MIG NOT_TIP 17 1 1 16 2
3 3 mind NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
4 4 nigrilore NOT_ROOT NOT_MIG TIP 17 NA NA NA NA
5 15 luz NOT_ROOT NOT_MIG TIP 16 NA NA NA NA
6 16 <NA> NOT_ROOT NOT_MIG NOT_TIP 2 3 1 15 1
7 17 <NA> ROOT NOT_MIG NOT_TIP 17 4 1 2 3
8 18 <NA> NOT_ROOT MIG NOT_TIP 16 15 NA NA NA
V11
1 zam:0.019439
2 (zam:0.019439,(mind:0.00142987,luz:0.00242581):0.00220331):0.0770478
3 mind:0.00142987
4 nigrilore:0.0770478
5 luz:0.00242581
6 (mind:0.00142987,luz:0.00242581):0.00220331
7 (nigrilore:0.0770478,(zam:0.019439,(mind:0.00142987,luz:0.00242581):0.00220331):0.0770478);
8 <NA>
x y ymin ymax
1 0.09279405 0.625000 0.50 0.75
2 0.07704780 0.500000 0.00 0.75
3 0.08068098 0.375000 0.25 0.50
4 0.07704780 0.875000 0.75 1.00
5 0.08167692 0.125000 0.00 0.25
6 0.07925111 0.250000 0.00 0.50
7 0.00000000 0.750000 0.00 1.00
8 0.08016236 0.203044 0.00 0.25
$e
V1 V2 V3 V4 V5 V6 V7
1 2 1 0.01574625 0.9000190 NOT_MIG 1.000000 0.810034
2 2 16 0.00220331 1.0000000 NOT_MIG 1.000000 1.000000
3 17 4 0.07704780 1.0000000 NOT_MIG 1.000000 1.000000
4 17 2 0.07704780 1.0000000 NOT_MIG 1.000000 1.000000
5 16 3 0.00142987 1.0000000 NOT_MIG 1.000000 1.000000
6 16 18 0.00091125 1.0000000 NOT_MIG 1.000000 1.000000
7 18 15 0.00151456 1.0000000 NOT_MIG 0.375648 1.000000
8 18 1 0.00000000 0.0999812 MIG 0.375648 0.0000000.9 Evaluate the optimal number of migration edges
#plot to see how much variance is explained by each edge
m=NULL
for(i in 0:3){
m[i+1] <- get_f(paste0("~/Desktop/phil.dicaeum/treemix.nomtbusa/treem",i))
}
#print variance explained by each tree with 0,1,2,3 mig edges
m[1] 0.9999993 1.0000000 1.0000000 1.0000000#plot
plot(seq(0,3),m,pch="*",cex=2,col="blue", type="b",xlab="migration edge number", ylab="% explained variance")
0.10 Bootstrap the optimal tree with all samples included
#
#SBATCH \--job-name=treemix \# Job Name
#SBATCH \--nodes=1 \# 40 nodes
#SBATCH \--ntasks-per-node=1 \# 40 CPU allocation per Task
#SBATCH \--partition=bi \# Name of the Slurm partition used
#SBATCH \--chdir=/home/d669d153/work/phil.dicaeum/treemix/boots \# Set working d\$
#SBATCH \--mem-per-cpu=1000 \# memory requested
#SBATCH \--time=200
#100 bootstraps over 100 SNP blocks with 1 migration edge
for i in {1..100}; do
/panfs/pfs.local/work/bi/bin/treemix-1.13/src/treemix -i /home/d669d153/work/phil.dicaeum/treemix/filtered.85.unlinked.p.treemix.gz -m 1 -g /home/d669d153/work/phil.dicaeum/treemix/treem1.vertices.gz /home/d669d153/work/phil.dicaeum/treemix/treem1.edges.gz -bootstrap -k 100 -o \$i.treemix
done;
# unzip the tree files
for i in { ls \*treeout.gz }; do
gzip -d \$i
done;
# in R:
#summarize bootstraps into a single file called 'bootstraps.100.trees'
module load R
R \--no-save
setwd("/home/d669d153/work/phil.dicaeum/treemix/boots/")
x \<- list.files(pattern="\*treeout")
for(a in 1:length(x)) {
if (a==1) {
output \<- scan(x\[a\], what="character")\[1\]
} else {
output \<- c(output, scan(x\[a\], what="character")\[1\])
}
}
write(output, file="bootstraps.100.trees", ncolumns=1)
quit()
# in bash
# summarize bootstrap support info from the file '100.bootstraps.trees' as internal branch annotations on a summary tree called 'boots.summed.tre'
/home/d669d153/work/DendroPy/applications/sumtrees/sumtrees.py \--output=boots.summed.tre \--min-clade-freq=0.05 bootstraps.100.trees