python/salus/alignbam.py

   1 #!/usr/bin/env python3
   2
   3 import sys
   4 import pysam
   5 from skbio import DNA
   6 from skbio.alignment import global_pairwise_align_nucleotide
   7 import alignfq
   8
   9 def main():
  10     if len(sys.argv) < 2 :
  11         print('Usage:',sys.argv[0],'<bam file> >alignments.out',file=sys.stderr,flush=True);
  12         exit(0);
  13     bamFile = sys.argv[1]
  14     Adapters = alignfq.getAdapters()
  15
  16     with pysam.AlignmentFile(bamFile, "rb") as inbam:
  17         for aln in inbam:
  18             # aln.query_name
  19             fqSeq = aln.get_forward_sequence()
  20             refNfo = ''.join((
  21                 str(aln.flag),'|',
  22                 aln.reference_name,':',
  23                 str(aln.reference_start),'+',
  24                 str(aln.reference_length),
  25                 ':',aln.cigarstring
  26             ))
  27             faID = ' '.join((
  28                 aln.query_name,
  29                 refNfo
  30             ))
  31             """
  32             alnSeq = aln.seq
  33             if aln.is_reverse:
  34                 alnSeq = str(DNA(alnSeq).reverse_complement())
  35             if fqSeq != alnSeq:
  36                 eprint('-'.join('[x]',fqSeq,alnSeq))
  37             """
  38             maxAlignments = alignfq.aliOneSeq(Adapters,fqSeq)
  39             alignfq.printAli(maxAlignments,fqSeq,faID)
  40
  41 if __name__ == "__main__":
  42     main()  # time ./alignbam.py t.multi.bam 2>/dev/null >t.multi.bam.adapterali