SNP FILTERING
- Ischnura
- Raw data stats
- Step 1: Remove genotypes based on sequencing depth:
- Step 2: On each previous dataset we will filter out SNPs based on missing data
- Step 3: On each previous dataset we will filter samples based on missing data
- Step 4: Finally applying minimum minor allele count
- Running SNP counting script
Run this script while knitting it in Rmarkdown! This will create a report with the quality of each filtering protocol, that can be used as our output results file to choose the best filtering scheme. While this script is running it will jump from R to BASH. Be sure to have the needed software installed, mainly: src directory with auxiliary SNPs, vcftools and plink. This script was adapted to be used on Ischnura using as base O’Leary et al. (2018). These aren’t the loci you’e looking for: Principles of effective SNP filtering for molecular ecologists. Molecular Ecology. (Repository)[https://github.com/sjoleary/SNPFILT].
Ischnura
In this script we will filter first genotypes by sequencing depth, then loci by missing data, samples by missing data and apply a minor allele count in all finished filters. Add to the followig chunk the filters you wish to apply!
write(depths <- c("5","10","15"), "depths.tmp", ncolumns = 1)
write(loci <- c("90", "95"), "loci.tmp", ncolumns = 1)
write(samples <- c("10","15"), "samples.tmp", ncolumns = 1)
write(mac <- c("3"), "mac.tmp", ncolumns = 1)Raw data stats
Query stats
Query raw stats using vcftools. Filtering for bi-allelic snps only
# Testing if duplicated SNPs were called in the same position (this was a problem before using stacks populations --ordered-export flag). See for more info: https://groups.google.com/g/stacks-users/c/Ag8YyEFe7z0
vcftools --vcf ../data/04_Populations_Filter/populations.snps.vcf --hardy --out ../data/04_Populations_Filter/raw
# Estimating summary statistics of quality
out=../data/05_FilteringTests/00_raw_data
mkdir -p $out
vcftools --vcf ../data/04_Populations_Filter/populations.snps.vcf --out $out/RawData --min-alleles 2 --max-alleles 2 --remove-indels --recode --recode-INFO-all
for i in het depth site-mean-depth missing-indv missing-site plink
do
vcftools --vcf $out/RawData.recode.vcf --out $out/RawData --$i
done
# Estimating heterozygosity
plink --file $out/RawData --hardy --out $out/RawDataVisualize
## Is the database free from duplicated SNPs?## [1] TRUE
## [1] "Data set contains 305 individuals and 2177571 loci"Step 1: Remove genotypes based on sequencing depth:
** While keeping again biallelic snps only***
# Removing all not vcf files from raw data
find ../data/05_FilteringTests/00_raw_data/ -type f ! -name '*.vcf' -delete
# Filtering by sequencing depth
dir=../data/05_FilteringTests
for d in $(cat depths.tmp)
do
out=$dir/DP$d
mkdir -p $out
vcftools --vcf $dir/00_raw_data/RawData.recode.vcf --out $out/DP$d --minDP $d --recode --recode-INFO-all
# Summarizing
for i in het depth site-mean-depth missing-indv missing-site plink
do
vcftools --vcf $out/DP$d.recode.vcf --out $out/DP$d --$i
done
# Estimating heterozygosity
plink --file $out/DP$d --hardy --out $out/DP$d
donePlot stats for filtered data set
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP5"
## [1] "Data set contains 305 individuals and 2177571 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP10"
## [1] "Data set contains 305 individuals and 2177571 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP15"
## [1] "Data set contains 305 individuals and 2177571 loci"Step 2: On each previous dataset we will filter out SNPs based on missing data
dir=../data/05_FilteringTests
for d in $(cat depths.tmp)
do
# Removing all not vcf files from previous filters
find $dir/DP$d/ -type f ! -name '*.vcf' -delete
# Filtering by loci missing data
for l in $(cat loci.tmp)
do
out=$dir/DP${d}geno${l}
mkdir -p $out
vcftools --vcf $dir/DP$d/DP$d.recode.vcf --out $out/DP${d}geno${l} --max-missing .$l --recode --recode-INFO-all
# Summarizing
for i in het depth site-mean-depth missing-indv missing-site plink
do
vcftools --vcf $out/DP${d}geno${l}.recode.vcf --out $out/DP${d}geno${l} --$i
done
# Estimating heterozygosity
plink --file $out/DP${d}geno${l} --hardy --out $out/DP${d}geno${l}
done
donePlot stats for filtered data set
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP5geno90"
## [1] "Data set contains 305 individuals and 10093 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP5geno95"
## [1] "Data set contains 305 individuals and 133 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP10geno90"
## [1] "Data set contains 305 individuals and 2524 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP10geno95"
## [1] "Data set contains 305 individuals and 6 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP15geno90"
## [1] "Data set contains 305 individuals and 728 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP15geno95"
## [1] "Data set contains 305 individuals and 4 loci"Step 3: On each previous dataset we will filter samples based on missing data
for(d in depths){
for(l in loci){
dir <- paste0("../data/05_FilteringTests/DP",d,"geno",l)
vcf <- paste0("DP",d,"geno",l)
imissing <- read.table(paste0(dir,"/",vcf,".imiss"), header = TRUE, stringsAsFactors = FALSE)
for(s in samples){
LQ_indv <- imissing %>%
filter(F_MISS > as.numeric(paste0("0.",s))) %>%
select(INDV)
p <- ggplot(imissing, aes(x = F_MISS)) +
geom_histogram(binwidth = .01, color = "black", fill = "grey95") +
geom_vline(aes(xintercept = mean(F_MISS, na.rm = TRUE)),
color = "red", linetype = "dashed", size = 1) +
geom_vline(aes(xintercept = as.numeric(paste0("0.",s))),
color = "darkblue", linetype = "dashed", size = 1) +
labs(title = vcf, x = "missing data per indv") +
scale_x_continuous(limits = c(0, 1)) +
theme_standard
print(p)
write.table(LQ_indv, paste0(dir,"/IND_FIL_",s,".txt"), col.names = FALSE, row.names = FALSE, quote = FALSE)
}
}
}
Applying samples filtering
dir=../data/05_FilteringTests
for d in $(cat depths.tmp)
do
for l in $(cat loci.tmp)
do
# Filtering samples by missing data
for s in $(cat samples.tmp)
do
out=$dir/DP${d}geno${l}IND${s}
mkdir -p $out
vcftools --vcf $dir/DP${d}geno${l}/DP${d}geno${l}.recode.vcf --out $out/DP${d}geno${l}IND${s} --remove $dir/DP${d}geno${l}/IND_FIL_$s.txt --recode --recode-INFO-all
# Summarizing
for i in het depth site-mean-depth missing-indv missing-site plink
do
vcftools --vcf $out/DP${d}geno${l}IND${s}.recode.vcf --out $out/DP${d}geno${l}IND${s} --$i
done
# Estimating heterozygosity
plink --file $out/DP${d}geno${l}IND${s} --hardy --out $out/DP${d}geno${l}IND${s}
done
# Removing all not vcf files from previous filters
find $dir/DP${d}geno${l}/ -type f ! -name '*.vcf' -delete
done
donePlot stats for filtered data set
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP5geno90IND10"
## [1] "Data set contains 259 individuals and 10093 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP5geno90IND15"
## [1] "Data set contains 270 individuals and 10093 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP5geno95IND10"
## [1] "Data set contains 276 individuals and 133 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP5geno95IND15"
## [1] "Data set contains 280 individuals and 133 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP10geno90IND10"
## [1] "Data set contains 255 individuals and 2524 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP10geno90IND15"
## [1] "Data set contains 267 individuals and 2524 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP10geno95IND10"
## [1] "Data set contains 266 individuals and 6 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP10geno95IND15"
## [1] "Data set contains 266 individuals and 6 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP15geno90IND10"
## [1] "Data set contains 248 individuals and 728 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP15geno90IND15"
## [1] "Data set contains 259 individuals and 728 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP15geno95IND10"
## [1] "Data set contains 278 individuals and 4 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP15geno95IND15"
## [1] "Data set contains 278 individuals and 4 loci"Step 4: Finally applying minimum minor allele count
dir=../data/05_FilteringTests
for d in $(cat depths.tmp)
do
for l in $(cat loci.tmp)
do
for s in $(cat samples.tmp)
do
# Removing all not vcf files from previous filters
find $dir/DP${d}geno${l}IND${s}/ -type f ! -name '*.vcf' -delete
for m in $(cat mac.tmp)
do
out=$dir/DP${d}geno${l}IND${s}mac${m}
mkdir -p $out
vcftools --vcf $dir/DP${d}geno${l}IND${s}/DP${d}geno${l}IND${s}.recode.vcf --out $out/DP${d}geno${l}IND${s}mac${m} --mac $m --recode --recode-INFO-all
# Summarizing
for i in het depth site-mean-depth missing-indv missing-site plink
do
vcftools --vcf $out/DP${d}geno${l}IND${s}mac${m}.recode.vcf --out $out/DP${d}geno${l}IND${s}mac${m} --$i
done
# Estimating heterozygosity
plink --file $out/DP${d}geno${l}IND${s}mac${m} --hardy --out $out/DP${d}geno${l}IND${s}mac${m}
done
done
done
donePlot stats for filtered data set
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP5geno90IND10mac3"
## [1] "Data set contains 259 individuals and 10093 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP5geno90IND15mac3"
## [1] "Data set contains 270 individuals and 10093 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP5geno95IND10mac3"
## [1] "Data set contains 276 individuals and 133 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP5geno95IND15mac3"
## [1] "Data set contains 280 individuals and 133 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP10geno90IND10mac3"
## [1] "Data set contains 255 individuals and 2524 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP10geno90IND15mac3"
## [1] "Data set contains 267 individuals and 2524 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP10geno95IND10mac3"
## [1] "Data set contains 266 individuals and 6 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP10geno95IND15mac3"
## [1] "Data set contains 266 individuals and 6 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP15geno90IND10mac3"
## [1] "Data set contains 248 individuals and 728 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP15geno90IND15mac3"
## [1] "Data set contains 259 individuals and 728 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP15geno95IND10mac3"
## [1] "Data set contains 278 individuals and 4 loci"
## [1] ""
## [1] ""
## [1] "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
## [1] "DP15geno95IND15mac3"
## [1] "Data set contains 278 individuals and 4 loci"Running SNP counting script
Comparing SNPs/contigs/indv between filtering schemes.
## FILTER SNP CONTIG INDV
## 1 00_raw_data/RawData 2177571 65 305
## 2 DP10/DP10 2177571 65 305
## 3 DP10geno90/DP10geno90 2524 16 305
## 4 DP10geno90IND10/DP10geno90IND10 2524 16 255
## 5 DP10geno90IND10mac3/DP10geno90IND10mac3 2524 16 255
## 6 DP10geno90IND15/DP10geno90IND15 2524 16 267
## 7 DP10geno90IND15mac3/DP10geno90IND15mac3 2524 16 267
## 8 DP10geno95/DP10geno95 6 3 305
## 9 DP10geno95IND10/DP10geno95IND10 6 3 266
## 10 DP10geno95IND10mac3/DP10geno95IND10mac3 6 3 266
## 11 DP10geno95IND15/DP10geno95IND15 6 3 266
## 12 DP10geno95IND15mac3/DP10geno95IND15mac3 6 3 266
## 13 DP15/DP15 2177571 65 305
## 14 DP15geno90/DP15geno90 728 14 305
## 15 DP15geno90IND10/DP15geno90IND10 728 14 248
## 16 DP15geno90IND10mac3/DP15geno90IND10mac3 728 14 248
## 17 DP15geno90IND15/DP15geno90IND15 728 14 259
## 18 DP15geno90IND15mac3/DP15geno90IND15mac3 728 14 259
## 19 DP15geno95/DP15geno95 4 1 305
## 20 DP15geno95IND10/DP15geno95IND10 4 1 278
## 21 DP15geno95IND10mac3/DP15geno95IND10mac3 4 1 278
## 22 DP15geno95IND15/DP15geno95IND15 4 1 278
## 23 DP15geno95IND15mac3/DP15geno95IND15mac3 4 1 278
## 24 DP5/DP5 2177571 65 305
## 25 DP5geno90/DP5geno90 10093 17 305
## 26 DP5geno90IND10/DP5geno90IND10 10093 17 259
## 27 DP5geno90IND10mac3/DP5geno90IND10mac3 10093 17 259
## 28 DP5geno90IND15/DP5geno90IND15 10093 17 270
## 29 DP5geno90IND15mac3/DP5geno90IND15mac3 10093 17 270
## 30 DP5geno95/DP5geno95 133 13 305
## 31 DP5geno95IND10/DP5geno95IND10 133 13 276
## 32 DP5geno95IND10mac3/DP5geno95IND10mac3 133 13 276
## 33 DP5geno95IND15/DP5geno95IND15 133 13 280
## 34 DP5geno95IND15mac3/DP5geno95IND15mac3 133 13 280