Hello, I'm new to variant detection pipeline and came across a problem. I have amplicon sequencing data from control and trasngenic sample. The gene we are targeting is a multi-copy one and exact number of alleles possible are not know, but Sanger seqeuncing showed it should be between 6-9 copies. We have done Oxford Nanopore Sequencing, and found few of the variants matched whatever was known from Sanger. The reference panel was made with variants detected from both Sanger and ONT. We now have total 9 reference panel to begin with, Now, for my control and transgenic sample I'm using similar pipeline of QC, trimming, alignment witht the 9 reference panel, index, sort, marking duplicates, again sort. After that I ran freebayes for variant calling, indexing the vcf file and finally generated the consensus. However, control sample didnt have any indels as was expected.
For trasngenic samples, running freebayes with default didn't give any meaningful result, Then I ran with "--pooled-continuous", "--min-alternate-count 1", "--haplotype-length 0", "--report-monomorphic" flag. And finally the consesus generated using bcftools consensus with the "--haplotype A" flag. When aligned with the reference panel (same specific haplotype) and the control (same specific haplotype) showed +A insertion around the cut site. But we are expecting a +C insertion. My logic was, since the edits were heterozygous, allele frequency may be quite low to be detected, so dropping the parameters to the lowest value could be helpful. But why the expected indels cannot be detected? From my pipeline do you see any major methodological flaws? What
Two things in that pipeline are probably discarding your reads before freebayes ever sees them.
MarkDuplicates on amplicon data flags nearly everything, since every read starts and ends at the same place by design, and freebayes skips duplicate-flagged reads unless you pass --use-duplicate-reads. Then the 9-copy panel: reads that map equally well to several near-identical haplotypes come out at MAPQ 0, and freebayes drops those too by default. I'd check the depth freebayes actually used at the cut site rather than the depth sitting in the BAM.
Separate problem - freebayes isn't built for nanopore error profiles, and +A versus +C at a cut site is exactly where that bites, because homopolymer indels are ONT's signature error. For CRISPR amplicons I'd use CRISPResso2 instead. It quantifies indels around the cut site directly and skips the consensus step entirely.
Also, bcftools consensus -H A doesn't pick a haplotype from your panel, it just takes the first ALT at each site.
Thanks for yoru reply. This is actually 150bp illumina paired-end short-read data, not nanopore. Nanopore reads are used for building the reference panel. I think Crispresso can handle data where the reads span whole of the amplicon length, not in my case where the reads barely span.
~99% of the reads spanning the putative cut site have MAPQ=0, further supporting the possibility of extensive cross-mapping between highly similar haplotypes which you mentioned. I tried aligning reads independently to each haplotype and extracting uniquely mapped reads, but neither approach substantially improved mapping specificity.
the cut site also shows approximately half the overall BAM depth and corresponds to the insertion region present in at least 4 of the 9 haplotypes.
Would you suggest a better strategy to distinguish true haplotype-specific/cut-site reads from cross-mapped reads?
The broader aim of this experiement was to assess whether Illumina short-read sequencing can provide sufficient resolution to partially substitute for ONT, particularly as a cost-effective alternative in our case, since we already found this working for a comparatively simpler gene with known number of haplotypes.
Log in to answer this question.