Get length of the alignment from bam file

Question

[David]

Hi, I would like to know if it is it possible to get the read length of the higher alignment in a bam file having more than 85% identity.

For example an output with this fields:

• Read_length
• Alignment_length
• Percent_identity

thanks

[WouterDeCoster]

You can do this with pysam. Do you have any (Python) coding experience?

[David]

I´m more of a perl guy but can give it a try

Answer 1

[Pierre Lindenbaum]

using bioalcidaejdk : http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

• stream each record
• discard unmapped record or record without cigar

• print : the read length using https://samtools.github.io/htsjdk/javadoc/htsjdk/htsjdk/samtools/Cigar.html#getReadLength-- , the reference length using https://samtools.github.io/htsjdk/javadoc/htsjdk/htsjdk/samtools/Cigar.html#getPaddedReferenceLength-- and the % of identify using the NM Attribute (must be present) https://samtools.github.io/hts-specs/SAMtags.pdf

  java -jar dist/bioalcidaejdk.jar -e 'stream().filter(R->!R.getReadUnmappedFlag() && R.getCigar()!=null).forEach(R->println(R.getReadLength()+" "+R.getCigar().getPaddedReferenceLength()+" "+((double)(R.getReadLength()- R.getIntegerAttribute("NM")))/R.getReadLength()));' in.bam
  

[vinaykusuma]

Hi, Does bioaclidaejdk get the output of those three column in same order as present in sam file. So can i apped read id information row by row to the output.?