I am trying to subset SNPs from 32 cultivars. I have the following exclusion criteria:

  1. >20% missing data
  2. SNPs present in less than 50% of genotypes
  3. Each genotype must have at least 50% of included SNPs
  4. >=60% heterozygous calls
  5. Bi-allelic SNPs only
  6. No very rare variants (i.e. MAF<0.05)

    bcftools view --output-type u --min-alleles 2 --max-alleles 2 \ 
      --types snps --exclude MAF[0]<0.05 --exclude ** \
        --known variants_raw_sorted.bcf > SNP_filtered_sorted_kkf.bcf

Where the ** is, I am trying to implement the remainder my exclusion criteria, (criteria 1-5), however I am having a great deal of difficulty understanding how to use the expression in bcftools. I am looking to understand how to use the information in the info fields (AN,AC, AF/MAF, NS) to implement my exclusion criteria. Any help is very much appreciated.

Further information : This is the upstream code that led to the bcf file.

fastp -i ACDC_R1.fastq.gz -I ACDC_R2.fastq.gz \ 
  -o fastp_ACDC_R1P.fastq.gz fastp_ACDC_R2P.fastq.gz \ 
  -R "ACDC_fastp_QC" -h "Fastp_ACDC" -j "ACDC_fastp" \ 
  -e 20 -c -a "auto" -L -r \ 
  -m --merged_out "ACDC_fastp_Merged" --out1 "ACDC_R1_fastp_Unmerged" \
  --out2 "ACDC_R2_fastp_Unmerged" \ 
  --unpaired1 "ACDC_R1Pass_R2Fail" --unpaired2 "ACDC_R2Pass_R1Fail"

bowtie2 -p 20 -x cs10 -U ACDC_fastp_Merged.fastq.gz -S ACDC_aligned.sam

samtools view -bS ACDC.sam > ACDC.bam

bcftools mpileup --output-type u -f cs10.fasta -b bam_list.txt | 
    bcftools call -vmO z -o all_raw_condition_1.vcf.gz --threads=20

bcftools view --output-type u --output-file variants_raw.bcf variants_raw.vcf.gz

bcftools sort --max-mem 14000 --output-type u --output-file variants_raw_sorted.bcf
  • 1
    $\begingroup$ Your criteria 2, 3 and 5 don't seem to make sense to me - could you be a bit more clear about the terminology, and what you mean by SNP in this case? Generally, 'SNP' is another word used for genetic marker, or position, so it doesn't make sense for it to be present in 50% of genotypes. $\endgroup$
    – user438383
    Apr 6 '20 at 10:44
  • $\begingroup$ Here I mean SNP in the usual sense - i.e. a mutation of a single nucleotide at a set position. Essentially, the idea is to use SNPs to discriminate between different cultivars. E.g. with a perfect set of 5 SNPs, you could discriminate between 2^5 cultivars. To clarify the criteria, I will remove criteria 5 - misleading. For the others : Individually, we would want each SNP (i.e. mutation) in our set to be present in less than 50% of genotypes (criteria 2). Across our whole set of SNPs, we would want at least 50% of the included SNPs to be present in each genotype (crit 3). $\endgroup$
    – Alex Cull
    Apr 6 '20 at 18:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.