changeset 8:74ae6f6a3383 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Iterative_mapping/iterative_map.py	Mon Oct 20 01:33:11 2014 -0400
@@ -0,0 +1,103 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import sys
+import os
+from read_file import *
+from read_s_file import *
+
+type_input = sys.argv[1]
+seq_file = sys.argv[2]
+ref_file = sys.argv[3]
+shift = sys.argv[4]
+length = sys.argv[5]
+map_type = sys.argv[6]
+output_file = sys.argv[7]
+
+
+if map_type!="default":
+    s = ""
+    s = s+"-v "+sys.argv[8]
+    s = s+" -5 "+sys.argv[9]
+    s = s+" -3 "+sys.argv[10]
+    s = s+" -k "+sys.argv[11]
+    if sys.argv[12]:
+        s = s+" -a"
+    if int(sys.argv[13])>=1:
+        s = s+" -m "+sys.argv[13]
+    if sys.argv[14]:
+        s = s+" --best --strata "
+    
+else:
+    s = "-v 3 -a --best --strata "
+
+ospath = os.path.realpath(sys.argv[0])
+ost = ospath.split('/')
+syspath = ""
+for i in range(len(ost)-1):
+    syspath = syspath+ost[i].strip()
+    syspath = syspath+'/'
+
+os.system("bowtie-build -f "+ref_file+" "+syspath+"ref > "+syspath+"log.txt")
+
+os.system("cp "+seq_file+" "+syspath+"seq0.fa")
+
+if type_input == "fasta":
+    tp = 'fasta'
+if type_input == "fastq":
+    tp = 'fastq'
+
+k = 0
+print(type_input)
+while(True):
+    if type_input == "fasta":
+        os.system("bowtie "+s+"-f "+syspath+"ref"+" "+syspath+"seq"+str(k)+".fa --quiet -S > "+syspath+"map"+str(k)+".sam")
+    if type_input == "fastq":
+        os.system("bowtie "+s+"-q "+syspath+"ref"+" "+syspath+"seq"+str(k)+".fa --quiet -S > "+syspath+"map"+str(k)+".sam")
+    os.system("samtools view -Sb -F 0xfff "+syspath+"map"+str(k)+".sam > "+syspath+"mapped"+str(k)+".bam 2>"+syspath+"log.txt") #get mapped reads
+    os.system("samtools view -Sb -f 0x4 "+syspath+"map"+str(k)+".sam > "+syspath+"umapped"+str(k)+".bam 2>"+syspath+"log.txt") #get unmapped reads
+    os.system("samtools view -Sb -f 0x10 "+syspath+"map"+str(k)+".sam > "+syspath+"rmapped"+str(k)+".bam 2>"+syspath+"log.txt") #get reversed mapped reads
+    os.system("samtools merge -f "+syspath+"unmapped"+str(k)+".bam "+syspath+"umapped"+str(k)+".bam "+syspath+"rmapped"+str(k)+".bam") #get reversed mapped reads
+    os.system("samtools view -h -o "+syspath+"unmapped"+str(k)+".sam "+syspath+"unmapped"+str(k)+".bam") #get reversed mapped reads
+    if k>0:
+        os.system("samtools view -h -o "+syspath+"mapped"+str(k)+".sam "+syspath+"mapped"+str(k)+".bam") #get reversed mapped reads
+        os.system("cut -f 1 "+syspath+"unmapped"+str(k)+".sam > "+syspath+"unmapped"+str(k)+".txt")
+        os.system("cut -f 1 "+syspath+"mapped"+str(k)+".sam > "+syspath+"mapped"+str(k)+".txt")
+        os.system("python "+syspath+"remove_map.py "+syspath+"unmapped"+str(k)+".txt "+syspath+"mapped"+str(k)+".txt "+syspath+"runmapped"+str(k)+".txt")
+        os.system("rm "+syspath+"mapped"+str(k)+".sam")
+        os.system("rm "+syspath+"mapped"+str(k)+".txt")
+        os.system("rm "+syspath+"unmapped"+str(k)+".txt")
+    else:
+        os.system("cut -f 1 "+syspath+"unmapped"+str(k)+".sam > "+syspath+"runmapped"+str(k)+".txt")
+    
+    os.system("rm "+syspath+"unmapped"+str(k)+".bam")
+    os.system("rm "+syspath+"umapped"+str(k)+".bam")
+    os.system("rm "+syspath+"rmapped"+str(k)+".bam")
+    os.system("python "+syspath+"seq_track.py "+syspath+"runmapped"+str(k)+".txt "+syspath+"seq"+str(k)+".fa "+syspath+"unmap_seq"+str(k)+".fa "+tp) #get unmapped sequence
+    os.system("python "+syspath+"truncate.py "+syspath+"unmap_seq"+str(k)+".fa "+shift+" "+syspath+"seq"+str(k+1)+".fa "+length) #truncate unmapped sequence
+    os.system("rm "+syspath+"seq"+str(k)+".fa") #Remove sequences being mapped
+    os.system("rm "+syspath+"map"+str(k)+".sam") #Remove mapping file
+    os.system("rm "+syspath+"unmap_seq"+str(k)+".fa") #Remove unmapped sequnce
+    os.system("rm "+syspath+"runmapped"+str(k)+".txt")
+    os.system("rm "+syspath+"unmapped"+str(k)+".sam")
+    
+    os.system("wc -l "+syspath+"seq"+str(k+1)+".fa > "+syspath+"count"+str(k+1)+".txt")
+    c = read_sp_file(syspath+"count"+str(k+1)+".txt")
+    if c[0][0] == '0': #If no reads is in the sequence file, stop
+        os.system("rm "+syspath+"count"+str(k+1)+".txt")
+        os.system("rm "+syspath+"seq"+str(k+1)+".fa")
+        break
+    os.system("rm "+syspath+"count"+str(k+1)+".txt")
+    k = k+1
+
+ss = ""
+for i in range(0,k+1):
+    ss = ss+" "+syspath+"mapped"+str(i)+".bam"
+
+
+os.system("samtools merge -f "+output_file+" "+ss)
+#print("samtools merge mapped_all.bam"+ss)
+os.system("rm "+syspath+"mapped*.bam")
+os.system("rm "+syspath+"ref*")
+
+