Interpreting Results¶

The pipeline's product is a pair of tab-separated allele frequency tables. This page describes their format precisely, works through an example, and covers the mistakes that are easy to make when reading them.

The frequency tables¶

Two files are written to Output/Frequencies/, named after vcf.fileName:

File	Contents
`<name>_snp_freq.tsv`	Every surviving SNP site
`<name>_indel_freq.tsv`	Every surviving insertion or deletion site
`<name>_snp_depth.tsv`	The read counts those SNP frequencies were computed from
`<name>_indel_depth.tsv`	The same, for indels

With the default vcf.fileName = 'Test' those are Test_snp_freq.tsv, Test_indel_freq.tsv and their two _depth.tsv counterparts.

Columns¶

CHROM   POS   REF   ALLELE   TOTAL_AD   <sample 1>   <sample 2>   …

Column	Contents
`CHROM`	Reference sequence name
`POS`	1-based position
`REF`	The reference allele for the site, repeated on every row of that site
`ALLELE`	The allele this row reports on
`TOTAL_AD`	Frequency of `ALLELE` across all samples combined
sample columns	Frequency of `ALLELE` in that sample

Sample columns appear in metadata.csv row order — see Row order decides column order.

One row per allele¶

This is the part that surprises people. A site does not occupy one row; it occupies one row per allele, including the reference allele. A biallelic SNP produces two rows, a triallelic site three.

REF is constant within a site and ALLELE varies. The reference row is the one where the two are equal.

The header says `TOTAL_AD`, the column holds a frequency¶

The fifth column is derived from INFO/AD, the cohort-wide allelic depth, and the conversion divides it by the row total exactly as it does for the sample columns. The header label is carried over from the intermediate depth file and was not renamed. Read it as overall frequency, not as a depth.

Worked example¶

Take a triallelic site with two samples:

CHROM  POS   REF  ALT    INFO/AD        sample1 AD    sample2 AD
chr1   1000  A    G,T    800,150,50     400,100,0     400,50,50

Cohort total is 800 + 150 + 50 = 1000. Sample 1 totals 500, sample 2 totals 500. The table contains:

CHROM  POS   REF  ALLELE  TOTAL_AD  sample1  sample2
chr1   1000  A    A       0.8       0.8      0.8
chr1   1000  A    G       0.15      0.2      0.1
chr1   1000  A    T       0.05      0        0.1

Every column within one site sums to 1. A zero means the allele was not observed in that sample, not that the site was missing there.

The depth tables¶

<name>_snp_depth.tsv and <name>_indel_depth.tsv are the input that conversion was applied to: the same sites, holding read counts instead of frequencies. They are what the worked example above calls the AD columns, exactly as they came out of the VCF.

They are not the same shape as the frequency tables, and reading them as if they were is the mistake to avoid. A depth table has one row per site, and each cell holds a comma-separated list of counts in REF-then-ALT order. The frequency table expands that into one row per allele. The worked example's site is a single row here:

CHROM  POS   REF  ALT  TOTAL_AD     sample1    sample2
chr1   1000  A    G,T  800,150,50   400,100,0  400,50,50

Two consequences. A depth table has fewer rows than its frequency table — one per site rather than one per allele. And its fourth column is headed ALT, listing only the alternates, where the frequency table's is ALLELE and names one allele per row including REF.

Use these when a frequency alone is not enough: a frequency of 0.5 from 400 reads and one from 2 reads are the same number and very different evidence. They are also what the analysis layer reads to compute each site's effective sample size.

Reading the tables correctly¶

REF is the major allele, not the reference genome's base. Step 7 re-encodes each site so the most-read allele becomes REF (why). This makes rows comparable across samples and runs, but it means REF will often disagree with the FASTA you supplied. If you need the assembly's base, take it from the assembly.

"Most-read" is decided across the whole cohort. A sample where the cohort-minor allele dominates locally will show a REF frequency below 0.5. That is a real signal, not a defect.

A missing row means the variant did not survive the chain, not that it is absent. Sites are removed at nine points between the FASTQ and the table. Before concluding a variant is absent from your population, check it against The Filter Chain — in particular the cross-sample requirement at stage 6, which discards alleles seen in too few pools regardless of how frequent they are in those pools.

Frequencies are read proportions, not estimates with error bars. The table reports the fraction of reads carrying each allele. Sampling error in that fraction depends on depth at the position and on the pool size, and the pipeline does not propagate it. If your analysis needs uncertainty, compute it from the depth — which is why the VCFs retain AD and DP.

There is no genotype to read. FORMAT/GT is ./. throughout, deliberately. Genotype-based tools pointed at these VCFs will find nothing.

The VCF files¶

Output/VCF/ holds the call sets. After a complete run:

File	What it is
`<name>.vcf`	Raw output of step 6 — every called site, no step 7 filtering, original reference encoding
`<name>_annotated.vcf`	SnpEff annotation of the raw call set, if `annotate = true`

Those two are the whole of it. Every VCF step 7 produces — _sort, _sort_fp, _sort_fp_dq, and the split _snp and _indel files — is consumed by the next process and deleted, so none of them survives a finished run. See Steps delete their own inputs.

The filtered call set is therefore not kept as a VCF. What it becomes is the frequency tables. If you need the filtered sites in VCF form, the raw call set plus the positions in the tables is what you have to work from.

Annotation is applied to the unfiltered call set

Step 8 runs on the output of step 6, in parallel with the frequency branch rather than after it. So <name>_annotated.vcf contains sites that the step 7 filters removed, and its allele encoding is the original reference-based one, not the major-allele normalized one. To attach annotations to your frequency tables, join on CHROM/POS and expect unmatched rows on the annotation side.

The reports¶

Output/Reports/ collects everything generated along the way.

Path	Produced by	Useful for
`Alignment/<sample>_alignment_report.txt`	`bamtools stats`	Mapping rate, duplicate rate, paired-end statistics
`Coverage/<sample>_coverage_report.txt`	`samtools coverage`	Per-contig depth and breadth
`Depth/<sample>_depth_histogram.tsv`	`samtools stats`	The sample's whole depth distribution, one line per depth
`Depth/<sample>_depth_report.txt`	PoolSeqFlow	The ceiling put on that sample, and why
`Fastqc/<sample>/`	FastQC	Raw, trimmed and clipped read quality
`Trimming/<sample>/`	Trim Galore	How much was removed, and which adapter was detected
`snpeff_summary.html`	SnpEff	Variant effect summary, if annotation ran
`snpeff_summary.genes.txt`	SnpEff	The same counts per gene and transcript, tab separated
`PoolSeqFlow_pipeline_report.html`	Nextflow	Per-task resource usage
`PoolSeqFlow_pipeline_timeline.html`	Nextflow	Where wall-clock time went
`PoolSeqFlow_pipeline_trace.txt`	Nextflow	Machine-readable task trace
`PoolSeqFlow_pipeline_dag.html`	Nextflow	Workflow graph

run_parameters.txt is not in here — it sits one level up, at the root of Output/, beside the results rather than among the reports about them. It is a readable record of the analysis parameters these outputs were built from.

The three worth reading on every run are the coverage report, the depth report and run_parameters.txt. The depth report is the one that is easy to skip and should not be: the capped BAM it describes is transient, so this file is the only record of which reads reached bcftools for that sample.

If you ran a table of several runs, the four Nextflow reports describe the whole invocation rather than any one run, and are written once under Output/All_Runs/Reports/.

A quick sanity pass¶

After a run finishes, four checks catch most problems:

Sample columns. Does the table have the number of columns you expect, in the order you laid out in metadata.csv? A count lower than expected means rows were merged by a shared RG_Sample (why).
What was capped. grep -H 'ceiling applied' Output/Reports/Depth/*_depth_report.txt. A ceiling far below the sample's typical depth, or a sample left uncapped when you expected otherwise, are both worth a look at its histogram before you trust the frequencies.
Row counts. Compare the site count in <name>.vcf with the distinct positions that reached the tables — tail -n +2 <name>_snp_freq.tsv | cut -f1,2 | sort -u | wc -l, and the same for the indel table. A very large drop points at the cross-sample filter; check sampleThreshold against your sample count in the table above.
Column sums. Frequencies within a site should sum to 1 in every column.