Extract RP transcripts for provided organism about htsinfer HOT 15 CLOSED

Hey, I'd like to work on this with Daniel Batora.

from htsinfer.

Hey, I'd like to take this issue.

from htsinfer.

Sure, @rockdeme & @danielbatora

from htsinfer.

This is the code that generated the error :

from Bio import SeqIO

reads = []
with open("transcripts.fasta", "rU") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        reads.append(record.description.split("|")[3])

from htsinfer.

@danielbatora: unfortunately I couldn't reproduce the error. However, the file obviously has some issues. Most of the times it has Unix-style line breaks, but sometimes Windows-style line breaks. And there's a really weird block of characters (looks like binary content) starting in line 166572. I have to get back to the person who created this file, but for now I have removed the "binary" characters and changed all line breaks to Unix style. As I couldn't reproduce the error, I can't verify that the new version of the file is okay. Let's try this together tomorrow.

from htsinfer.

@uniqueg I pushed the pilot function to the fork repo. As of now, the function accepts SeqRecord objects, so it does not run into error with the corrupted file. When reading and cleaning the file happened inside the function (in that case the function acceps a path to the FASTA file as an argument, the error did not let the function to return. Other than that, which file type should I use when creating the temporary directory ? Is pickle good?

from htsinfer.

Thanks @danielbatora: can you try with the updated file? I think it's best if the function takes and returns a FASTA file path, otherwise we would need an extra function to read it and just outsource the problem :)

from htsinfer.

@uniqueg Okay, should the output file be a FASTA as well, or pickle ?

from htsinfer.

Yup, FASTA. The STAR alignment program used in #26 and #27 is not a Python package and thus cannot deal with a SeqRecord object or a Pickle file.

from htsinfer.

Okay, Done ! Should the temporary directory created in another issue be an argument of this function ? I do not understand it completely.

from htsinfer.

Yes, good point, the temporary directory needs to be passed as an argument (e.g., tmp_dir). I will implement the code to create the temp directory, you can add the parameter to the function signature (remember to add a type hint and a description in the docstring) and then create your output file as something like:

import os

out_file = os.path.join(tmp_dir, "transcripts_filtered.fasta")

(untested)

from htsinfer.

As discussed in the chat:

if (
    isinstance(organism, str) and 
    record.description.split("|")[3] == organism
):
    data.update({organism + "_" + str(counter) : str(record.seq)})
    counter += 1 
elif (
    isinstance(organism, int) and
    record.description.split("|")[4] == str(organism)
):
    data.update({str(organism) + "_" + str(counter) : str(record.seq)})
    counter += 1 

When merging all conditions into a single if this would become something like:

if (
    # check for organism short name
    (
        isinstance(organism, str) and 
        record.description.split("|")[3] == organism
    ) or
    # check for taxon identifier
    (
        isinstance(organism, int) and
        record.description.split("|")[4] == str(organism)
    )
):
    data.update({str(organism) + "_" + str(counter) : str(record.seq)})
    counter += 1 

Note that this last solution, while having a relatively hard to read conditional section (which is why I added comments), avoids code duplication, which is always a source of error (note that the data.update({str(organism) + "_" + str(counter) : str(record.seq)}) call with the extra str() should work for both int and str organism arguments)

from htsinfer.

Hi @danielbatora & @rockdeme , actually forgot to insist on unit tests for the function. But seeing that the file tests/test_infer_read_orientation.py where these tests should go is still empty and that it's rather hard to start on an empty file, and there are currently only unit tests for the wrapper function infer() in tests/test_infer_single_paired.py either, but not tests for each individual function, as examples, I think I can take care of this for this PR. Benefit of the first PR ;-)

That is, unless you would like to try it yourself. In that case, please let me know. In general, if you make changes to your feature branch and push them, the PR will be updated automatically, so it's easy to add tests (or any other changes) to the existing PR.

from htsinfer.

Hi @uniqueg , are we having one more hackaton lecture? In that case, we could do the unit tests or start a new issue (if there is any that noone is assigned to ). What do you suggest ?

from htsinfer.

Yes, we will have one more, so feel free to do it then. I will see that I'll add some more examples for the htsinfer/infer_single_paired.py module until then. I'd prefer to keep the code and corresponding tests in one PR

from htsinfer.