Take advantage of newer pysam version

lib/mutfilter/realignment_filter.py

@@ -46,41 +46,6 @@ def math_log_fisher_pvalue(self,fisher_pvalue):
         return val
 
 
-    ############################################################
-    def makeTwoReference(self, chr,start,end,ref,alt, output):
-
-        hOUT = open(output, 'w')
-
-        seq = ""
-        label = ','.join([chr, str(start), str(end), ref, alt])
-        range = chr + ":" + str(int(start) - self.window + 1) +"-"+ str(int(end) + self.window)
-        for item in pysam.faidx(self.reference_genome, range):
-            # if item[0] == ">": continue
-            seq = seq + item.rstrip('\n').upper()
-        seq = seq.replace('>', '')
-        seq = seq.replace(range.upper(), '')
-
-        if re.search(r'[^ACGTUWSMKRYBDHVN]', seq) is not None:
-            print >> sys.stderr, "The return value in get_seq function includes non-nucleotide characters:"
-            print >> sys.stderr, seq
-            sys.exit(1)
-
-        print >> hOUT, '>' + label + "_ref"
-        print >> hOUT, seq
-
-        # for insertion
-        if ref == "-":   seq = seq[0:(self.window + 1)] + alt + seq[-self.window:]
-        # for deletion
-        elif alt == "-": seq = seq[0:self.window] + seq[-self.window:]
-         # for SNV
-        else:            seq = seq[0:self.window] + alt + seq[-self.window:]
-
-        print >> hOUT, '>' + label + "_alt"
-        print >> hOUT, seq
-
-        hOUT.close()
-
-
     ############################################################
     def extractRead(self, bamfile, chr,start,end, output):
 
@@ -158,10 +123,15 @@ def makeTwoReference(self, chr,start,end,ref,alt, output):
         label = ','.join([chr, str(start), str(end), ref, alt])
         range = chr + ":" + str(int(start) - self.window + 1) +"-"+ str(int(end) + self.window)
         for item in pysam.faidx(self.reference_genome, range):
-            if item[0] == ">": continue
+            # if item[0] == ">": continue
             seq = seq + item.rstrip('\n').upper()
-            seq = seq.replace('>', '')
-            seq = seq.replace(range, '')
+        seq = seq.replace('>', '')
+        seq = seq.replace(range.upper(), '')
+
+        if re.search(r'[^ACGTUWSMKRYBDHVN]', seq) is not None:
+            print >> sys.stderr, "The return value in get_seq function includes non-nucleotide characters:"
+            print >> sys.stderr, seq
+            sys.exit(1)
 
         print >> hOUT, '>' + label + "_ref"
         print >> hOUT, seq