Samtools pileup extract read length for position

Question

[rah]

I am working with bam files and vcf files, where i wish to extract the read length for those reads which overlaps with the positions on the vcf file. To do this i have simply created a txt file with the positions and then i have tried with samtools.

I have tried with samtools mpileup, but i can only get it to give me the position in the read using --output-BP. I can also get it to output additional fields or tags, such as FLAG,NM,SQ

samtools mpileup file.bam -a -l pos.txt --output-BP --output-extra FLAG,NM,SQ

And if i look into https://samtools.github.io/hts-specs/SAMtags.pdf to see other tags which i might find useful, such as the cigar (CIGAR tag) or actual alignment (OA tag) which i then can just print the length of it fails.

So my question is: 1) Can anyone suggest a way to extract the read length for reads in a bam file matching either the positions in a .txt, .bed file or .vcf file?

Answer 1

[Pierre Lindenbaum]

samtools view -M -L <(bcftools query -f '%CHROM\t%POS0\t%END'  in.vcf) in.bam | awk '{print length($10);}'

[swbarnes2]

Another way; use bedtools to get the intersect of your bam with your vcf, then use the awk above to get the read lengths.