Question: Plotting contig length after a MIRA run / iteration
0
gravatar for Burgenix
19 months ago by
Burgenix0
Burgenix0 wrote:

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.

matplotlib python • 484 views
ADD COMMENTlink written 19 months ago by Burgenix0

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 REPLYlink written 19 months ago by mobiusklein160
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 803 users visited in the last hour