#Download genomics general (run the following in bash)
#download package
#wget https://github.com/simonhmartin/genomics_general/archive/v0.4.tar.gz
#extract files from zipped archive
#tar -xzf v0.4.tar.gz
#delete zipped file
#rm v0.4.tar.gz
#ensure libraries are recognizable by python
#export PYTHONPATH=$PYTHONPATH:genomics_general-0.4
library(vcfR)
##
## ***** *** vcfR *** *****
## This is vcfR 1.12.0
## browseVignettes('vcfR') # Documentation
## citation('vcfR') # Citation
## ***** ***** ***** *****
#in R run the following code to generate a vcf file with only the samples of interest
#read in vcf as vcfR
vcfR <- read.vcfR("~/Desktop/aph.data/unzipped.filtered.vcf")
## Scanning file to determine attributes.
## File attributes:
## meta lines: 14
## header_line: 15
## variant count: 16307
## column count: 104
##
Meta line 14 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 16307
## Character matrix gt cols: 104
## skip: 0
## nrows: 16307
## 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: 16307
## All variants processed
dim(vcfR)
## variants fix_cols gt_cols
## 16307 8 96
colnames(vcfR@gt)
## [1] "FORMAT" "A_californica_333849" "A_californica_333854"
## [4] "A_californica_333855" "A_californica_333857" "A_californica_333860"
## [7] "A_californica_333862" "A_californica_333871" "A_californica_333873"
## [10] "A_californica_333874" "A_californica_333914" "A_californica_333917"
## [13] "A_californica_333931" "A_californica_333932" "A_californica_333934"
## [16] "A_californica_334002" "A_californica_334012" "A_californica_334015"
## [19] "A_californica_334017" "A_californica_334019" "A_woodhouseii_334171"
## [22] "A_californica_342037" "A_californica_342048" "A_californica_342050"
## [25] "A_californica_342066" "A_californica_342072" "A_californica_343428"
## [28] "A_californica_343430" "A_californica_343432" "A_californica_343438"
## [31] "A_californica_343451" "A_californica_393615" "A_californica_393721"
## [34] "A_coerulescens_396251" "A_coerulescens_396254" "A_coerulescens_396256"
## [37] "A_coerulescens_396263" "A_coerulescens_396264" "A_coerulescens_396265"
## [40] "A_insularis_334025" "A_insularis_334029" "A_insularis_334031"
## [43] "A_insularis_334032" "A_insularis_334033" "A_insularis_334034"
## [46] "A_insularis_334037" "A_insularis_334038" "A_sumichrasti_343502"
## [49] "A_sumichrasti_343503" "A_sumichrasti_343510" "A_sumichrasti_343511"
## [52] "A_sumichrasti_343512" "A_sumichrasti_343514" "A_sumichrasti_343515"
## [55] "A_sumichrasti_393635" "A_sumichrasti_393638" "A_sumichrasti_393640"
## [58] "A_woodhouseii_334062" "A_woodhouseii_334063" "A_woodhouseii_334086"
## [61] "A_woodhouseii_334088" "A_woodhouseii_334096" "A_woodhouseii_334132"
## [64] "A_woodhouseii_334133" "A_woodhouseii_334134" "A_woodhouseii_334142"
## [67] "A_woodhouseii_334153" "A_woodhouseii_334161" "A_woodhouseii_334188"
## [70] "A_woodhouseii_334190" "A_woodhouseii_334196" "A_woodhouseii_334210"
## [73] "A_woodhouseii_334211" "A_woodhouseii_334217" "A_woodhouseii_334230"
## [76] "A_woodhouseii_334241" "A_woodhouseii_334242" "A_woodhouseii_334243"
## [79] "A_woodhouseii_334244" "A_woodhouseii_334247" "A_woodhouseii_334250"
## [82] "A_woodhouseii_343458" "A_woodhouseii_343461" "A_woodhouseii_343476"
## [85] "A_woodhouseii_343480" "A_woodhouseii_343481" "A_woodhouseii_343483"
## [88] "A_woodhouseii_343497" "A_woodhouseii_393605" "A_woodhouseii_393606"
## [91] "A_woodhouseii_393697" "A_woodhouseii_393698" "A_woodhouseii_393702"
## [94] "A_woodhouseii_393712" "A_woodhouseii_393713" "A_woodhouseii_395768"
#retain only ((sumi,int.wood),cali,fla) as ((P1,P2),P3,P4)
vcfR@gt<-vcfR@gt[,c(1:39,48:74)]
dim(vcfR)
## variants fix_cols gt_cols
## 16307 8 66
#retain only SNPs that still have a maf > 0
source("~/Desktop/snpfiltR/min_mac.R")
vcfR<-min_mac(vcfR, min.mac = 1)
## [1] "14.83% of SNPs fell below a minor allele count of 1 and were removed from the VCF"
dim(vcfR)
## variants fix_cols gt_cols
## 13888 8 66
write.vcf(vcfR, file="~/Desktop/aph.data/dfs/dfs.sumi.wood.cali.fla.vcf.gz")
#generate popmap
colnames(vcfR@gt)
## [1] "FORMAT" "A_californica_333849" "A_californica_333854"
## [4] "A_californica_333855" "A_californica_333857" "A_californica_333860"
## [7] "A_californica_333862" "A_californica_333871" "A_californica_333873"
## [10] "A_californica_333874" "A_californica_333914" "A_californica_333917"
## [13] "A_californica_333931" "A_californica_333932" "A_californica_333934"
## [16] "A_californica_334002" "A_californica_334012" "A_californica_334015"
## [19] "A_californica_334017" "A_californica_334019" "A_woodhouseii_334171"
## [22] "A_californica_342037" "A_californica_342048" "A_californica_342050"
## [25] "A_californica_342066" "A_californica_342072" "A_californica_343428"
## [28] "A_californica_343430" "A_californica_343432" "A_californica_343438"
## [31] "A_californica_343451" "A_californica_393615" "A_californica_393721"
## [34] "A_coerulescens_396251" "A_coerulescens_396254" "A_coerulescens_396256"
## [37] "A_coerulescens_396263" "A_coerulescens_396264" "A_coerulescens_396265"
## [40] "A_sumichrasti_343502" "A_sumichrasti_343503" "A_sumichrasti_343510"
## [43] "A_sumichrasti_343511" "A_sumichrasti_343512" "A_sumichrasti_343514"
## [46] "A_sumichrasti_343515" "A_sumichrasti_393635" "A_sumichrasti_393638"
## [49] "A_sumichrasti_393640" "A_woodhouseii_334062" "A_woodhouseii_334063"
## [52] "A_woodhouseii_334086" "A_woodhouseii_334088" "A_woodhouseii_334096"
## [55] "A_woodhouseii_334132" "A_woodhouseii_334133" "A_woodhouseii_334134"
## [58] "A_woodhouseii_334142" "A_woodhouseii_334153" "A_woodhouseii_334161"
## [61] "A_woodhouseii_334188" "A_woodhouseii_334190" "A_woodhouseii_334196"
## [64] "A_woodhouseii_334210" "A_woodhouseii_334211" "A_woodhouseii_334217"
names<-colnames(vcfR@gt)[-1]
pops<-c(rep("P3", times=19),"P2",rep("P3", times=12),rep("OG", times=6),rep("P1",times=10),rep("P2",times=17))
d<-data.frame(names=names,pops=pops)
table(d$pops)
##
## OG P1 P2 P3
## 6 10 18 31
#write popmap
write.table(d, file="~/Desktop/aph.data/dfs/pops.sumi.wood.cali.fla.txt", sep = "\t", row.names = F, col.names = F, quote = F)
#try flipping the nodes and doing ((island,cali),int.wood,fla) as ((P1,P2),P3,P4)
#read in vcf as vcfR
vcfR <- read.vcfR("~/Desktop/aph.data/unzipped.filtered.vcf")
## Scanning file to determine attributes.
## File attributes:
## meta lines: 14
## header_line: 15
## variant count: 16307
## column count: 104
##
Meta line 14 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 16307
## Character matrix gt cols: 104
## skip: 0
## nrows: 16307
## 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: 16307
## All variants processed
dim(vcfR)
## variants fix_cols gt_cols
## 16307 8 96
colnames(vcfR@gt)
## [1] "FORMAT" "A_californica_333849" "A_californica_333854"
## [4] "A_californica_333855" "A_californica_333857" "A_californica_333860"
## [7] "A_californica_333862" "A_californica_333871" "A_californica_333873"
## [10] "A_californica_333874" "A_californica_333914" "A_californica_333917"
## [13] "A_californica_333931" "A_californica_333932" "A_californica_333934"
## [16] "A_californica_334002" "A_californica_334012" "A_californica_334015"
## [19] "A_californica_334017" "A_californica_334019" "A_woodhouseii_334171"
## [22] "A_californica_342037" "A_californica_342048" "A_californica_342050"
## [25] "A_californica_342066" "A_californica_342072" "A_californica_343428"
## [28] "A_californica_343430" "A_californica_343432" "A_californica_343438"
## [31] "A_californica_343451" "A_californica_393615" "A_californica_393721"
## [34] "A_coerulescens_396251" "A_coerulescens_396254" "A_coerulescens_396256"
## [37] "A_coerulescens_396263" "A_coerulescens_396264" "A_coerulescens_396265"
## [40] "A_insularis_334025" "A_insularis_334029" "A_insularis_334031"
## [43] "A_insularis_334032" "A_insularis_334033" "A_insularis_334034"
## [46] "A_insularis_334037" "A_insularis_334038" "A_sumichrasti_343502"
## [49] "A_sumichrasti_343503" "A_sumichrasti_343510" "A_sumichrasti_343511"
## [52] "A_sumichrasti_343512" "A_sumichrasti_343514" "A_sumichrasti_343515"
## [55] "A_sumichrasti_393635" "A_sumichrasti_393638" "A_sumichrasti_393640"
## [58] "A_woodhouseii_334062" "A_woodhouseii_334063" "A_woodhouseii_334086"
## [61] "A_woodhouseii_334088" "A_woodhouseii_334096" "A_woodhouseii_334132"
## [64] "A_woodhouseii_334133" "A_woodhouseii_334134" "A_woodhouseii_334142"
## [67] "A_woodhouseii_334153" "A_woodhouseii_334161" "A_woodhouseii_334188"
## [70] "A_woodhouseii_334190" "A_woodhouseii_334196" "A_woodhouseii_334210"
## [73] "A_woodhouseii_334211" "A_woodhouseii_334217" "A_woodhouseii_334230"
## [76] "A_woodhouseii_334241" "A_woodhouseii_334242" "A_woodhouseii_334243"
## [79] "A_woodhouseii_334244" "A_woodhouseii_334247" "A_woodhouseii_334250"
## [82] "A_woodhouseii_343458" "A_woodhouseii_343461" "A_woodhouseii_343476"
## [85] "A_woodhouseii_343480" "A_woodhouseii_343481" "A_woodhouseii_343483"
## [88] "A_woodhouseii_343497" "A_woodhouseii_393605" "A_woodhouseii_393606"
## [91] "A_woodhouseii_393697" "A_woodhouseii_393698" "A_woodhouseii_393702"
## [94] "A_woodhouseii_393712" "A_woodhouseii_393713" "A_woodhouseii_395768"
#retain only ((island,cali),int.wood,fla) as ((P1,P2),P3,P4)
vcfR@gt<-vcfR@gt[,c(1:47,58:74)]
dim(vcfR)
## variants fix_cols gt_cols
## 16307 8 64
#retain only SNPs that still have a maf > 0
source("~/Desktop/snpfiltR/min_mac.R")
vcfR<-min_mac(vcfR, min.mac = 1)
## [1] "21.28% of SNPs fell below a minor allele count of 1 and were removed from the VCF"
dim(vcfR)
## variants fix_cols gt_cols
## 12837 8 64
write.vcf(vcfR, file="~/Desktop/aph.data/dfs/dfs.isl.cali.wood.fla.vcf.gz")
#generate popmap
colnames(vcfR@gt)
## [1] "FORMAT" "A_californica_333849" "A_californica_333854"
## [4] "A_californica_333855" "A_californica_333857" "A_californica_333860"
## [7] "A_californica_333862" "A_californica_333871" "A_californica_333873"
## [10] "A_californica_333874" "A_californica_333914" "A_californica_333917"
## [13] "A_californica_333931" "A_californica_333932" "A_californica_333934"
## [16] "A_californica_334002" "A_californica_334012" "A_californica_334015"
## [19] "A_californica_334017" "A_californica_334019" "A_woodhouseii_334171"
## [22] "A_californica_342037" "A_californica_342048" "A_californica_342050"
## [25] "A_californica_342066" "A_californica_342072" "A_californica_343428"
## [28] "A_californica_343430" "A_californica_343432" "A_californica_343438"
## [31] "A_californica_343451" "A_californica_393615" "A_californica_393721"
## [34] "A_coerulescens_396251" "A_coerulescens_396254" "A_coerulescens_396256"
## [37] "A_coerulescens_396263" "A_coerulescens_396264" "A_coerulescens_396265"
## [40] "A_insularis_334025" "A_insularis_334029" "A_insularis_334031"
## [43] "A_insularis_334032" "A_insularis_334033" "A_insularis_334034"
## [46] "A_insularis_334037" "A_insularis_334038" "A_woodhouseii_334062"
## [49] "A_woodhouseii_334063" "A_woodhouseii_334086" "A_woodhouseii_334088"
## [52] "A_woodhouseii_334096" "A_woodhouseii_334132" "A_woodhouseii_334133"
## [55] "A_woodhouseii_334134" "A_woodhouseii_334142" "A_woodhouseii_334153"
## [58] "A_woodhouseii_334161" "A_woodhouseii_334188" "A_woodhouseii_334190"
## [61] "A_woodhouseii_334196" "A_woodhouseii_334210" "A_woodhouseii_334211"
## [64] "A_woodhouseii_334217"
names<-colnames(vcfR@gt)[-1]
pops<-c(rep("P2", times=19),"P3",rep("P2", times=12),rep("OG", times=6),rep("P1",times=8),rep("P3",times=17))
d<-data.frame(names=names,pops=pops)
table(d$pops)
##
## OG P1 P2 P3
## 6 8 31 18
#write popmap
write.table(d, file="~/Desktop/aph.data/dfs/pops.isl.cali.wood.fla.txt", sep = "\t", row.names = F, col.names = F, quote = F)
#run the following in bash to generate your 3D sfs for the 3 focal species (plus outgroup): ((sumi,int.wood),cali,fla) as ((P1,P2),P3,P4)
#convert vcf to .geno
#cd ~/Desktop/aph.data/dfs
#python /Users/devder/Downloads/genomics_general-0.4/VCF_processing/parseVCF.py -i dfs.sumi.wood.cali.fla.vcf.gz -o sumi.wood.cali.fla.geno.gz
#compute basecounts for all sites, popsfile must identify each pop
#python /Users/devder/Downloads/genomics_general-0.4/freq.py -g sumi.wood.cali.fla.geno.gz -p P1 -p P2 -p P3 -p OG --popsFile pops.sumi.wood.cali.fla.txt | gzip > sumi.wood.cali.fla.basecounts.tsv.gz
#now run sfs.py to generate sfs. must have identical number of samples in p1 and p2
#python /Users/devder/Downloads/genomics_general-0.4/sfs.py -i sumi.wood.cali.fla.basecounts.tsv.gz --inputType baseCounts --outgroup OG --FSpops P1 P2 P3 --subsample 20 20 20 --pref mydata. --suff sumi.wood.cali.fla.sfs
#run the following in bash to generate your 3D sfs for the 3 focal species (plus outgroup): ((island,cali),int.wood,fla) as ((P1,P2),P3,P4)
#convert vcf to .geno
#cd ~/Desktop/aph.data/dfs
#python /Users/devder/Downloads/genomics_general-0.4/VCF_processing/parseVCF.py -i dfs.isl.cali.wood.fla.vcf.gz -o isl.cali.wood.fla.geno.gz
#compute basecounts for all sites, popsfile must identify each pop
#python /Users/devder/Downloads/genomics_general-0.4/freq.py -g isl.cali.wood.fla.geno.gz -p P1 -p P2 -p P3 -p OG --popsFile pops.isl.cali.wood.fla.txt | gzip > isl.cali.wood.fla.basecounts.tsv.gz
#now run sfs.py to generate sfs. must have identical number of samples in p1 and p2
#python /Users/devder/Downloads/genomics_general-0.4/sfs.py -i isl.cali.wood.fla.basecounts.tsv.gz --inputType baseCounts --outgroup OG --FSpops P1 P2 P3 --subsample 16 16 16 --pref mydata. --suff isl.cali.wood.fla.sfs
#####
####
###
##
#
# Simon H. Martin 2020
# simon.martin@ed.ac.uk
################################ Start Here ####################################
# This script accompanies the paper:
# "Signatures of introgression across the allelel frequency spectrum"
# by Simon H. Martin and William Amos
# Each section below computes and plots the D frequency spectrum (DFS)
# from a different empirical site frequency spectrum (SFS).
# The input SFS is provided in tabular format:
# A tab-delimited table in which the first three columns (or 4 in the case of a 4D SFS)
# give the allele count in each population (equivalent to the indices of the multidimensional SFS).
# The final column gives the corresponding number of sites.
# In most cases the input SFS is 3D, with the first three columns corresponding to
# populations P1, P2 and P3 (see the paper for details). This means that the SFS is
# polarized, and the outgroup is assumed to carry the ancestral allele at these sites.
# In the case of Helcionius butterflies, we provide an example of using a 4D SFS.
# In this case, the input SFS is not polarzed. This means that the frequencies provided
# correspond to minor allele counts, and the fourth column gives the count in teh outgroup.
# Before DFS can be computed from a 4D SFS, it must first be polarized, and sites at
# which the outgroup is polymorpic ust be discarded.
#The functions to compute DFS and related statistics are provided in the accompanything script DFS.R
#First import these functions
source("~/Downloads/DFS.R")
################################################################################
############################## Arabidopsis ####################################
################################################################################
### import the frequency spectrum
FS <- read.table("/Users/devder/Desktop/aph.data/dfs/mydata.P1_P2_P3sumi.wood.cali.fla.sfs")
### get Dfs
dfs_data <- get.DFS(base_counts=FS[,-4], #base counts are the first three columns (i.e everything minus column 4)
site_counts=FS[,4], # site counts are column 4
Ns = c(20,20,20)) # Ns provide the haploid sample sizes of each population (1 and 2 must always be equal)
### plot
#png("/Users/devder/Desktop/aph.data/dfs/DFS_isl_ca_intwood_fla.png", width = 2000, height = 1000, res=300)
par(mar=c(1,4,1,1))
plotDFS(dfs_data$DFS, dfs_data$weights, method="lines", col_D="red", no_xlab=T)
#dev.off()
# We can also compute related statistics
D <- get.D.from.base.counts(base_counts=FS[,-4], site_counts=FS[,4], Ns = c(100,100,2)) #overall D
f <- get.f.from.base.counts(base_counts=FS[,-4], site_counts=FS[,4], Ns = c(100,100,2)) # overall fraction of introgression
dcfs <- get.dcfs(base_counts=FS[,-4], site_counts=FS[,4], Ns = c(100,100,2)) # dcfs from Yang et al. 2012
### code for exporting a table of plotted values
# write.table(data.frame(D=round(dfs_data$DFS,4),
# weight=round(dfs_data$weights,4)),
# file="empirical_data/Arabidopsis/DFS_arabidopsis.lyr2_lyr4_arn4.csv",
# sep=",",row.names=FALSE,quote=FALSE)
#next orientation
### import the frequency spectrum
FS <- read.table("/Users/devder/Desktop/aph.data/dfs/mydata.P1_P2_P3isl.cali.wood.fla.sfs")
### get Dfs
dfs_data <- get.DFS(base_counts=FS[,-4], #base counts are the first three columns (i.e everything minus column 4)
site_counts=FS[,4], # site counts are column 4
Ns = c(16,16,16)) # Ns provide the haploid sample sizes of each population (1 and 2 must always be equal)
### plot
#png("/Users/devder/Desktop/aph.data/dfs/DFS_isl_ca_intwood_fla.png", width = 2000, height = 1000, res=300)
par(mar=c(1,4,1,1))
plotDFS(dfs_data$DFS, dfs_data$weights, method="lines", col_D="red", no_xlab=T)
#dev.off()
# We can also compute related statistics
D <- get.D.from.base.counts(base_counts=FS[,-4], site_counts=FS[,4], Ns = c(100,100,2)) #overall D
f <- get.f.from.base.counts(base_counts=FS[,-4], site_counts=FS[,4], Ns = c(100,100,2)) # overall fraction of introgression
dcfs <- get.dcfs(base_counts=FS[,-4], site_counts=FS[,4], Ns = c(100,100,2)) # dcfs from Yang et al. 2012