@@ -0,0 +1,1575 @@
+%0A# given input from compare_bams, generate a sam file with matching mismatches based on specified criteria%0A# e.g. python extract_mismatched_reads.py --min_distance 50 --mappable True ../../data/e-coli-mg1655_snp_ins_del_50_x1_bwa.bam %3C ../../data/compare_bams3.log %3E mappable_test.sam%0A%0Aimport argparse%0Aimport sys%0A%0Aimport bio%0Aimport config%0A%0Aparser = argparse.ArgumentParser(description='Find mismatched reads matching criteria')%0Aparser.add_argument('bam', help='bam files to extract reads')%0Aparser.add_argument('--min_distance', dest='min_distance', type=int, default=0, help='only consider alignments with at least this distance from the true location')%0Aparser.add_argument('--max_distance', dest='max_distance', type=int, default=1e12, help='only consider alignments with at most this distance from the true location')%0Aparser.add_argument('--mappable', dest='mappable', type=bool, default=False, help='only consider alignments that have a mappable location')%0A%0Aargs = parser.parse_args()%0A%0Aallowed_tags = set()%0Afor line in sys.stdin:%0A  fields = line.strip().split(',')%0A  if len(fields) == 5 and fields%5B0%5D.isdigit():%0A    distance = int(fields%5B3%5D)%0A    mappable = int(fields%5B4%5D) == 0%0A    if distance %3E= args.min_distance and distance %3C= args.max_distance:%0A      if mappable and args.mappable or not mappable and not args.mappable:%0A        allowed_tags.add( int( fields%5B2%5D ) )%0Abio.log_stderr( '%25i allowed tags' %25 len(allowed_tags) )%0A%0Abio.SamFilter( sam_fh=bio.BamReaderExternal( config.BAM_TO_SAM, args.bam ), target_fh=sys.stdout, allowed_tags=allowed_tags, log=bio.log_stderr )%0A

build filtered sam file from compare bam output