Mercurial > repos > tyty > structurefold

comparison Iterative_mapping/iterative_map.py @ 61:563b7b27cb52 draft

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