I have four vcf files, SNPs_s1.vcf, SNPs_s2.vcf, SNPs_s3.vcf, and SNPs_s4.vcf, which contain information about SNPs. These vcf files were obtained by using the following methods: the initial input files were short-paired reads

  1. I did mapping with minimap2

    ./minimap2 -ax sr ref.fa read1.fq.gz read2.fq.gz > aln.sam
  2. converted to bam file

    samtools view -Sb aln.sam > aln.bam
  3. Sorted the bam file

    samtools sort aln.bam > aln_sorted.bam
  4. created an index file

    samtools index aln_sorted.bam
  5. created an index file of reference

    samtools faidx ref.fa
  6. extracted the raw variants using freebayes

    freebayes -f ref.fa aln_sorted.bam > out_var_s1.vcf
  7. extracted SNPs using vcf tools

    vcftools --vcf out_var_s1.vcf --remove-indels --recode \
       --recode-INFO-all --out SNPs_s1.vcf

I would like to calculate the fixation index (Fst) between SNPs_s1.vcf and SNPs_s2.vcf, SNPs_s1.vcf and SNPs_s3.vcf etc. I attempted to calculate Fst using bcftools, but I have to provide the id list. I tried the code below for that.

bcftools query -l snps_s1.vcf

The output was "unknown".

How do I calculate Fst using the above vcf files? Your suggestions would be appreciated.

  • $\begingroup$ Can you please explain a bit more about Fst, and why you need to calculate it? I'm used to Fst being interpreted as a population statistic, and it doesn't [currently] make sense to me to calculate it for only a few individuals. $\endgroup$
    – gringer
    Dec 9, 2023 at 9:29
  • $\begingroup$ @gringer I would like to use Fst to find the genetic diversity between populations. Which method is better for few individuals? $\endgroup$
    – hina
    Dec 9, 2023 at 23:44
  • $\begingroup$ Who is asking you to calculate Fst, and why? If you're doing population genetics, you need populations. A few individuals do not make a population, unless it's kākāpō or kakī, or some other similarly endangered species. $\endgroup$
    – gringer
    Dec 10, 2023 at 6:11
  • 1
    $\begingroup$ The first issue is easy, the name of your sample in your snps_s1.vcf is set to unknown. You can fix that bit with sed -i '/^#CHROM/s/unknown/sample1/' snps_s1.vcf. But gringer's points are the main issue. What are you really trying to do here with only 4 individuals? $\endgroup$
    – terdon
    Dec 10, 2023 at 14:15


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.