AI Analysis ‐ Pipeline Audit - rsalani13/EPP622-2026-SNPCalling GitHub Wiki
Evaluation of the pipeline generated at the beginning of the learning module.
Model used: OpenAI - ChatGPT, free version
Date generated: September 14, 2026
Exact prompt:
write me a a complete short-read SNP calling pipeline for a haploid fungal genome, starting from paired FASTQ and ending at a filtered VCF
Evaluation
Three things it got right.
The first thing the AI pipeline did right was create an organized directory structure before beginning. While I did this as I went through the steps of my pipeline, it may have been helpful for me to create directories beforehand to familiarize myself with my workspace. Another thing AI got right was filtering by a variety of characteristics. It filtered for variant quality, read depth, mapping quality, allele balance, missingness, and MAF. While too stringent of filtering may lead to false negatives, gentle filtering for many trait will likely help it achieve accurate SNP calling. The AI model also performed quality control at frequently throughout the pipeline. This is important to identify errors as they come up so they don't get lost in the pipeline.
Three things you don't understand and would want to learn more about before using. Or do you understand all of it? If so be ready to answer questions about it.
One thing that I don't understand in the AI pipeline is the inclusion of a section to account for ploidy. While I understand ploidy and realize that this may make it necessary for data to be handled differently, I do not know how this changes the code or what the exact issues would be if this was not accounted for. The AI pipeline also referenced allele balance. I do not know if this was included in our pipeline that we discussed in class and I do not understand how this would be used as a tool for filtering. A third thing I didn't understand in the Ai pipeline was "Check transition/transversion ratio". I don't know what this is at all, but I would like to learn more about it.
How specifically would you assess the results of this pipeline for accuracy?
In order to assess the results of this pipeline for accuracy, I would look at each Multiqc report generated as it moves through the pipeline. The graphs and charts offer insight into the data that can be compared to the previous Multiqc report to confirm that that particular step seemed to function correctly. For the final results, I would look at the stats generated in the final VCF statistics report. I'd also add the step in our pipeline that the AI did not include: performing both gentle and stringent filtering. I think this is a helpful way to see how your data behaves; odd results are clear indicators of something being wrong with the code.
What did the labs teach you that turned out to matter for evaluating the AI-generated code?
Before we completed module 2 labs in class, I tried to read the AI pipeline and got lost at the first step. The labs taught me the names of file types (like .fastq, .fasta, .bam, etc.) and tools (like BWA-MEM2, fastQC/MultiQC, SAM/SAM, etc.) that are commonly used for SNP calling. The AI pipeline used the same tools that we learned how to work with in class. I also learned the steps necessary before the actual calling of SNPs. I have experience in extracting DNA for SNP calling and have worked with SNPs, but was never involved with any of the steps in between. Learning the steps for the bioinformatics side of this process was necessary for me to know enough to evaluate AI code.
What is the value of learning about a computational genomics task before using AI for it — or is there none?
It's very valuable to learn about a computational genomics task before using AI for it for many reasons. Being able to troubleshoot your own data is a vital part of being a researcher; this includes both issues caused by the AI and issues caused by working with imperfect, real-life data. AI models have been known to make mistakes and, as the person responsible for the data and the communication of it, those mistakes reflect on the bioinformatician. Presenting information means the presenter must be able to defend it, and relying completely on AI to do the work for you will make it impossible to do so under scrutiny.
AI Response:
for reference
Absolutely. Below is a complete, practical short-read SNP-calling pipeline for a haploid fungal genome, starting with paired-end FASTQ files and ending with a filtered, analysis-ready VCF.
I’ll use a standard workflow:
FASTQ → QC → adapter/quality trimming → reference indexing → read alignment → BAM processing → mapping/duplicate QC → variant calling → SNP filtering → final VCF
This assumes:
paired-end Illumina reads
a haploid fungal genome
a reference genome in FASTA format
one or more samples
Linux/HPC environment
conda/environment modules are available
1. Pipeline overview
Sample_R1.fastq.gz ─┐
├─ FastQC
Sample_R2.fastq.gz ─┘
↓
MultiQC (optional)
↓
fastp / trimming
↓
trimmed R1/R2
↓
BWA-MEM2
↓
SAM → sorted BAM
↓
Mark duplicates
↓
BAM indexing + QC
↓
bcftools
variant calling
↓
raw VCF
↓
SNP-only filtering
↓
quality/depth filtering
↓
filtered SNP VCF
For a fungal genome, I would generally use BWA-MEM2 + bcftools for a straightforward short-read pipeline.
2. Set up your directory structure
For example:
mkdir -p fungal_snp_pipeline
cd fungal_snp_pipeline
mkdir -p \
reference \
raw_fastq \
trimmed \
qc \
alignments \
variants \
logs
Your directory would look like:
fungal_snp_pipeline/
├── reference/
├── raw_fastq/
├── trimmed/
├── qc/
├── alignments/
├── variants/
└── logs/
Put your reference genome here:
reference/fungus_reference.fasta
and your reads here:
raw_fastq/
├── sample1_R1.fastq.gz
├── sample1_R2.fastq.gz
├── sample2_R1.fastq.gz
├── sample2_R2.fastq.gz
...
3. Create a software environment
One option is:
conda create -n fungal_snp \
fastqc \
multiqc \
fastp \
bwa-mem2 \
samtools \
bcftools \
-c conda-forge -c bioconda
Activate it:
conda activate fungal_snp
Check that everything is available:
fastqc --version
multiqc --version
fastp --version
bwa-mem2 version
samtools --version
bcftools --version
4. Reference genome preparation
First make sure your reference is a proper FASTA.
grep "^>" reference/fungus_reference.fasta
You can get chromosome/contig lengths with:
samtools faidx reference/fungus_reference.fasta
This creates:
fungus_reference.fasta.fai
Create the BWA-MEM2 index:
bwa-mem2 index reference/fungus_reference.fasta
You should now have several index files in reference/.
5. Initial FASTQ quality control
Run FastQC on all raw FASTQ files:
fastqc \
raw_fastq/*.fastq.gz \
-o qc/
Then summarize the results with MultiQC:
multiqc qc/ -o qc/multiqc_raw
The important things to look at include:
per-base sequence quality
adapter contamination
sequence duplication
overrepresented sequences
GC content
read length
number of reads
For Illumina data, you generally want high-quality bases toward the beginning of the reads and to identify whether adapters are present.
6. Trim adapters and low-quality bases
I recommend fastp because it performs adapter detection and quality trimming while producing useful reports.
For one sample:
fastp \
-i raw_fastq/sample1_R1.fastq.gz \
-I raw_fastq/sample1_R2.fastq.gz \
-o trimmed/sample1_R1.trimmed.fastq.gz \
-O trimmed/sample1_R2.trimmed.fastq.gz \
--detect_adapter_for_pe \
--trim_poly_g \
--cut_front \
--cut_tail \
--cut_mean_quality 20 \
--length_required 50 \
--html qc/sample1.fastp.html \
--json qc/sample1.fastp.json \
--thread 8
The important parameters here are:
--detect_adapter_for_pe
detect adapters in paired-end data.
--cut_mean_quality 20
trim low-quality regions.
--length_required 50
discard reads shorter than 50 bp after trimming.
You may want to adjust these depending on your read length and sequencing quality.
7. Run trimming for all samples
Instead of manually repeating the command, you can use a loop.
for R1 in raw_fastq/*_R1.fastq.gz
do
sample=$(basename "$R1" _R1.fastq.gz)
R2="raw_fastq/${sample}_R2.fastq.gz"
fastp \
-i "$R1" \
-I "$R2" \
-o "trimmed/${sample}_R1.trimmed.fastq.gz" \
-O "trimmed/${sample}_R2.trimmed.fastq.gz" \
--detect_adapter_for_pe \
--trim_poly_g \
--cut_front \
--cut_tail \
--cut_mean_quality 20 \
--length_required 50 \
--html "qc/${sample}.fastp.html" \
--json "qc/${sample}.fastp.json" \
--thread 8
done
Then summarize:
multiqc qc/ -o qc/multiqc_trimmed
8. Align reads to the fungal reference
Now align each sample using BWA-MEM2.
For example:
bwa-mem2 mem \
-t 8 \
-R '@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA' \
reference/fungus_reference.fasta \
trimmed/sample1_R1.trimmed.fastq.gz \
trimmed/sample1_R2.trimmed.fastq.gz \
> alignments/sample1.sam
The read group is important.
In particular:
ID = sequencing/read-group ID
SM = sample name
PL = sequencing platform
The SM field should uniquely identify your biological sample.
9. Convert SAM → sorted BAM
You generally don't want to retain the huge SAM file after successful conversion.
samtools view \
-@ 8 \
-b \
alignments/sample1.sam \
| samtools sort \
-@ 8 \
-o alignments/sample1.sorted.bam
Index it:
samtools index alignments/sample1.sorted.bam
You now have:
sample1.sorted.bam
sample1.sorted.bam.bai
Remove the SAM if disk space matters:
rm alignments/sample1.sam
10. Add read groups correctly for all samples
You can combine alignment and sorting into one command:
bwa-mem2 mem \
-t 8 \
-R '@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA' \
reference/fungus_reference.fasta \
trimmed/sample1_R1.trimmed.fastq.gz \
trimmed/sample1_R2.trimmed.fastq.gz \
| samtools sort \
-@ 8 \
-o alignments/sample1.sorted.bam
Then:
samtools index alignments/sample1.sorted.bam
For multiple samples:
for R1 in trimmed/*_R1.trimmed.fastq.gz
do
sample=$(basename "$R1" _R1.trimmed.fastq.gz)
R2="trimmed/${sample}_R2.trimmed.fastq.gz"
echo "Aligning $sample"
bwa-mem2 mem \
-t 8 \
-R "@RG\tID:${sample}\tSM:${sample}\tPL:ILLUMINA" \
reference/fungus_reference.fasta \
"$R1" \
"$R2" \
| samtools sort \
-@ 8 \
-o "alignments/${sample}.sorted.bam"
samtools index "alignments/${sample}.sorted.bam"
done
11. Check alignment statistics
For each sample:
samtools flagstat alignments/sample1.sorted.bam
Also:
samtools stats alignments/sample1.sorted.bam > \
qc/sample1.samtools.stats.txt
You can look at mapping quality:
samtools view alignments/sample1.sorted.bam \
| awk '{print $5}' \
| sort -n \
| uniq -c
A particularly useful statistic is the percentage of reads mapped to the reference.
You should investigate samples with unexpectedly low mapping rates.
12. Mark PCR duplicates
PCR duplicates can create artificial evidence for variants.
Use samtools markdup.
First, the BAM needs to have appropriate read-coordinate information. A common workflow is:
samtools fixmate \
-m \
alignments/sample1.sorted.bam \
alignments/sample1.fixmate.bam
Sort again:
samtools sort \
-o alignments/sample1.positionsort.bam \
alignments/sample1.fixmate.bam
Then mark duplicates:
samtools markdup \
-r \
alignments/sample1.positionsort.bam \
alignments/sample1.dedup.bam
Index:
samtools index alignments/sample1.dedup.bam
The -r option removes duplicates.
If you want to retain duplicates but mark them:
samtools markdup \
alignments/sample1.positionsort.bam \
alignments/sample1.dedup.bam
For variant calling, either approach can be appropriate, but removing/marking duplicates should be based on the library type and experimental design. For standard Illumina genomic DNA libraries, duplicate removal is commonly used.
13. Calculate sequencing depth
For example:
samtools depth \
alignments/sample1.dedup.bam \
> qc/sample1.depth.txt
Average depth:
awk '{sum += $3; n++} END {print "Mean depth:", sum/n}' \
qc/sample1.depth.txt
You can also calculate the fraction of reference positions covered at least 10×:
awk '
{
total++;
if ($3 >= 10) covered++;
}
END {
print "Fraction >=10x:", covered/total
}' qc/sample1.depth.txt
This is important because variant filtering thresholds should depend partly on your actual sequencing depth.
For example, a sample with:
30× mean coverage
can reasonably be treated differently from one with:
3× mean coverage
14. Generate the BAM list
If you have multiple samples:
ls alignments/*.dedup.bam > variants/bam_list.txt
Check it:
cat variants/bam_list.txt
It might look like:
alignments/sample1.dedup.bam
alignments/sample2.dedup.bam
alignments/sample3.dedup.bam
alignments/sample4.dedup.bam
15. Call variants with bcftools
For a multi-sample haploid dataset, one convenient approach is:
bcftools mpileup \
-f reference/fungus_reference.fasta \
-q 20 \
-Q 20 \
-Ou \
-b variants/bam_list.txt \
| \
bcftools call \
-mv \
-Ov \
-o variants/raw.vcf
Here:
-f
specifies the reference.
-q 20
requires mapping quality ≥20.
-Q 20
requires base quality ≥20.
-m
uses the multiallelic caller.
-v
outputs variant sites.
16. Prefer compressed VCF/BCF for large datasets
For a real fungal genome, I'd generally use compressed output:
bcftools mpileup \
-f reference/fungus_reference.fasta \
-q 20 \
-Q 20 \
-Ou \
-b variants/bam_list.txt \
| \
bcftools call \
-mv \
-Oz \
-o variants/raw.vcf.gz
Index:
bcftools index \
-t \
variants/raw.vcf.gz
Now you have:
raw.vcf.gz
raw.vcf.gz.tbi
17. Because this is haploid, specify ploidy
This is an important part of your pipeline.
Your fungus is haploid, so variant calling should reflect haploid genotypes.
Create:
cat > variants/haploid.ploidy << 'EOF'
## This file specifies haploid ploidy
EOF
However, rather than relying on a generic ploidy file, you can explicitly tell bcftools call to use haploid ploidy:
bcftools call \
--ploidy 1 \
-mv \
-Oz \
-o variants/raw.vcf.gz
So the complete call becomes:
bcftools mpileup \
-f reference/fungus_reference.fasta \
-q 20 \
-Q 20 \
-Ou \
-b variants/bam_list.txt \
| \
bcftools call \
--ploidy 1 \
-mv \
-Oz \
-o variants/raw.vcf.gz
This is particularly important because a haploid organism should not be interpreted using diploid genotype probabilities.
18. Inspect the raw VCF
First look at the header:
bcftools view -h variants/raw.vcf.gz | head
Look at the first variants:
bcftools view variants/raw.vcf.gz | head -30
Count variants:
bcftools view \
-H \
variants/raw.vcf.gz \
| wc -l
19. Keep only SNPs
The raw VCF may contain:
SNPs
insertions
deletions
complex variants
If your goal is specifically SNPs:
bcftools view \
-v snps \
-Oz \
-o variants/snps.vcf.gz \
variants/raw.vcf.gz
Index:
bcftools index -t variants/snps.vcf.gz
20. Filter by variant quality
A basic first filter is:
bcftools filter \
-i 'QUAL>=30' \
-Oz \
-o variants/snps.Q30.vcf.gz \
variants/snps.vcf.gz
This retains variants with:
QUAL ≥ 30
21. Filter by read depth
Depth filtering is more complicated because you don't want to simply require enormous depth.
For a single sample, you might consider something like:
bcftools filter \
-i 'QUAL>=30 && DP>=10' \
-Oz \
-o variants/snps.filtered.vcf.gz \
variants/snps.vcf.gz
But for multiple samples, DP is generally the combined site depth.
You may instead want to filter individual genotypes using FORMAT/DP.
For example:
bcftools filter \
-S . \
-e 'FMT/DP<10' \
variants/snps.vcf.gz \
-Oz \
-o variants/snps.genotypeDP10.vcf.gz
This turns low-depth genotype calls into missing genotypes rather than automatically removing the entire SNP.
That distinction is important.
22. Filter mapping quality and allele balance
For a haploid fungal genome, I'd also pay attention to:
mapping quality
base quality
depth
genotype quality
allele balance
missingness
repetitive regions
For example, you can require genotype quality:
bcftools filter \
-S . \
-e 'FMT/GQ<20' \
variants/snps.genotypeDP10.vcf.gz \
-Oz \
-o variants/snps.DP10_GQ20.vcf.gz
23. Filter missingness
Suppose you have 20 fungal isolates and you only want variants present in at least 90% of samples.
You can use:
bcftools view \
-i 'F_MISSING<=0.10' \
variants/snps.DP10_GQ20.vcf.gz \
-Oz \
-o variants/snps.missing10.vcf.gz
This means:
≤10% missing genotypes
or approximately:
≥90% genotype completeness
24. Filter minor allele frequency
For population-genetic analyses, you may also want a MAF filter.
For example:
bcftools view \
-i 'MAF>=0.05' \
variants/snps.missing10.vcf.gz \
-Oz \
-o variants/snps.MAF05.vcf.gz
However, do not automatically use MAF ≥ 0.05.
If you're studying rare variants or looking for mutations unique to particular isolates, this would eliminate biologically interesting variants.
For GWAS/population analyses, MAF filtering is more commonly appropriate.
25. Final filtered VCF
A reasonable starting final filter for a typical high-coverage haploid fungal dataset might therefore be:
bcftools view \
-v snps \
-i 'QUAL>=30 && F_MISSING<=0.10' \
-Oz \
-o variants/final.filtered.snps.vcf.gz \
variants/raw.vcf.gz
Then index:
bcftools index \
-t \
variants/final.filtered.snps.vcf.gz
If you also want minimum genotype depth and GQ:
bcftools filter \
-S . \
-e 'FMT/DP<10 || FMT/GQ<20' \
variants/final.filtered.snps.vcf.gz \
-Oz \
-o variants/final.filtered.genotypes.vcf.gz
Index:
bcftools index \
-t \
variants/final.filtered.genotypes.vcf.gz
26. Check the final number of SNPs
bcftools view \
-H \
variants/final.filtered.genotypes.vcf.gz \
| wc -l
You can compare:
echo "Raw:"
bcftools view -H variants/raw.vcf.gz | wc -l
echo "SNPs:"
bcftools view -H variants/snps.vcf.gz | wc -l
echo "Final:"
bcftools view -H variants/final.filtered.genotypes.vcf.gz | wc -l
27. Generate a final VCF statistics report
This is extremely useful:
bcftools stats \
variants/final.filtered.genotypes.vcf.gz \
> variants/final.stats.txt
Then:
grep "^SN" variants/final.stats.txt
You'll get information such as:
number of SNPs
number of transitions
number of transversions
Ts/Tv ratio
number of samples
singleton counts
etc.
28. Check transition/transversion ratio
You can inspect:
grep "TSTV" variants/final.stats.txt
The Ts/Tv ratio can be useful as a rough diagnostic of variant quality.
29. Look at genotype calls
You can extract a simple genotype matrix:
bcftools query \
-f '%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n' \
variants/final.filtered.genotypes.vcf.gz \
> variants/genotype_matrix.txt
You'll get something resembling:
chr1 12345 A G 1 1 0 1
chr1 12891 C T 0 0 1 1
chr2 45213 G A 1 0 0 1
For haploid samples, the genotype representation should be:
0
1
.
rather than diploid:
0/0
0/1
1/1
depending on how the caller represents the haploid genotype.
30. Important: check that the organism is actually behaving as haploid
This matters biologically.
If you see many heterozygous calls such as:
0/1
in a supposedly haploid isolate, don't immediately assume the caller is wrong.
Potential explanations include:
contamination
mixed culture
aneuploidy
diploid/heterokaryotic regions
collapsed duplicated regions
repetitive sequences
paralogous mapping
sequencing artifacts
For a genuinely haploid isolate, widespread heterozygosity can be a useful QC signal.
31. Check contamination / mixed isolates
For each sample, calculate the proportion of heterozygous sites.
If the organism is strictly haploid, substantial heterozygosity is suspicious.
You can inspect genotype counts with:
bcftools query \
-f '[%GT\n]' \
variants/raw.vcf.gz
For more detailed analysis, I'd use a custom script to calculate per-sample:
total variant sites
called genotypes
missing genotypes
heterozygous genotypes
homozygous alternate genotypes
This is especially useful for fungal isolates because biological assumptions about ploidy can be complicated by heterokaryosis.
32. Mask repetitive/low-complexity regions
This is one of the biggest things I'd add to a real fungal SNP pipeline.
Variants in:
transposable elements
repetitive DNA
low-complexity regions
duplicated genes
poorly assembled regions
can have unreliable mapping.
If you have a BED file containing regions to exclude:
reference/masked_regions.bed
you can exclude them with:
bcftools view \
-T ^reference/masked_regions.bed \
variants/final.filtered.genotypes.vcf.gz \
-Oz \
-o variants/final.masked.filtered.vcf.gz
Then:
bcftools index \
-t \
variants/final.masked.filtered.vcf.gz
The ^ means exclude those regions.
33. A more complete final filtering strategy
For a high-quality haploid fungal population-genomics dataset, I would think about filtering in roughly this order:
Raw variants
↓
SNPs only
↓
QUAL ≥ 30
↓
minimum mapping/base quality during calling
↓
minimum genotype depth
↓
minimum genotype quality
↓
remove highly missing sites
↓
remove problematic genomic regions
↓
optional MAF filter
↓
FINAL SNP VCF
The exact numerical thresholds should be based on your sequencing depth and dataset rather than blindly copying values.
34. One complete bash script
If you want something you can actually put into a file and run, here is a compact version.
Save as:
fungal_snp_pipeline.sh
#!/bin/bash
set -euo pipefail
# ============================================================
# Haploid fungal short-read SNP calling pipeline
# FASTQ -> filtered SNP VCF
# ============================================================
# -----------------------------
# Settings
# -----------------------------
THREADS=8
REF="reference/fungus_reference.fasta"
RAW="raw_fastq"
TRIM="trimmed"
QC="qc"
BAM="alignments"
VAR="variants"
# -----------------------------
# Create directories
# -----------------------------
mkdir -p "$TRIM" "$QC" "$BAM" "$VAR" logs
# ============================================================
# 1. Reference indexing
# ============================================================
echo "Indexing reference..."
samtools faidx "$REF"
bwa-mem2 index "$REF"
# ============================================================
# 2. FASTQ QC
# ============================================================
echo "Running FastQC..."
fastqc \
"$RAW"/*.fastq.gz \
-o "$QC"
multiqc \
"$QC" \
-o "$QC/multiqc_raw"
# ============================================================
# 3. Adapter and quality trimming
# ============================================================
for R1 in "$RAW"/*_R1.fastq.gz
do
SAMPLE=$(basename "$R1" _R1.fastq.gz)
R2="$RAW/${SAMPLE}_R2.fastq.gz"
echo "Trimming $SAMPLE..."
fastp \
-i "$R1" \
-I "$R2" \
-o "$TRIM/${SAMPLE}_R1.trimmed.fastq.gz" \
-O "$TRIM/${SAMPLE}_R2.trimmed.fastq.gz" \
--detect_adapter_for_pe \
--trim_poly_g \
--cut_front \
--cut_tail \
--cut_mean_quality 20 \
--length_required 50 \
--html "$QC/${SAMPLE}.fastp.html" \
--json "$QC/${SAMPLE}.fastp.json" \
--thread "$THREADS"
done
# ============================================================
# 4. QC after trimming
# ============================================================
echo "Running MultiQC after trimming..."
multiqc \
"$QC" \
-o "$QC/multiqc_trimmed"
# ============================================================
# 5. Alignment
# ============================================================
for R1 in "$TRIM"/*_R1.trimmed.fastq.gz
do
SAMPLE=$(basename "$R1" _R1.trimmed.fastq.gz)
R2="$TRIM/${SAMPLE}_R2.trimmed.fastq.gz"
echo "Aligning $SAMPLE..."
bwa-mem2 mem \
-t "$THREADS" \
-R "@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tPL:ILLUMINA" \
"$REF" \
"$R1" \
"$R2" \
| samtools sort \
-@ "$THREADS" \
-o "$BAM/${SAMPLE}.sorted.bam"
samtools index \
"$BAM/${SAMPLE}.sorted.bam"
done
# ============================================================
# 6. Alignment QC
# ============================================================
for FILE in "$BAM"/*.sorted.bam
do
SAMPLE=$(basename "$FILE" .sorted.bam)
echo "QC for $SAMPLE..."
samtools flagstat \
"$FILE" \
> "$QC/${SAMPLE}.flagstat.txt"
samtools stats \
"$FILE" \
> "$QC/${SAMPLE}.stats.txt"
done
# ============================================================
# 7. Mark/remove duplicates
# ============================================================
for FILE in "$BAM"/*.sorted.bam
do
SAMPLE=$(basename "$FILE" .sorted.bam)
echo "Marking duplicates for $SAMPLE..."
samtools fixmate \
-m \
"$FILE" \
"$BAM/${SAMPLE}.fixmate.bam"
samtools sort \
-o "$BAM/${SAMPLE}.positionsort.bam" \
"$BAM/${SAMPLE}.fixmate.bam"
samtools markdup \
-r \
"$BAM/${SAMPLE}.positionsort.bam" \
"$BAM/${SAMPLE}.dedup.bam"
samtools index \
"$BAM/${SAMPLE}.dedup.bam"
done
# ============================================================
# 8. Generate BAM list
# ============================================================
ls "$BAM"/*.dedup.bam > "$VAR/bam_list.txt"
# ============================================================
# 9. Variant calling
# ============================================================
echo "Calling haploid variants..."
bcftools mpileup \
-f "$REF" \
-q 20 \
-Q 20 \
-Ou \
-b "$VAR/bam_list.txt" \
| bcftools call \
--ploidy 1 \
-mv \
-Oz \
-o "$VAR/raw.vcf.gz"
bcftools index \
-t \
"$VAR/raw.vcf.gz"
# ============================================================
# 10. Keep SNPs only
# ============================================================
echo "Keeping SNPs..."
bcftools view \
-v snps \
-Oz \
-o "$VAR/snps.vcf.gz" \
"$VAR/raw.vcf.gz"
bcftools index \
-t \
"$VAR/snps.vcf.gz"
# ============================================================
# 11. Basic site filtering
# ============================================================
echo "Filtering variants..."
bcftools view \
-i 'QUAL>=30 && F_MISSING<=0.10' \
-Oz \
-o "$VAR/snps.sitefiltered.vcf.gz" \
"$VAR/snps.vcf.gz"
bcftools index \
-t \
"$VAR/snps.sitefiltered.vcf.gz"
# ============================================================
# 12. Genotype filtering
# ============================================================
bcftools filter \
-S . \
-e 'FMT/DP<10 || FMT/GQ<20' \
"$VAR/snps.sitefiltered.vcf.gz" \
-Oz \
-o "$VAR/final.filtered.snps.vcf.gz"
bcftools index \
-t \
"$VAR/final.filtered.snps.vcf.gz"
# ============================================================
# 13. VCF statistics
# ============================================================
bcftools stats \
"$VAR/final.filtered.snps.vcf.gz" \
> "$VAR/final.filtered.stats.txt"
# ============================================================
# 14. Genotype matrix
# ============================================================
bcftools query \
-f '%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n' \
"$VAR/final.filtered.snps.vcf.gz" \
> "$VAR/final.genotype_matrix.txt"
# ============================================================
# Done
# ============================================================
echo ""
echo "======================================"
echo "SNP calling pipeline complete!"
echo "======================================"
echo ""
echo "Final VCF:"
echo "$VAR/final.filtered.snps.vcf.gz"
echo ""
echo "Number of final SNPs:"
bcftools view \
-H \
"$VAR/final.filtered.snps.vcf.gz" \
| wc -l
echo ""
echo "VCF statistics:"
echo "$VAR/final.filtered.stats.txt"
Run it with:
chmod +x fungal_snp_pipeline.sh
then:
./fungal_snp_pipeline.sh
35. What you ultimately get
The most important output is:
variants/final.filtered.snps.vcf.gz
with its index:
variants/final.filtered.snps.vcf.gz.tbi
Conceptually, the final VCF looks like:
CHROM POS REF ALT SAMPLE1 SAMPLE2 SAMPLE3
chr1 10543 A G 0 1 0
chr1 20591 C T 1 1 0
chr1 39182 G A 0 0 1
chr2 58201 T C 1 0 1
where, for a haploid organism, the genotypes represent a single allele per sample.
A few things I would change before you run this on real data
The script above is a good general pipeline, but I would not blindly treat QUAL >= 30, DP >= 10, GQ >= 20, and 10% missingness as universally correct.
For your actual fungal dataset, the next step should be to determine:
Genome size
Reference assembly quality
Illumina read length
Number of isolates
Expected sequencing depth
Whether the isolates are strictly haploid vs. potentially heterokaryotic
Whether the reference contains repeats/low-complexity regions
Whether you want SNPs for GWAS, population genetics, phylogeny, or identifying mutations
Those choices substantially affect the final filtering criteria.