***** *** vcfR *** *****
This is vcfR 1.14.0
browseVignettes('vcfR') # Documentation
citation('vcfR') # Citation
***** ***** ***** *****
Code
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
Code
library(StAMPP)
Loading required package: pegas
Loading required package: ape
Attaching package: 'pegas'
The following object is masked from 'package:ape':
mst
The following objects are masked from 'package:vcfR':
getINFO, write.vcf
Registered S3 method overwritten by 'ade4':
method from
print.amova pegas
Code
library(adegenet)
Loading required package: ade4
Attaching package: 'ade4'
The following object is masked from 'package:pegas':
amova
pop(gen)<-gen@ind.names#assign populations (a StaMPP requirement)gen@pop<-as.factor(gen@ind.names)#generate pairwise divergence matrixsample.div <-stamppNeisD(gen, pop =FALSE)#export for splitstree#stamppPhylip(distance.mat=sample.div, file="~/Desktop/grosbeak.rad/grosbeak.90.splits.txt")knitr::include_graphics(c("/Users/devonderaad/Desktop/grosbeak.rad/90.splits.png"))
3 re-make splitstree with mito and locality info listed in each sample name
Code
#convert to genlightgen<-vcfR2genlight(v)#read in locality and mito infosamps<-read.csv("~/Desktop/grosbeak.data.csv")samps<-samps[samps$passed.genomic.filtering =="TRUE",] #retain only samples that passed filteringsamps$sample_id == gen@ind.names #check if sample info table order matches the vcf
#make splitstree with these updated labelspop(gen)<-gen@ind.names#assign populations (a StaMPP requirement)gen@pop<-as.factor(gen@ind.names)#generate pairwise divergence matrixsample.div <-stamppNeisD(gen, pop =FALSE)#export for splitstree#stamppPhylip(distance.mat=sample.div, file="~/Desktop/grosbeak.rad/grosbeak.mito.site.splits.txt")
4 Calculate Fst between the ends of the transect
Code
#isolate only the parental ends of the transectsv.sub<-v[,c(TRUE,samps$site ==0| samps$site ==12)]#make sure this workedv.sub
***** Object of Class vcfR *****
17 samples
1027 CHROMs
51,595 variants
Object size: 88.7 Mb
4.652 percent missing data
***** ***** *****
Code
#convert vcfR to genlightgen<-vcfR2genlight(v.sub)#assign samples to the three groups shown abovegen@pop<-as.factor(samps$site[samps$site ==0| samps$site ==12])#calculate pairwise Fst using the stampp packagestamppFst(gen)
---title: "Grosbeak.splitstree"format: html: code-fold: show code-tools: truetoc: truetoc-title: Document Contentsnumber-sections: trueembed-resources: true---## load libraries and read in data```{r, results=FALSE}library(vcfR)library(SNPfiltR)library(StAMPP)library(adegenet)library(ggplot2)#read in filev<-read.vcfR("~/Desktop/grosbeak.rad/grosbeak.filtered.snps.vcf.gz")```## make splitstree```{r}#get info for this datasetv#convert to genlightgen<-vcfR2genlight(v)#fix sample names to fit in <= 10 charactersgen@ind.namesgen@ind.names<-gsub("P_hybrid_","hyb", gen@ind.names)gen@ind.names<-gsub("P_ludovicianus_","lud", gen@ind.names)gen@ind.names<-gsub("P_melanocephalus_","mel", gen@ind.names)gen@ind.namespop(gen)<-gen@ind.names#assign populations (a StaMPP requirement)gen@pop<-as.factor(gen@ind.names)#generate pairwise divergence matrixsample.div <-stamppNeisD(gen, pop =FALSE)#export for splitstree#stamppPhylip(distance.mat=sample.div, file="~/Desktop/grosbeak.rad/grosbeak.90.splits.txt")knitr::include_graphics(c("/Users/devonderaad/Desktop/grosbeak.rad/90.splits.png"))```## re-make splitstree with mito and locality info listed in each sample name```{r}#convert to genlightgen<-vcfR2genlight(v)#read in locality and mito infosamps<-read.csv("~/Desktop/grosbeak.data.csv")samps<-samps[samps$passed.genomic.filtering =="TRUE",] #retain only samples that passed filteringsamps$sample_id == gen@ind.names #check if sample info table order matches the vcfsamps<-samps[match(gen@ind.names,samps$sample_id),] #use 'match' to match orderssamps$sample_id == gen@ind.names #check if this worked#fix sample names to fit in <= 10 characters and have locality and mito info and still be uniquegen@ind.namesgen@ind.names<-gsub("P_hybrid_","", gen@ind.names)gen@ind.names<-gsub("P_ludovicianus_","", gen@ind.names)gen@ind.names<-gsub("P_melanocephalus_","", gen@ind.names)gen@ind.names<-paste(samps$site,gen@ind.names,sep="_")gen@ind.names<-paste(samps$mtDNA,gen@ind.names,sep="_")gen@ind.names<-gsub("NA","N", gen@ind.names)gen@ind.names#make splitstree with these updated labelspop(gen)<-gen@ind.names#assign populations (a StaMPP requirement)gen@pop<-as.factor(gen@ind.names)#generate pairwise divergence matrixsample.div <-stamppNeisD(gen, pop =FALSE)#export for splitstree#stamppPhylip(distance.mat=sample.div, file="~/Desktop/grosbeak.rad/grosbeak.mito.site.splits.txt")```## Calculate Fst between the ends of the transect```{r}#isolate only the parental ends of the transectsv.sub<-v[,c(TRUE,samps$site ==0| samps$site ==12)]#make sure this workedv.sub#convert vcfR to genlightgen<-vcfR2genlight(v.sub)#assign samples to the three groups shown abovegen@pop<-as.factor(samps$site[samps$site ==0| samps$site ==12])#calculate pairwise Fst using the stampp packagestamppFst(gen)```
