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

<ol>
<li>>20% missing data</li>
<li>SNPs present in less than 50% of genotypes</li>
<li>Each genotype must have at least 50% of included SNPs</li>
<li>>=60% heterozygous calls</li>
<li>Bi-allelic SNPs only</li>
<li><p>No very rare variants (i.e. MAF&lt;0.05)</p>

<pre><code>bcftools view --output-type u --min-alleles 2 --max-alleles 2 \ 
  --types snps --exclude MAF[0]&lt;0.05 --exclude ** \
    --known variants_raw_sorted.bcf &gt; SNP_filtered_sorted_kkf.bcf
</code></pre></li>
</ol>

<p>Where the <code>**</code> 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.</p>

<hr>

<p>Further information : This is the upstream code that led to the bcf file.</p>

<pre><code>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 &gt; 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
</code></pre>

CC BY-SA 4.0
2020-04-06T10:44:31.783Z
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.
user438383
CC BY-SA 4.0
2020-04-06T18:32:53.653Z
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).
Alex Cull
<snp><bcftools><subset>
Subsetting SNPs with specific exclusion criteria using bcftools
# Subsetting SNPs with specific exclusion criteria using bcftools

(Asked by: username_0 on 2020-04-05)

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

>20% missing data
SNPs present in less than 50% of genotypes
Each genotype must have at least 50% of included SNPs
>=60% heterozygous calls
Bi-allelic SNPs only
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

**Tags:** snp, bcftools, subset
> **Comments:**
> **username_1**: 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.  
> **username_0**: 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).  

---

## 0 answers

