call variants for every read

Question

[mimakaev]

My goal is to align Illumina reads to the reference genome and look at variants within each read. The eventual goal is to look at co-occurrence of variants on the two side of the paired end read, in a metagenomic sample.

I found many tools that would call variants from the whole sample: iVar, bcftools+samtools, etc. However, I need variants called for every read

I also found that minimap2 can output the "--cs" tag in bam files. But that seems to be far from the proper variant calls (A25101G would be what I need).

[Pierre Lindenbaum]

so you want to look at each paired read and detect each SNV occuring in the overlapping interval ?

[mimakaev]

I want to look at each paired read and detect each SNV occurring in the read. Paired end sequencing is used just to increase a chance of a read hitting several SNPs with one read.

Answer 1

[Pierre Lindenbaum]

I wrote http://lindenb.github.io/jvarkit/Biostar489074.html it's a SNV caller which only use the overlapping part of the paired-end reads. It only report the number of time the ALT allele was found (but doesn't report the number of time the REF allele was found)

input must be sorted on read name using samtools sort -n or samtools collate

it will be much faster is the reads belong to one chromosome

Example

$ samtools view -O BAM --reference "ref.fasta" in.cram "chr22:41201525-41490147" |\
    samtools collate -O -u - |\
    java -jar dist/biostar489074.jar --reference "ref.fasta"