Mercurial > repos > tyty > structurefold

comparison Iterative_mapping/iterative_map.py @ 0:d56631911cc1 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author tyty
date Mon, 15 Sep 2014 14:41:13 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:d56631911cc1
1 #!/usr/bin/env python
2 # -*- coding: utf-8 -*-
3
4 import sys
5 import os
6 from read_file import *
7 from read_s_file import *
8
9 type_input = sys.argv[1]
10 seq_file = sys.argv[2]
11 ref_file = sys.argv[3]
12 shift = sys.argv[4]
13 length = sys.argv[5]
14 map_type = sys.argv[6]
15 output_file = sys.argv[7]
16
17
18 if map_type!="default":
19 s = ""
20 s = s+"-v "+sys.argv[8]
21 s = s+" -5 "+sys.argv[9]
22 s = s+" -3 "+sys.argv[10]
23 s = s+" -k "+sys.argv[11]
24 if sys.argv[12]:
25 s = s+" -a"
26 if int(sys.argv[13])>=1:
27 s = s+" -m "+sys.argv[13]
28 if sys.argv[14]:
29 s = s+" --best --strata "
30
31 else:
32 s = "-v 3 -a --best --strata "
33
34 ospath = os.path.realpath(sys.argv[0])
35 ost = ospath.split('/')
36 syspath = ""
37 for i in range(len(ost)-1):
38 syspath = syspath+ost[i].strip()
39 syspath = syspath+'/'
40
41 os.system("bowtie-build -f "+ref_file+" "+syspath+"ref > "+syspath+"log.txt")
42
43 os.system("cp "+seq_file+" "+syspath+"seq0.fa")
44
45 if type_input == "fasta":
46 tp = 'fasta'
47 if type_input == "fastq":
48 tp = 'fastq'
49
50 k = 0
51 print(type_input)
52 while(True):
53 if type_input == "fasta":
54 os.system("bowtie "+s+"-f "+syspath+"ref"+" "+syspath+"seq"+str(k)+".fa --quiet -S > "+syspath+"map"+str(k)+".sam")
55 if type_input == "fastq":
56 os.system("bowtie "+s+"-q "+syspath+"ref"+" "+syspath+"seq"+str(k)+".fa --quiet -S > "+syspath+"map"+str(k)+".sam")
57 os.system("samtools view -Sb -F 0xfff "+syspath+"map"+str(k)+".sam > "+syspath+"mapped"+str(k)+".bam 2>"+syspath+"log.txt") #get mapped reads
58 os.system("samtools view -Sb -f 0x4 "+syspath+"map"+str(k)+".sam > "+syspath+"umapped"+str(k)+".bam 2>"+syspath+"log.txt") #get unmapped reads
59 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
60 os.system("samtools merge -f "+syspath+"unmapped"+str(k)+".bam "+syspath+"umapped"+str(k)+".bam "+syspath+"rmapped"+str(k)+".bam") #get reversed mapped reads
61 os.system("samtools view -h -o "+syspath+"unmapped"+str(k)+".sam "+syspath+"unmapped"+str(k)+".bam") #get reversed mapped reads
62 if k>0:
63 os.system("samtools view -h -o "+syspath+"mapped"+str(k)+".sam "+syspath+"mapped"+str(k)+".bam") #get reversed mapped reads
64 os.system("cut -f 1 "+syspath+"unmapped"+str(k)+".sam > "+syspath+"unmapped"+str(k)+".txt")
65 os.system("cut -f 1 "+syspath+"mapped"+str(k)+".sam > "+syspath+"mapped"+str(k)+".txt")
66 os.system("python "+syspath+"remove_map.py "+syspath+"unmapped"+str(k)+".txt "+syspath+"mapped"+str(k)+".txt "+syspath+"runmapped"+str(k)+".txt")
67 os.system("rm "+syspath+"mapped"+str(k)+".sam")
68 os.system("rm "+syspath+"mapped"+str(k)+".txt")
69 os.system("rm "+syspath+"unmapped"+str(k)+".txt")
70 else:
71 os.system("cut -f 1 "+syspath+"unmapped"+str(k)+".sam > "+syspath+"runmapped"+str(k)+".txt")
72
73 os.system("rm "+syspath+"unmapped"+str(k)+".bam")
74 os.system("rm "+syspath+"umapped"+str(k)+".bam")
75 os.system("rm "+syspath+"rmapped"+str(k)+".bam")
76 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
77 os.system("python "+syspath+"truncate.py "+syspath+"unmap_seq"+str(k)+".fa "+shift+" "+syspath+"seq"+str(k+1)+".fa "+length) #truncate unmapped sequence
78 os.system("rm "+syspath+"seq"+str(k)+".fa") #Remove sequences being mapped
79 os.system("rm "+syspath+"map"+str(k)+".sam") #Remove mapping file
80 os.system("rm "+syspath+"unmap_seq"+str(k)+".fa") #Remove unmapped sequnce
81 os.system("rm "+syspath+"runmapped"+str(k)+".txt")
82 os.system("rm "+syspath+"unmapped"+str(k)+".sam")
83
84 os.system("wc -l "+syspath+"seq"+str(k+1)+".fa > "+syspath+"count"+str(k+1)+".txt")
85 c = read_sp_file(syspath+"count"+str(k+1)+".txt")
86 if c[0][0] == '0': #If no reads is in the sequence file, stop
87 os.system("rm "+syspath+"count"+str(k+1)+".txt")
88 os.system("rm "+syspath+"seq"+str(k+1)+".fa")
89 break
90 os.system("rm "+syspath+"count"+str(k+1)+".txt")
91 k = k+1
92
93 ss = ""
94 for i in range(0,k+1):
95 ss = ss+" "+syspath+"mapped"+str(i)+".bam"
96
97
98 os.system("samtools merge -f "+output_file+" "+ss)
99 #print("samtools merge mapped_all.bam"+ss)
100 os.system("rm "+syspath+"mapped*.bam")
101 os.system("rm "+syspath+"ref*")
102
103