Unconfigured Ad

Collapse
X
 
  • Time
  • Show
Clear All
new posts
  • IonTom
    Member
    • Apr 2014
    • 32

    #1

    Python samtools BytesIO

    I currently have a a problem with samtools and python.
    Iḿ trying to write a sam or bam file into a file like object
    which works well using the following code

    import sh
    import io

    #define samtools as a python functions
    samtools = sh.Command("samtools")

    #where is the input sam
    sam_input = "test.sam"

    #generate a file like object
    sam_IO2 = io.BytesIO()

    #let samtools view write in the file like object
    samtools.view(sam_output,S=True,h=True,_out = sam_IO2)

    #show the sam file
    sam_IO2.getvalue()

    or

    bam_IO2 = io.BytesIO()
    samtools.view(bam_input,_out = bam_IO2)

    #show the bam file
    bam_IO2.getvalue()


    but now i try to read from the file like object for both cases:

    samtools.view("/dev/stdin",S=True,_in = sam_IO2)

    samtools.view("/dev/stdin",_in = bam_IO2)

    it fails in both tcases.


    Does anyone have an idea how to do this correctly ?
    Would be really helpful.
    Last edited by IonTom; 05-05-2014, 12:28 PM.
  • dpryan
    Devon Ryan
    • Jul 2011
    • 3478

    #2
    Without actually going through your code, is there a reason you're not just using pysam? That would seem simpler.

    Also, it's generally better to give the actual error message produced instead of just writing "it fails".

    Comment

    • IonTom
      Member
      • Apr 2014
      • 32

      #3
      Sorry forgot tom post this.
      Here is the error message:

      samtools.view("/dev/stdin",_in = bam_IO2)

      Traceback (most recent call last):
      File "<stdin>", line 1, in <module>
      File "/usr/local/lib/python2.7/dist-packages/sh.py", line 769, in __call__
      return RunningCommand(cmd, call_args, stdin, stdout, stderr)
      File "/usr/local/lib/python2.7/dist-packages/sh.py", line 330, in __init__
      self.wait()
      File "/usr/local/lib/python2.7/dist-packages/sh.py", line 334, in wait
      self._handle_exit_code(self.process.wait())
      File "/usr/local/lib/python2.7/dist-packages/sh.py", line 348, in _handle_exit_code
      self.process.stderr
      sh.ErrorReturnCode_1:

      RAN: '/usr/bin/samtools view /dev/stdin'

      STDOUT:

      STDERR:
      [bam_header_read] EOF marker is absent. The input is probably truncated.
      [bam_header_read] invalid BAM binary header (this is not a BAM file).
      [main_samview] fail to read the header from "/dev/stdin".

      Isn't pysam more or less doing the same a the sh solution ?
      Basically calling the c implementation ?
      So would it help to interact with the file like objects ?

      Comment

      • dpryan
        Devon Ryan
        • Jul 2011
        • 3478

        #4
        Pysam actually wraps the C API, so it's a bit more useful. Why bother reinventing the wheel.

        Comment

        • maubp
          Peter (Biopython etc)
          • Jul 2009
          • 1544

          #5
          I would also recommend trying pysam instead, but it could be your code is breaking due to a simple user error:

          Where is the BytesIO handle position within the data stream? Does a
          Code:
          .seek(0)
          to ensure it is at the start help? My guess is the handle position is at the end of the data, and thus when samtools tries to read from it there is no (more) data.

          Comment

          • IonTom
            Member
            • Apr 2014
            • 32

            #6
            Here is a solution to the problem above using pipes in python

            bam_input = "input.bam"
            bam_output = "output.bam"

            outputter = open(bam_output,mode = "w")

            p1 = subprocess.Popen(['samtools',"view","-u",bam_output],stdout=subprocess.PIPE,bufsize=1) #Set up the samtools view command to output to pipe
            p2 = subprocess.Popen(["samtools","view","-b","-"], stdin=p1.stdout,stdout=outputter,bufsize = 1) #send p1's output to p2 and p2 to bam file
            p1.stdout.close() #make sure we close the output so p2 doesn't hang waiting for more input
            output = p2.communicate()[0] #run our commands

            outputter.close()

            This example is just a toy example but the principle can also be used for real tasks.
            Last edited by IonTom; 05-13-2014, 03:19 PM.

            Comment

            Latest Articles

            Collapse

            • SEQadmin2
              Multiomic Liquid Biopsy: What DNA/RNA Panels Reveal That DNA Alone Can’t
              by SEQadmin2



              Tumors contain complex mixtures of cell populations that evolve in response to treatment and host-specific factors. Tissue biopsies provide only static snapshots of a tumor’s state. But what if we could follow that evolution over time using a minimally invasive technique?

              Liquid biopsies can detect tumor-derived material circulating in the bloodstream, providing a way to study treatment response, detect emerging resistance, and discover new biomarkers. Circulating
              ...
              10-14-2026, 09:55 AM
            • SEQadmin2
              New Genomics Technologies Take Aim at Long-Standing Limits
              by SEQadmin2


              Researchers using sequencing and genomics tools often have to make trade-offs. They can choose between speed or scale, short reads or long-range information, or targeted panels or a view of the whole transcriptome. New technologies that have been released this year are built to address those tough choices.

              We asked six companies the same four questions to learn about their latest products. The new technologies bring a lot to the table, including rethinking sequencing
              ...
              09-28-2026, 10:25 AM

            ad_right_rmr

            Collapse

            News

            Collapse

            Topics Statistics Last Post
            Started by SEQadmin2, 10-15-2026, 10:52 AM
            0 responses
            26 views
            0 reactions
            Last Post SEQadmin2  
            Started by SEQadmin2, 10-14-2026, 05:39 AM
            0 responses
            33 views
            0 reactions
            Last Post SEQadmin2  
            Started by SEQadmin2, 10-13-2026, 09:56 AM
            0 responses
            31 views
            0 reactions
            Last Post SEQadmin2  
            Started by SEQadmin2, 10-09-2026, 11:14 AM
            0 responses
            28 views
            0 reactions
            Last Post SEQadmin2  
            Working...