Tutorial:How To Create A Bioinformatics Pipeline Using Spotify !

Question

[Rad]

A lot of big streaming companies like Spotify, Pandora, etc are more and more pushing towards a better and more stable frameworks and the best thing to do that is to go open source and get useful feedback and keep working on making the platform better.

Spotify developed a platform called Luigi, a python framework to handle users logs, and mine them intuitively by plugging several machine learning algorithms to improve their recommendation systems and their suggestions to clients.

Luigi works almost like any make-like python framework for pipeline development, like Ruffus or Snakemake etc.., but it has a plus over these solutions, it is designed to create Hadoop friendly pipelines and also comes with a visual diagnostic of each part of your pipeline while it is running. Another feature I like is that it notifies you via email when a task fails.

Here is a simple adaptation of Luigi for Bioinformatics. This pipeline:

• takes a fastq samples list,
• align them via bwa_mem
• Convert SAM files to BAM files
• Sort the bam files
• Index them
• Call variants using samtools mpileup
• Convert bcf to vcf

Code is viewable and editable here http://coderscrowd.com/app/public/codes/view/229

[pinn]

why its showing error ? how to sort out.

CMD

python Luigipipeline.py --local-scheduler Custom_Genome_Pipeline

====== Running BWA ======
INFO: [pid 18146] Worker Worker(salt=235547989, workers=1, host=Curium, username=likith_reddy, pid=18146) done Bwa_Mem(sample=SRR098401testsickle_1.fastq)
DEBUG: 1 running tasks, waiting for next task to finish
INFO: Informed scheduler that task Bwa_Mem_SRR098401testsic_c845cfe3a6 has status DONE
DEBUG: Asking scheduler for work...
DEBUG: Pending tasks: 11
INFO: [pid 18146] Worker Worker(salt=235547989, workers=1, host=Curium, username=likith_reddy, pid=18146) running Convert_Sam_Bam(sample=SRR098401testsickle_1.fastq)
ERROR: [pid 18146] Worker Worker(salt=235547989, workers=1, host=Curium, username=likith_reddy, pid=18146) failed Convert_Sam_Bam(sample=SRR098401testsickle_1.fastq)
Traceback (most recent call last):
File "/home/likith_reddy/.local/lib/python2.7/site-packages/luigi/worker.py", line 194, in run
new_deps = self._run_get_new_deps()
File "/home/likith_reddy/.local/lib/python2.7/site-packages/luigi/worker.py", line 131, in _run_get_new_deps
task_gen = self.task.run()
File "Luigipipeline.py", line 60, in run
"sam/"+self.sample+".sam"])
File "Luigipipeline.py", line 15, in run_cmd
p = subprocess.Popen(cmd, shell=False, universal_newlines=True, stdout=subprocess.PIPE)
File "/usr/lib/python2.7/subprocess.py", line 711, in init
errread, errwrite)
File "/usr/lib/python2.7/subprocess.py", line 1343, in _execute_child
raise child_exception
OSError: [Errno 2] No such file or directory
DEBUG: 1 running tasks, waiting for next task to finish
INFO: Informed scheduler that task Convert_Sam_Bam_SRR098401testsic_c845cfe3a6 has status FAILED
DEBUG: Asking scheduler for work...
DEBUG: Done
DEBUG: There are no more tasks to run at this time
DEBUG: There are 11 pending tasks possibly being run by other workers
DEBUG: There are 11 pending tasks unique to this worker
DEBUG: There are 11 pending tasks last scheduled by this worker
INFO: Worker Worker(salt=235547989, workers=1, host=Curium, username=likith_reddy, pid=18146) was stopped. Shutting down Keep-Alive thread
INFO:
===== Luigi Execution Summary =====

Scheduled 13 tasks of which:

    2 ran successfully:
        2 Bwa_Mem(sample=SRR098401testsickle_1.fastq,SRR098401testsickle_2.fastq)
    2 failed:
        2 Convert_Sam_Bam(sample=SRR098401testsickle_1.fastq,SRR098401testsickle_2.fastq)
    9 were left pending, among these:
        9 had failed dependencies:
            2 Call_Variant(sample=SRR098401testsickle_1.fastq,SRR098401testsickle_2.fastq)
            2 Convert_Bcf_Vcf(sample=SRR098401testsickle_1.fastq,SRR098401testsickle_2.fastq)
            1 Custom_Genome_Pipeline()
            2 Index_Bam(sample=SRR098401testsickle_1.fastq,SRR098401testsickle_2.fastq)
            2 Sort_Bam(sample=SRR098401testsickle_1.fastq,SRR098401testsickle_2.fastq)

This progress looks :( because there were failed tasks

===== Luigi Execution Summary =====

[Ram]

Read: http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002202

Your question makes no sense as you provide neither context nor steps to reproduce what you tried. If you do not invest effort into making it easier for people that volunteer their time to help you, it will be immensely difficult for you to get help in a timely manner.

[pinn]

Hi, can you provide the Luigi original Code. thanks!

Answer 1

[valentine]

Also see Ratatosk, a pipeline framework for bioinformatics tasks built on Luigi: http://ratatosk.readthedocs.org/en/latest/index.html

[Rad]

Great ! Didn't know about that, thanks for sharing

Answer 2

[Samuel Lampa]

For some bioinformatics tasks, the extra boilerplate needed with luigi can be a bit hampering. This lead me to write a little (~50 lines of code) helper library that lets you define inputs and outputs inline in the command pattern, and connect inputs and outputs between tasks using single-assignment syntax, similar to when you set variables.

So for anyone interested in checking it out, the library is available at github:

• https://github.com/samuell/luigis_monkey_wrench

... and a somewhat realistic NGS bioinformatics example can be found in this gist:

• 	import luigi
  	from luigis_monkey_wrench import *
  	REF='human_17_v37.fasta'
  	INDIVIDUALS=['NA06984','NA07000']
  	SAMPLES=['1','2']
  	BASENAME='.ILLUMINA.low_coverage.17q_'
  	PICARDDIR='/sw/apps/bioinfo/picard/1.69/kalkyl/'
  	KNOWNSITES='/proj/g2014207/labs/gatk/ALL.chr17.phase1_integrated_calls.20101123.snps_indels_svs.genotypes.vcf'
  	class GATKWorkflow(WorkflowTask):
  	def requires(self):
  	for i in INDIVIDUALS:
  	# Workflow definition
  	# ---------------------------------------------------------------------------------
  	# files() will return a pseudo task that just outputs an existing file,
  	# while not running anything.
  	# shell() will create a new task with a command that can take inputs
  	# and outputs.
  	fq1 = file('fastq:{i}/{i}{b}1.fq'.format(i=i,b=BASENAME))
  	fq2 = file('fastq:{i}/{i}{b}2.fq'.format(i=i,b=BASENAME))
  	# Step 2 in [1]--------------------------------------------------------------------
  	aln1 = shell('bwa aln {ref} <i:fastq> > <o:sai:<i:fastq>.sai>'.format(ref=REF))
  	aln1.inports['fastq'] = fq1.outport('fastq')
  	aln2 = shell('bwa aln {ref} <i:fastq> > <o:sai:<i:fastq>.sai>'.format(ref=REF))
  	aln2.inports['fastq'] = fq2.outport('fastq')
  	# Step 3 in [1]--------------------------------------------------------------------
  	merg = shell('bwa sampe {ref} <i:sai1> <i:sai2> <i:fq1> <i:fq2> > <o:merged:{i}/{i}{b}.merged.sam>'.format(
  	ref=REF,
  	i=i,
  	b=BASENAME))
  	merg.inports['sai1'] = aln1.outport('sai')
  	merg.inports['sai2'] = aln2.outport('sai')
  	merg.inports['fq1'] = fq1.outport('fastq')
  	merg.inports['fq2'] = fq2.outport('fastq')
  	# Step 4a in [1]------------------------------------------------------------------
  	mergbam = shell('''
  	java -Xmx2g -jar {p}/AddOrReplaceReadGroups.jar
  	INPUT=<i:sam>
  	OUTPUT=<o:bam:<i:sam>.bam>
  	SORT_ORDER=coordinate
  	RGID={sample}-id
  	RGLB={sample}-lib
  	RGPL=ILLUMINA
  	RGPU={sample}-01
  	RGSM={sample}
  	'''.format(
  	p=PICARDDIR,
  	sample=i))
  	mergbam.inports['sam'] = merg.outport('merged')
  	# Step 4b in [1] -----------------------------------------------------------------
  	index_mergbam = shell('''
  	java -Xmx2g -jar
  	/sw/apps/bioinfo/picard/1.69/kalkyl/BuildBamIndex.jar
  	INPUT=<i:bamr
  	# <o:bai:<i:bam:.bam|.bai>>
  	''')
  	index_mergbam.inports['bam'] = mergbam.outport('bam')
  	# Step 5a in [1]------------------------------------------------------------------
  	local_realign = shell('''
  	java -Xmx2g -jar /sw/apps/bioinfo/GATK/1.5.21/GenomeAnalysisTK.jar
  	-I <i:bam>
  	-R {ref}
  	-T RealignerTargetCreator
  	-o <o:intervals:<i:bam>.intervals>
  	'''.format(
  	ref=REF))
  	local_realign.inports['bam'] = mergbam.outport('bam')
  	# Step 5b in [1]-----------------------------------------------------------------
  	actual_realign = shell('''
  	java -Xmx2g -jar /sw/apps/bioinfo/GATK/1.5.21/GenomeAnalysisTK.jar
  	-I <i:bam>
  	-R {ref}
  	-T IndelRealigner
  	-o <o:realigned_bam:<i:intervals>.realign.bam>
  	-targetIntervals <i:intervals>
  	# <o:realigned_bai:<i:intervals>.realign.bai>
  	'''.format(
  	ref=REF))
  	actual_realign.inports['bam'] = mergbam.outport('bam')
  	actual_realign.inports['intervals'] = local_realign.outport('intervals')
  	# Step 5c in [1]-----------------------------------------------------------------
  	mark_dupes = shell('''
  	java -Xmx2g -jar /sw/apps/bioinfo/picard/1.69/kalkyl/MarkDuplicates.jar '
  	INPUT=<i:bam> '
  	OUTPUT=<o:marked_bam:<i:bam>.marked.bam> '
  	METRICS_FILE=<o:metrics:<i:bam>.marked.metrics>
  	''')
  	mark_dupes.inports['bam'] = actual_realign.outport('realigned_bam')
  	# Step 5d in [1], Index bam (picard does not do that automatically)---------------
  	index_marked_bam = shell(('java -Xmx2g -jar /sw/apps/bioinfo/picard/1.69/kalkyl/BuildBamIndex.jar '
  	'INPUT=<i:bam> '
  	'# <o:bai:<i:bam:.bam|.bai>>'))
  	index_marked_bam.inports['bam'] = mark_dupes.outport('marked_bam')
  	# Step 5e in [1], quality recalibration with GATK---------------------------------
  	count_covar = shell('''
  	java -Xmx2g -jar /sw/apps/bioinfo/GATK/1.5.21/GenomeAnalysisTK.jar
  	-T CountCovariates -I <i:bam>
  	-R {ref}
  	-knownSites {sites}
  	-cov ReadGroupCovariate
  	-cov CycleCovariate
  	-cov DinucCovariate
  	-cov QualityScoreCovariate
  	-recalFile <o:covariate:<i:bam>.covar>
  	'''.format(
  	ref=REF,
  	sites=KNOWNSITES))
  	count_covar.inports['bam'] = mark_dupes.outport('marked_bam')
  	# Step 5f in [1]-----------------------------------------------------------------
  	table_recalib = shell('''
  	java -Xmx2g -jar /sw/apps/bioinfo/GATK/1.5.21/GenomeAnalysisTK.jar
  	-T TableRecalibration
  	-I <i:bam>
  	-R {ref}
  	-recalFile <i:calib>
  	-o <o:calibrated_bam:<i:calib>.bam>
  	'''.format(ref=REF))
  	table_recalib.inports['bam'] = mark_dupes.outport('marked_bam')
  	table_recalib.inports['calib'] = count_covar.outport('covariate')
  	yield table_recalib
  	if __name__ == '__main__':
  	luigi.run()
  	# REFERENCES
  	# ----------
  	# [1] http://uppnex.se/twiki/do/view/Courses/NgsIntro1502/ResequencingAnalysis

  view raw reseqtut_lmw.py hosted with ❤ by GitHub

And, it has a page here on BioStars too:

• Luigi's Monkey Wrench - A small helper library for commandline heavy bioinformatics workflows in Spotify's Luigi workflow tool

[Samuel Lampa]

Just for reference: Luigi's Monkey Wrench is nowadays deprecated in favor of SciLuigi

Answer 3

[ostrokach]

The best pipelining tool I have used is snakemake. Make can work OK as well (+qmake if you are using SGE).

For more info, see here: Workflow management software for pipeline development in NGS