Dicaeum splitstree
Devon DeRaad
2022-11-30
read in your filtered vcf file
library(vcfR)##
## ***** *** vcfR *** *****
## This is vcfR 1.12.0
## browseVignettes('vcfR') # Documentation
## citation('vcfR') # Citation
## ***** ***** ***** *****library(adegenet)## Loading required package: ade4##
## /// adegenet 2.1.5 is loaded ////////////
##
## > overview: '?adegenet'
## > tutorials/doc/questions: 'adegenetWeb()'
## > bug reports/feature requests: adegenetIssues()library(StAMPP)## Loading required package: pegas## Loading required package: ape## Registered S3 method overwritten by 'pegas':
## method from
## print.amova ade4##
## Attaching package: 'pegas'## The following object is masked from 'package:ape':
##
## mst## The following object is masked from 'package:ade4':
##
## amova## The following objects are masked from 'package:vcfR':
##
## getINFO, write.vcf#read vcf
v<-read.vcfR("~/Desktop/phil.dicaeum/filtered.85.vcf")## Scanning file to determine attributes.
## File attributes:
## meta lines: 14
## header_line: 15
## variant count: 18499
## column count: 69
##
Meta line 14 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 18499
## Character matrix gt cols: 69
## skip: 0
## nrows: 18499
## row_num: 0
##
Processed variant 1000
Processed variant 2000
Processed variant 3000
Processed variant 4000
Processed variant 5000
Processed variant 6000
Processed variant 7000
Processed variant 8000
Processed variant 9000
Processed variant 10000
Processed variant 11000
Processed variant 12000
Processed variant 13000
Processed variant 14000
Processed variant 15000
Processed variant 16000
Processed variant 17000
Processed variant 18000
Processed variant: 18499
## All variants processedprepare input matrix
#convert vcfR to genlight
gen<-vcfR2genlight(v)
#check sample names (remember the splitstree gui won't let you visualize more than 10 characters in a sample name, and all sample names must be unique)
gen@ind.names## [1] "D_hypoleucum_18070" "D_hypoleucum_18191"
## [3] "D_hypoleucum_14037" "D_hypoleucum_14079"
## [5] "D_hypoleucum_14065" "D_hypoleucum_14120"
## [7] "D_hypoleucum_18159" "D_hypoleucum_14075"
## [9] "D_hypoleucum_20218" "D_hypoleucum_19177"
## [11] "D_hypoleucum_19178" "D_hypoleucum_20213"
## [13] "D_hypoleucum_FMNH454949" "D_hypoleucum_25637"
## [15] "D_hypoleucum_1271" "D_hypoleucum_1273"
## [17] "D_hypoleucum_1275" "D_hypoleucum_25921"
## [19] "D_hypoleucum_3274" "D_hypoleucum_26984"
## [21] "D_hypoleucum_3208" "D_hypoleucum_3158"
## [23] "D_hypoleucum_2253" "D_hypoleucum_3095"
## [25] "D_hypoleucum_462070" "D_hypoleucum_25675"
## [27] "D_hypoleucum_454950" "D_hypoleucum_25880"
## [29] "D_hypoleucum_3314" "D_hypoleucum_357608"
## [31] "D_hypoleucum_357615" "D_hypoleucum_357612"
## [33] "D_hypoleucum_19638" "D_hypoleucum_25868"
## [35] "D_hypoleucum_17969" "D_hypoleucum_25672"
## [37] "D_hypoleucum_26975" "D_hypoleucum_2229"
## [39] "D_hypoleucum_2067" "D_hypoleucum_28329"
## [41] "D_hypoleucum_28361" "D_hypoleucum_28376"
## [43] "D_hypoleucum_2208" "D_hypoleucum_28294"
## [45] "D_hypoleucum_1956" "D_hypoleucum_25670"
## [47] "D_hypoleucum_20921" "D_hypoleucum_20193"
## [49] "D_hypoleucum_28416" "D_hypoleucum_18193"
## [51] "D_hypoleucum_27450" "D_hypoleucum_19046"
## [53] "D_hypoleucum_19136" "D_hypoleucum_28676"
## [55] "D_nigrilore_KU28413" "D_nigrilore_KU28414"
## [57] "D_hypoleucum_29945" "D_hypoleucum_31636"
## [59] "D_hypoleucum_29951" "D_hypoleucum_31644"#edit your sample names to fit this criteria
gen@ind.names<-gsub("D_hypoleucum","hy", gen@ind.names)
gen@ind.names<-gsub("D_nigrilore","ny", gen@ind.names)
#double check sample names
gen@ind.names## [1] "hy_18070" "hy_18191" "hy_14037" "hy_14079"
## [5] "hy_14065" "hy_14120" "hy_18159" "hy_14075"
## [9] "hy_20218" "hy_19177" "hy_19178" "hy_20213"
## [13] "hy_FMNH454949" "hy_25637" "hy_1271" "hy_1273"
## [17] "hy_1275" "hy_25921" "hy_3274" "hy_26984"
## [21] "hy_3208" "hy_3158" "hy_2253" "hy_3095"
## [25] "hy_462070" "hy_25675" "hy_454950" "hy_25880"
## [29] "hy_3314" "hy_357608" "hy_357615" "hy_357612"
## [33] "hy_19638" "hy_25868" "hy_17969" "hy_25672"
## [37] "hy_26975" "hy_2229" "hy_2067" "hy_28329"
## [41] "hy_28361" "hy_28376" "hy_2208" "hy_28294"
## [45] "hy_1956" "hy_25670" "hy_20921" "hy_20193"
## [49] "hy_28416" "hy_18193" "hy_27450" "hy_19046"
## [53] "hy_19136" "hy_28676" "ny_KU28413" "ny_KU28414"
## [57] "hy_29945" "hy_31636" "hy_29951" "hy_31644"#assign sample names as populations (population assignments are a requirement for the stampp functions, but can be arbitrary here because we want a pairwise divergence matrix among all samples, not samples assigned populations)
pop(gen)<-gen@ind.names
#make pairwise divergence matrix among all samples
sample.div <- stamppNeisD(gen, pop = FALSE)
#export for splitstree
#stamppPhylip(distance.mat=sample.div, file="~/Desktop/phil.dicaeum/dicaeum.85.filtered.splits.txt")
#open this file using the splitstree GUI. Take a screenshot and use the knitr function 'include_graphics' to show the resulting tree in this document alongside the code used to generate it:
knitr::include_graphics("/Users/devder/Desktop/phil.dicaeum/phil.dicaeum.initial.splitstree.png")
fiddle with the results
#I don't like how the outgroup looks (weirdly nested within one of the clades on a super long branch)
#Let's remove the outgroup and then revisualize the splitstree
#remove two nigrilore samples from the vcf
v.sub<-v[,colnames(v@gt) != "D_nigrilore_KU28413" & colnames(v@gt) != "D_nigrilore_KU28414"]redo with no outgroup
#convert vcfR to genlight
gen<-vcfR2genlight(v.sub)
#fix sample names
gen@ind.names<-gsub("D_hypoleucum","hy", gen@ind.names)
#assign sample names as populations
pop(gen)<-gen@ind.names
#make pairwise divergence matrix among all samples
sample.div <- stamppNeisD(gen, pop = FALSE)
#export for splitstree
#stamppPhylip(distance.mat=sample.div, file="~/Desktop/phil.dicaeum/dicaeum.85.nooutgroup.filtered.splits.txt")
#open this file using the splitstree GUI. Take a screenshot and use the knitr function 'include_graphics' to show the resulting tree in this document alongside the code used to generate it:
knitr::include_graphics("/Users/devder/Desktop/phil.dicaeum/splitstree.nooutgroup.png")
Make this presentable for publication
#at this point, I will edit the tree in the splitstree GUI to remove sample names and add dots on the end of each branch, which can be color-coded according to sample locality, giving me a nice, publication quality product
#my final result:
knitr::include_graphics("/Users/devder/Desktop/phil.dicaeum/splitstree.png")