Frage zum Python Skript mit PySAM
Verfasst: Sonntag 21. August 2016, 19:12
Hallo,
ich soll innerhalb eines Studienprojektes den Workflow eines anderen Projektes nachvollziehen. In einem Schritt kommt ein Python Skript zum Einsatz. Darin wird der Input und Output so angegeben:
Ich hatte herausgefunden, dass der "-" für den Standard-Input bzw. Output steht. Bei den Optionen steht r für read, b für BAM Format, w für write und h für das TAM Format. Soweit so gut. Nun hatte ich das mit dem Befehl in der Form "python meinskript meininput.bam meinoutput.tam" gestartet. Nach dem Abschicken des Befehls sieht es so aus als ob er den Befehl ausführen würde, aber er kommt scheinbar nie zum Ende (auch nach mehreren Stunden nicht). Man kann das Ganze auch nicht mit Strg+C beenden. Nun zu meiner Frage: Muss der Befehl zum ausführen des Skriptes in der Form angegen werden wie ich es getan habe oder hab ich da etwas fehlinterpretiert? Ich hatte nämlich noch nie zuvor mit Python zu tun gehabt, von daher bin ich mir da nicht sicher. Ansonsten muss die Fehlerquelle wo anders liegen und ich muss woanders suchen. Vielen Dank schon mal. Hier ist nochmal der ganze Quellcode.
ich soll innerhalb eines Studienprojektes den Workflow eines anderen Projektes nachvollziehen. In einem Schritt kommt ein Python Skript zum Einsatz. Darin wird der Input und Output so angegeben:
Code: Alles auswählen
infile = pysam.AlignmentFile("-", "rb")
outfile = pysam.AlignmentFile("-", "wh", template=infile)Code: Alles auswählen
#!/usr/bin/env python
import pysam
import sys
from copy import copy
from vcftagprimersites import read_bed_file
def trim(s, start_pos, end):
if not end:
pos = s.pos
else:
pos = s.reference_end
eaten = 0
while 1:
## chomp stuff off until we reach pos
if end:
flag, length = cigar.pop()
else:
flag, length = cigar.pop(0)
#print >>sys.stderr, "Chomped a %s, %s" % (flag, length)
if flag == 0:
## match
#to_trim -= length
eaten += length
if not end:
pos += length
else:
pos -= length
if flag == 1:
## insertion to the ref
#to_trim -= length
eaten += length
if flag == 2:
## deletion to the ref
#eaten += length
if not end:
pos += length
else:
pos -= length
pass
if flag == 4:
eaten += length
if not end and pos >= start_pos and flag == 0:
break
if end and pos <= start_pos and flag == 0:
break
#print >>sys.stderr, "pos:%s %s" % (pos, start_pos)
extra = abs(pos - start_pos)
#print >> sys.stderr, "extra %s" % (extra)
if extra:
if flag == 0:
#print >>sys.stderr, "Inserted a %s, %s" % (0, extra)
if end:
cigar.append((0, extra))
else:
cigar.insert(0, (0, extra))
eaten -= extra
if not end:
s.pos = pos - extra
#print >>sys.stderr, "New pos: %s" % (s.pos)
if end:
cigar.append((4, eaten))
else:
cigar.insert(0, (4, eaten))
oldcigarstring = s.cigarstring
s.cigartuples = cigar
#print >>sys.stderr, s.query_name, oldcigarstring[0:50], s.cigarstring[0:50]
bed = read_bed_file('all')
def find_primer(pos, direction):
# {'Amplicon_size': '1874', 'end': 7651, '#Region': 'region_4', 'start': 7633, 'Coords': '7633', "Sequence_(5-3')": 'GCTGGCCCGAAATATGGT', 'Primer_ID': '16_R'}
from operator import itemgetter
closest = min([(abs(p['start'] - pos), p['start'] - pos, p) for p in bed if p['direction'] == direction], key=itemgetter(0))
return closest
infile = pysam.AlignmentFile("-", "rb")
outfile = pysam.AlignmentFile("-", "wh", template=infile)
for s in infile:
cigar = copy(s.cigartuples)
if len(sys.argv) > 1:
if not s.query_name.startswith(sys.argv[1]):
continue
## logic - if alignment start site is _before_ but within X bases of
## a primer site, trim it off
if s.is_unmapped:
continue
p1 = find_primer(s.reference_start, 'F')
p2 = find_primer(s.reference_end, 'R')
print >>sys.stderr, "%s\t%s\t%s_%s\t%s\t%s\t%s\t%s" % (s.reference_start, s.reference_end, p1[2]['Primer_ID'], p2[2]['Primer_ID'], p1[2]['Primer_ID'], abs(p1[1]), p2[2]['Primer_ID'], abs(p2[1]))
#amp_len = p1['end'] - p2['end']
## if the alignment starts before the end of the primer, trim to that position
#print >>sys.stderr, s.reference_start, p1[2]['end']
#print >>sys.stderr, s.reference_end, p2[2]['start']
if p1[0] < 50:
if s.reference_start < p1[2]['start'] - 1:
# trim(s, p1[2]['end']-1, 0)
trim(s, p1[2]['start']-1, 0)
else:
trim(s, s.reference_start + 20, 0)
if p2[0] < 50:
if s.reference_end > p2[2]['end'] - 1:
# trim(s, p2[2]['start']-1, 1)
trim(s, p2[2]['end']-1, 1)
else:
trim(s, s.reference_end - 20, 1)
# trim(s, 50, 0)
# trim(s, 50, 1)
outfile.write(s)