This is a test version of Biostars. For the public version, visit https://www.biostars.org.
call variants for every read

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).

alignment snp next-gen
• 2,042 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by mimakaev • mimakaev

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

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

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.

• 0 views • link
written by mimakaev • mimakaev

1 answer

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"
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this