Plotting contig length after a MIRA run / iteration
0
0
Entering edit mode
7.0 years ago
Burgenix • 0

Hi,

Running MIRA for a contig threshold here. I am trying to plot the contig length and the number of iterations, but somehow, the plotter stops the program after a iteration. Any idea why?

import os
import sys
import shutil
import subprocess
import matplotlib.pyplot as plt
import numpy as np

name = sys.argv[1] + '-0'
ref_loc = sys.argv[2]
read_loc = sys.argv[3]

manifest = open('manifest-0.conf', 'w')
manifest.write('project = ' + name + '\njob = genome,mapping,accurate\nparameters = -GE:not=6 -NW:mrnl=0 -AS:nop=1 SOLEXA_SETTINGS -CO:msr=no COMMON_SETTINGS -SB:tor=no\nreadgroup\nis_reference\ndata = ' + ref_loc + '\nstrain = ' + name + '\nreadgroup = reads\ndata = ' + read_loc + '\ntechnology=solexa')
manifest.close()

avg_cov = 0
count = 0
new_avg_cov = 0
contig_length = 0

while new_avg_cov <= avg_cov * 2 or avg_cov == 0: #set threshold for average coverage points

    command1 = 'mira manifest-' + str(count) + '.conf'
    print(command1)

    avg_cov = new_avg_cov

    subprocess.call(command1, shell=True)

    print('MIRA is really done!')

    with open(name + '_assembly/' + name +'_d_info/'+ name +'_info_contigstats.txt') as f:
        print('Fishing for coverage:')
        for line in f:
            if not line.startswith('#'):
                new_avg_cov = float(line.split('\t')[5])


    with open(name + '_assembly/' + name +'_d_info/'+ name +'_info_contigstats.txt') as f:
        if not line.startswith('#'):
            contig_length = int(line.split('\t')[1])

    plt.scatter(count, contig_length)
    plt.xlabel("Number of iterations", fontsize=14)
    plt.ylabel("Length of assembled contig", fontsize=14)
    plt.axis([0,100, 0,200000])
    plt.pause(0.001)
    plt.show()

    new_bait = name + '_assembly/' + name +'_d_results/'+ name +'_out_AllStrains.unpadded.fasta'
    print('new bait is ' + new_bait)
    count += 1

    name = name.split('-')[0] + '-' + str(count)
    manifest = open('manifest-' + str(count) +'.conf', 'w')
    print('new manifest is ' + 'manifest-' + str(count) +'.conf')
    manifest.write('project = ' + name +'\njob = genome,mapping,accurate\nparameters = -GE:not=6 -NW:mrnl=0 -AS:nop=1 SOLEXA_SETTINGS -CO:msr=no COMMON_SETTINGS -SB:tor=no\nreadgroup\nis_reference\ndata = /home/radu/Desktop/new_test' + new_bait + '\nstrain = ' + name + '\nreadgroup = reads\ndata = ' + read_loc + '\ntechnology=solexa')
    manifest.close()

    if os.path.exists('./manifest-' + str(count-2) + '.conf'):
        os.remove('./manifest-' + str(count-2) + '.conf')
        shutil.rmtree('./' + name.split('-')[0] + '-' + str(count-2) + '_assembly/' )

Tried to use a line instead of a scatterplot, same result. I think It might be a problem with the manifest file, but i cant figure it out.

python matplotlib • 1.3k views
ADD COMMENT
0
Entering edit mode

Have you tried running the program without plt.show or plt.pause? Both of these functions may cause the program to block until some GUI action is taken. Similarly, have you checked that the conditions in the while loop are not satisfied?

ADD REPLY

Login before adding your answer.

Traffic: 1571 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6