sc932 / ALE

Assembly Likelihood Estimator
Other
32 stars 7 forks source link

Plotting a single contig #6

Closed sylvainforet closed 10 years ago

sylvainforet commented 11 years ago

Thanks for developing ALE, very useful! Here is a trivial patch that make the plotting of single contigs considerably faster:

diff --git a/src/plotter3.py b/src/plotter3.py
index c1674e1..20ad276 100755
--- a/src/plotter3.py
+++ b/src/plotter3.py
@@ -773,7 +773,7 @@ def plot_histogram(input_data, save_figure=False, pdf_stream=None):
     else:
         plt.show()

-def read_in_info(placement_file):
+def read_in_info(placement_file, contig=None):
     """Reads in an ALE placement file, returns a list of Contigs.

     Args:
@@ -824,16 +824,22 @@ def read_in_info(placement_file):
     previous_line_one = ""
     previous_line_two = ""

+    hasContig=False
+
     for line in ale_placement_file:
         if line[0] == '#':
+            if hasContig:
+                break
             if previous_line_one == "":
                 previous_line_one = line
             else:
                 previous_line_two = previous_line_one
-                previous_line_one = line           
+                previous_line_one = line
         else:
             if previous_line_two != "":
                 tName = previous_line_two.split(' ')[2]              
+                if contig and tName == contig:
+                    hasContig = True
                 tLen = int(previous_line_two.split(' ')[-2])
                 contigs.append(Contig(tLen, name = tName))
                 place = 0
@@ -841,6 +847,8 @@ def read_in_info(placement_file):
                 print ""
                 bar = ProgressBar.progressBar(0, tLen, 42)
                 previous_line_two = ""
+            if contig and not hasContig:
+                continue
             data = line.split(' ')
             contigs[-1].depth[place] = numpy.double(data[2])
             for i in range(1,7):
@@ -925,7 +933,7 @@ def main():
     user_params.read_sys_args()

     # read in contigs
-    contigs = read_in_info(user_params.get_input("ale_file"))
+    contigs = read_in_info(user_params.get_input("ale_file"), user_params.get("specific_contig"))

     save_figure = user_params.get("save_figure")
     # open up a pdf_stream and file for output
sc932 commented 10 years ago

Thanks! I'll apply this patch in the next version release. In the meantime...

The authors recommend using IGV to view the output.

http://www.broadinstitute.org/igv/

Just import the assembly, bam and ALE scores. You can convert the .ale file to a set of .wig files with ale2wiggle.py and IGV can read those directly. Depending on your genome size you may want to convert the .wig files to the BigWig format.

The internal, debugging plotter is not in active development and is no longer supported.