This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Is vg giraffe appropriate for RADseq sequencing reads?

Hi vg team,

I am working with a plant pangenome generated from multiple genome assemblies. For our whole-genome resequencing data, we have been using vg giraffe successfully. The workflow is roughly:

vg giraffe \
  -Z graph.gbz \
  -x graph.xg \
  -m graph.shortread.withzip.min \
  -z graph.shortread.zipcodes \
  -d graph.dist \
  -f sample.R1.fq.gz \
  -f sample.R2.fq.gz \
  -t 8 \
  -N sample \
  -R sample \
  -o BAM \
  --ref-name Reference_name > sample.unsorted.bam

Then we sort the BAM, add read groups, mark duplicates, and run GATK HaplotypeCaller using the same linear reference FASTA corresponding to the reference path:

gatk HaplotypeCaller \
  -I sample.dedup.bam \
  -O sample.g.vcf \
  -R Reference_name.fasta \
  -ERC BP_RESOLUTION \
  --minimum-mapping-quality 30

This works for WGS data. However, we also have RADseq data from the same project. The RADseq data volume is much smaller, but when we tried to run it through the same vg giraffe workflow, the jobs ran for many days and appeared to be much slower than expected.

My questions are:

Is vg giraffe generally appropriate for RADseq or other reduced-representation sequencing data? Are there known issues when using giraffe with sparse, restriction-site-associated reads, compared with ordinary WGS short reads? If the downstream goal is SNP calling with GATK on the same reference coordinate system, would it be more appropriate to align the RADseq reads directly to the same linear reference FASTA using BWA-MEM, and then run the same GATK pipeline? If we use BWA-MEM for the RADseq samples and vg giraffe for the WGS samples, but both are ultimately aligned/projected to the same reference path and called with the same Reference_name.fasta, is that a reasonable strategy for producing VCFs on a consistent coordinate system? Are there recommended vg giraffe settings, alternative vg tools or strategies that would make RADseq alignment feasible, or is graph-based alignment usually not worth the cost for this data type?

Conceptually, I understand that for WGS the graph may improve read placement and reduce reference bias before projecting the alignments back to a linear reference path. But for RADseq, because the data only cover a small fraction of the genome, I am wondering whether the graph provides enough benefit to justify the computational cost.

Any suggestions or best practices would be appreciated.

Thanks!

vg radseq
• 65 views • link
updated by GenoMax GenoMax • written by RianeZ • RianeZ

1 answer

The thing I'd flag hardest is the mixed-aligner plan. Coordinates will match, so that part is fine, but reference bias won't. The whole reason you're running giraffe on the WGS is that the graph reduces reference bias, and BWA-MEM doesn't. So your RADseq samples would systematically under-call variants in regions that diverge from the reference path while your WGS samples don't. If those two sets ever get compared to each other -- population structure, allele frequencies, anything -- that difference is now confounded with library type. In a plant pangenome, where divergence from any single reference is the entire motivation, that isn't a small effect.

So if the two datasets get analysed together, use one aligner for both: either accept the reference bias and put everything through BWA, or pay the compute and put everything through giraffe. If they answer separate questions and never get pooled, mixing is fine.

On the speed, though -- RADseq running slower than WGS on far less data is itself odd, and I'd diagnose that rather than route around it. Couple of things worth checking. Are adapters and barcodes fully trimmed? Untrimmed reads soft-clip badly and can inflate alignment cost a lot. And RADseq reads stack up on the same restriction sites, so if some of those sites happen to sit in repetitive or bubble-dense parts of the graph you're paying that expensive placement over and over rather than once. Subsampling to ~100k reads and checking whether runtime scales linearly will tell you fairly quickly whether it's a uniform per-read cost or a handful of pathological reads stalling the run.

Given RADseq covers a small fraction of the genome, total work should be small. Days of runtime points at something specific being wrong rather than graph alignment being unsuited to the data type.

• 0 views • link
written by SeqBench • SeqBench

Log in to answer this question.

More posts like this