def go(args):
bed = read_bed_file(args.bedfile)
infile = pysam.AlignmentFile(args.alignment, "rb")
for s in infile:
#print s.get_aligned_pairs()
#print ">%s\n%s" % (s.query_name, s.query_alignment_sequence)
p1 = find_primer(bed, s.reference_start, '+')
p2 = find_primer(bed, s.reference_end, '-')
primer_start = p1[2]['start']
# start is the 5'
primer_end = p2[2]['start']
query_align_start = find_query_pos(s, primer_start)
query_align_end = find_query_pos(s, primer_end)
print >> sys.stderr, "%s\t%s\t%s\t%s" % (primer_start, primer_end, primer_end - primer_start, s.query_length)
startpos = max(0, query_align_start - 40)
endpos = min(query_align_end+40, s.query_length)
print ">%s\n%s" % (s.query_name, s.query_sequence[startpos:endpos])
#query_align_end + 30])
评论列表
文章目录