Answers.Body: <p>Taking a different tack from other answers, there's lots of tools for pipelines in Python. Note: there was a time when people would use "pipeline" to refer to a shell script. I'm talking about something more sophisticated that helps you decompose an analysis into parts and  runs it robustly.</p>
  
  <ul>
  <li><strong>Snakemake</strong> is my favourite. It's (nearly) pure Python and can generate reports.</li>
  <li><strong>Nextflow</strong> is growing in popularity and is pretty straightforward</li>
  <li><strong>Ruffus</strong> used to be reasonably popular, seemed fine when I used it</li>
  <li><strong>Bpipe</strong> is for bioinformatics</li>
  <li><strong>Airflow</strong> is a more industrial "big" solution but in wide use</li>
  <li>See this <a href="https://github.com/pditommaso/awesome-pipeline" rel="noreferrer">big list of pipeline systems</a></li>
  </ul>
  
  <p>I'm sure you can find something there.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-06T09:12:57.450Z
Answers.LastActivityDate: 2017-07-06T09:12:57.450Z
Answers.UserDisplayName: agapow
Answers.Body: <p>BioPython has some good tools for processing reads and alignments.
  <a href="http://biopython.org/DIST/docs/tutorial/Tutorial.html" rel="nofollow noreferrer">http://biopython.org/DIST/docs/tutorial/Tutorial.html</a></p>
  
  <p>There is a python library wrapping samtools so many of the samtools calls can be used directly as python objects and calls
  <a href="https://pysam.readthedocs.io/en/latest/" rel="nofollow noreferrer">https://pysam.readthedocs.io/en/latest/</a></p>
  
  <p>I would use subprocess to call the aligner and specify the output to a bam file that you have named and then read that bam file with pysam to do the analysis that you are interested in.  </p>
  
  <p>If you provide more specifics about the analysis or aligners I can help you more but you are in the right track just need to connect the two calls together.</p>
  
  <p>example:</p>
  
  <pre><code>import subprocess
  import pysam
  import os 
  
  for root, dirs, files in os.walk(rootPath):
      for filename in fnmatch.filter(files, ".fastq"):
  
          subprocess.run([`hisat`, `-f`, filename `-o`, `some.bam`])
  
          #get depth
          depth = pysam.depth('some.bam')
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-06T02:36:07.097Z
Answers.Text: edited code snippet based on OP comments about executing on multiple fastq files
Answers.UserDisplayName: Bioathlete
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-05T22:30:23.793Z
Answers.LastActivityDate: 2017-07-06T02:35:12.113Z
Answers.LastEditDate: 2017-07-06T02:35:12.113Z
Answers.UserDisplayName: Bioathlete
Answers.Body: <p>If you're just doing alignment and conversion to a sorted BAM file, there's no need to run it through python. A simple pipe on the unix command line works just as well (and probably runs faster):</p>
  
  <pre><code>hisat2 -x genome.index -1 reads_R1.fq.gz -2 reads_R2.fq.gz | 
    samtools sort &gt; reads_vs_genome.bam
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-06T01:38:01.643Z
Answers.Text: thank you @gringer, I am familiar with unix pipe. I am trying to learn to create a pipeline where when I give path of a directory which has fastqs I intend to submit .sh files to cluster and get the desired results . for example `python my_pipeline.py -dir /pathtofastq`  would submit jobs to cluster here i will use bunch of argparse and subprocess and predefinec bsub commands.
Answers.UserDisplayName: novicebioinforesearcher
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-06T02:31:45.050Z
Answers.Text: Look into python's os.path.walk.  This will allow you to take in the dir and then iterate over all the fastq files in the directory.  http://www.pythonforbeginners.com/fnmatch/os-walk-and-fnmatch-in-python
Answers.UserDisplayName: Bioathlete
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-06T00:25:58.437Z
Answers.LastActivityDate: 2017-07-06T02:51:18.897Z
Answers.LastEditDate: 2017-07-06T02:51:18.897Z
Answers.UserDisplayName: gringer
Answers.Body: <p>Have a look at <a href="https://github.com/spotify/luigi" rel="nofollow noreferrer">luigi</a> for writing pipelines in python. You can found BioLuigi if you like :-P</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-06T12:27:44.190Z
Answers.Text: could you please share link for BioLuigi I could not find it online
Answers.UserDisplayName: novicebioinforesearcher
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-07T08:46:04.183Z
Answers.Text: You need to start it ;-)
Answers.UserDisplayName: Dan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-07T13:58:17.883Z
Answers.Text: haha good one!!
Answers.UserDisplayName: novicebioinforesearcher
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-06T09:00:02.707Z
Answers.LastActivityDate: 2017-07-06T09:00:02.707Z
Answers.UserDisplayName: Dan
Answers.Body: <p>I <a href="https://bioinformatics.stackexchange.com/a/1008/191">agree</a> that using a specialized tool is probably a good idea.</p>
  
  <p>Nevertheless if you want to stick with Python, I suggest using <a href="https://plumbum.readthedocs.io/en/latest/" rel="nofollow noreferrer"><code>plumbum</code></a> instead of <code>subprocess</code>. It has a very nice syntax for this kind of problems.</p>
  
  <pre class="lang-python prettyprint-override"><code>from plumbum import local
  
  # load a command
  samtools = local['samtools']
  
  # you can then easily redirect
  (samtools['view', '-bS', '../some.sam'] &gt; '../some.bam')()
  
  # it is also easy to pipe 
  hisat2 = local['hisat2']
  chain = hisat2['-x', 'genome.index', '-1', 'reads_R1.fq.gz', '-2', 'reads_R2.fq.gz'] | samtools['sort'] &gt; 'reads_vs_genome.bam'
  
  chain()
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-06T09:53:12.087Z
Answers.LastActivityDate: 2017-07-06T11:16:07.240Z
Answers.LastEditDate: 2017-07-06T11:16:07.240Z
Answers.UserDisplayName: Iakov Davydov
Body: <p>I was wondering if there is a tutorial or a small code snippet to understand how to write bioinformatics pipeline using python, for example</p>
  
  <ol>
  <li>use a aligner (say hisat) </li>
  <li>get the output and process it using samtools </li>
  </ol>
  
  <p>I was able to use subprocess from python2.7 for this purpose using samtools but i am not able to link both the processes.i.e given path(which I can use argparse) for directory with fastq files output would be processed bam. </p>
  
  <p>sample code for samtools sam to bam :</p>
  
  <pre><code>import subprocess
  subprocess.run(['samtools','view', '-bS',
                  '../some.sam',
                  '&gt;',
                  '../some.bam'])
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-06T06:53:02.207Z
Comments.Text: If you ever need more elaborate pipelines in python, then I would strongly encourage you to look into [snakemake](https://snakemake.readthedocs.io/en/stable/), which is pretty commonly used for this purpose.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-06T11:58:16.583Z
Comments.Text: I agree that snakemake is worth learning. If you want to use subprocess, I found the following tutorial just yesterday, which looked quite clear: http://crashcourse.housegordon.org/python-subprocess.html
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-06T12:22:34.500Z
Comments.Text: @bli this was useful
Comments.UserDisplayName: novicebioinforesearcher
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-06T12:28:39.000Z
Comments.Text: @DevonRyan I looked snakemake it is available only in python3 is that correct?
Comments.UserDisplayName: novicebioinforesearcher
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-06T12:29:38.237Z
Comments.Text: @novicebioinforesearcher Yes, that's correct.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-06T15:02:49.093Z
Comments.Text: @novicebioinforesearcher Note that if you need to program in python2, you can still do it in external scripts and call them from snakemake like you would do for any other program.
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-06T15:31:30.560Z
Comments.Text: @bli how can I do that, sorry I am still learning stuff here could you share an example
Comments.UserDisplayName: novicebioinforesearcher
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-06T16:41:55.857Z
Comments.Text: @novicebioinforesearcher These comments are not adapted to give examples. What I mean is that you will still need python3 for snakemake, but if there are things you can only do in python2, you can write whatever computation you want to do in python2 in the form of an independent program (presumably accepting command-line arguments). Then, you will call this program from within snakemake (or whatever workflow management tool you happen to use, or a python script using subprocess) exactly like you would do for samtools, hisat or any other program.
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-06T18:28:52.330Z
Comments.Text: I'm curious as to why you are averse to writing your pipeline in bash? @novicebioinforesearcher
Comments.UserDisplayName: quantik
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-06T18:44:11.820Z
Comments.Text: @quantik because I am learning to write codes in python and also because (i may be wrong here) in case you have a greater number of samples bash would not be a good way to work(foreg 40 samples to be processed) i can use os.walk etc in python to submit jobs correct me if i am wrong
Comments.UserDisplayName: novicebioinforesearcher
Tags: <rna-seq><samtools><python><genomics>
Title: using python to write bioinformatics pipelines tutorial
ThreadText: # using python to write bioinformatics pipelines tutorial
  
  (Asked by: username_0 on 2017-07-05)
  
  I was wondering if there is a tutorial or a small code snippet to understand how to write bioinformatics pipeline using python, for example
  
  1. use a aligner (say hisat)
  2. get the output and process it using samtools
  
  I was able to use subprocess from python2.7 for this purpose using samtools but i am not able to link both the processes.i.e given path(which I can use argparse) for directory with fastq files output would be processed bam.
  
  sample code for samtools sam to bam :
  
  ```
  import subprocess
  subprocess.run(['samtools','view', '-bS',
                  '../some.sam',
                  '>',
                  '../some.bam'])
  
  ```
  
  **Tags:** rna-seq, samtools, python, genomics
  
  > **Comments:**
  > **username_1**: If you ever need more elaborate pipelines in python, then I would strongly encourage you to look into [snakemake](https://snakemake.readthedocs.io/en/stable/), which is pretty commonly used for this purpose.  
  > **username_2**: I agree that snakemake is worth learning. If you want to use subprocess, I found the following tutorial just yesterday, which looked quite clear: http://crashcourse.housegordon.org/python-subprocess.html  
  > **username_0**: @username_2 this was useful  
  > **username_0**: @username_1 I looked snakemake it is available only in python3 is that correct?  
  > **username_1**: @username_0 Yes, that's correct.  
  > **username_2**: @username_0 Note that if you need to program in python2, you can still do it in external scripts and call them from snakemake like you would do for any other program.  
  > **username_0**: @username_2 how can I do that, sorry I am still learning stuff here could you share an example  
  > **username_2**: @username_0 These comments are not adapted to give examples. What I mean is that you will still need python3 for snakemake, but if there are things you can only do in python2, you can write whatever computation you want to do in python2 in the form of an independent program (presumably accepting command-line arguments). Then, you will call this program from within snakemake (or whatever workflow management tool you happen to use, or a python script using subprocess) exactly like you would do for samtools, hisat or any other program.  
  > **username_3**: I'm curious as to why you are averse to writing your pipeline in bash? @username_0  
  > **username_0**: @username_3 because I am learning to write codes in python and also because (i may be wrong here) in case you have a greater number of samples bash would not be a good way to work(foreg 40 samples to be processed) i can use os.walk etc in python to submit jobs correct me if i am wrong  
  
  ---
  
  ## 5 answers
  
  ### Answer by username_8 on 2017-07-06. Score: 4 | Accepted answer
  
  I [agree](https://bioinformatics.stackexchange.com/a/1008/191) that using a specialized tool is probably a good idea.
  
  Nevertheless if you want to stick with Python, I suggest using [`plumbum`](https://plumbum.readthedocs.io/en/latest/) instead of `subprocess`. It has a very nice syntax for this kind of problems.
  
  ```python
  from plumbum import local
  
  # load a command
  samtools = local['samtools']
  
  # you can then easily redirect
  (samtools['view', '-bS', '../some.sam'] > '../some.bam')()
  
  # it is also easy to pipe 
  hisat2 = local['hisat2']
  chain = hisat2['-x', 'genome.index', '-1', 'reads_R1.fq.gz', '-2', 'reads_R2.fq.gz'] | samtools['sort'] > 'reads_vs_genome.bam'
  
  chain()
  
  ```
  
  
  ---
  
  ### Answer by username_4 on 2017-07-06. Score: 13
  
  Taking a different tack from other answers, there's lots of tools for pipelines in Python. Note: there was a time when people would use "pipeline" to refer to a shell script. I'm talking about something more sophisticated that helps you decompose an analysis into parts and runs it robustly.
  
  * **Snakemake** is my favourite. It's (nearly) pure Python and can generate reports.
  * **Nextflow** is growing in popularity and is pretty straightforward
  * **Ruffus** used to be reasonably popular, seemed fine when I used it
  * **Bpipe** is for bioinformatics
  * **Airflow** is a more industrial "big" solution but in wide use
  * See this [big list of pipeline systems](https://github.com/pditommaso/awesome-pipeline)
  
  I'm sure you can find something there.
  
  
  ---
  
  ### Answer by username_5 on 2017-07-05. Score: 5
  
  BioPython has some good tools for processing reads and alignments.
  <http://biopython.org/DIST/docs/tutorial/Tutorial.html>
  
  There is a python library wrapping samtools so many of the samtools calls can be used directly as python objects and calls
  <https://pysam.readthedocs.io/en/latest/>
  
  I would use subprocess to call the aligner and specify the output to a bam file that you have named and then read that bam file with pysam to do the analysis that you are interested in.
  
  If you provide more specifics about the analysis or aligners I can help you more but you are in the right track just need to connect the two calls together.
  
  example:
  
  ```
  import subprocess
  import pysam
  import os 
  
  for root, dirs, files in os.walk(rootPath):
      for filename in fnmatch.filter(files, ".fastq"):
  
          subprocess.run([`hisat`, `-f`, filename `-o`, `some.bam`])
  
          #get depth
          depth = pysam.depth('some.bam')
  
  ```
  
  > **Comments:**
  > **username_5**: edited code snippet based on OP comments about executing on multiple fastq files  
  
  ---
  
  ### Answer by username_6 on 2017-07-06. Score: 4
  
  If you're just doing alignment and conversion to a sorted BAM file, there's no need to run it through python. A simple pipe on the unix command line works just as well (and probably runs faster):
  
  ```
  hisat2 -x genome.index -1 reads_R1.fq.gz -2 reads_R2.fq.gz | 
    samtools sort > reads_vs_genome.bam
  
  ```
  
  > **Comments:**
  > **username_0**: thank you @username_6, I am familiar with unix pipe. I am trying to learn to create a pipeline where when I give path of a directory which has fastqs I intend to submit .sh files to cluster and get the desired results . for example `python my_pipeline.py -dir /pathtofastq`  would submit jobs to cluster here i will use bunch of argparse and subprocess and predefinec bsub commands.  
  > **username_5**: Look into python's os.path.walk.  This will allow you to take in the dir and then iterate over all the fastq files in the directory.  http://www.pythonforbeginners.com/fnmatch/os-walk-and-fnmatch-in-python  
  
  ---
  
  ### Answer by username_7 on 2017-07-06. Score: 0
  
  Have a look at [luigi](https://github.com/spotify/luigi) for writing pipelines in python. You can found BioLuigi if you like :-P
  
  > **Comments:**
  > **username_0**: could you please share link for BioLuigi I could not find it online  
  > **username_7**: You need to start it ;-)  
  > **username_0**: haha good one!!  
  
  ---
  
  