view Iterative_mapping/iterative_map.py @ 44:5d7cdbf181a3 draft

Uploaded

#!/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*")