6
$\begingroup$

I have an obitools script (de Barba et al. 2016) that I would like to run faster. How would you run it in parallel to cut down on time?

illuminapairedend -r rawdata_scandinavia_R2.fastq rawdata_scandinavia_R1.fastq | tee rawdata_scandinavia.fastq | obiannotate -S goodAli:'"Alignement" if score>40.00 else "Bad"' | obisplit -t goodAli -p rawdata_scandinavia.
touch rawdata_scandinavia.Bad.fastq
touch rawdata_scandinavia.Alignement.fastq
touch rawdata_scandinavia.fastq
ngsfilter -t rawdata_scandinavia.ngsfilter -u rawdata_scandinavia.unidentified.fastq rawdata_scandinavia.Alignement.fastq > rawdata_scandinavia.filtered.fastq
obisplit -p MICROSAT.PCR_ -t experiment rawdata_scandinavia.filtered.fastq

Can someone describe what tools and steps one needs to take to run this job in parallel, given that the functions are not designed to be run in parallel out of the box?

De Barba, M., Miquel, C., Lobréaux, S., Quenette, P. Y., Swenson, J. E., & Taberlet, P. (2016). High-throughput microsatellite genotyping in ecology: improved accuracy, efficiency, standardization and success with low-quantity and degraded DNA. Molecular Ecology Resources, 1–16. https://doi.org/10.1111/1755-0998.12594

$\endgroup$
  • 2
    $\begingroup$ I think this is a more general programming question and so would be more suited for stack overflow $\endgroup$ – Chris_Rands May 19 '17 at 14:02
  • 1
    $\begingroup$ Divide & conquer! $\endgroup$ – Kamil S Jaron May 19 '17 at 14:44
  • 4
    $\begingroup$ Bioinformatics encompasses biology, statistics, and computer science. Parallel processing is a common need in bioinformatics, and can be a challenging task to implement depending on the particular programs that are used. The programs used here in the example script are specific to bioinformatics, and answers in stack overflow might be lower quality due to the inexperience of the community in working with bioinformatics programs. $\endgroup$ – gringer Jun 1 '17 at 22:12
  • 2
    $\begingroup$ Pure programming questions (as was the impression from reading only the question title for this) are better placed in SO. However, this question is about a specific type of parallel processing that is associated with a particular bioinformatic tool. I can imagine that this question might be closed in SO for offtopic reasons. $\endgroup$ – gringer Jun 2 '17 at 8:25
  • 3
    $\begingroup$ @gringer I hope we can get one more vote. I've solved the problem in the mean time and would like to offer my fruits of labor to anybody who might come up with this kind of a problem in the future. $\endgroup$ – Roman Luštrik Jun 2 '17 at 8:30
3
$\begingroup$

Your only option here is to first split the fastq files into small chunks, process each and then merge them back together. Alternatively, find something else that runs faster.

$\endgroup$
3
$\begingroup$

Here is a bash script and GNU parallel along with some other tools to split-apply-combine the data. It may or may not work out of the box, but it's not far off and it should give the general idea of how one way of how to approach this.

#!/bin/bash

# ====== USER PROVIDED PARAMETERS ========
N=20 # number of cores to use
ngsfilter="rawdata_scandinavia.ngsfilter"
unidentified="unidentified_samples.fastq"

R1="rawdata_scandinavia_R1.fastq"
R2="rawdata_scandinavia_R2.fastq"
# ==== END USER PROVIDED PARAMETERS =====

# fastqutils available @ http://ngsutils.org/modules/fastqutils
fastqutils split $R1 to_illumina_R1_ $N &
fastqutils split $R2 to_illumina_R2_ $N &
wait

R1=$( ls | grep -P "to_illumina_R1_" )
R2=$( ls | grep -P "to_illumina_R2_" )

# Do alignment of two reads.
for i in $( seq $N )
do
    partR1=$( ls | grep -P "^to_illumina_R1_\\.$i\\.fastq$" )
    partR2=$( ls | grep -P "^to_illumina_R2_\\.$i\\.fastq$" )
    illuminapairedend -r $partR1 $partR2 | tee to_ngsfilter_$i.fastq | obiannotate -S goodAli:'"Align" if score>40.00 else "Bad"' | obisplit -t goodAli -p to_ngsfilter_$i. &
done

# Remove intermediate aligned files.
echo $R1 | xargs rm -r
echo $R2 | xargs rm -r
ls | grep -P "\\.seq$|\\.err$" | xargs rm -r

input=$( ls | grep -P "to_ngsfilter_[0-9]+\\.Align\\.fastq" )

# For some reasons, files need to be "touched" in order to work. Perhaps
# it's not closing them...
touch $unidentified

touchedbyangel=$( ls | grep -P "to_ngsfilter_"  )
for i in $touchedbyangel
do
    touch $i
done

parallel -j$N --result filtered_{} ngsfilter -t $ngsfilter -u $unidentified {} ::: $input

# Remove intermediate results.
ls | grep -P "\\.seq$|\\.err$" | xargs rm -r

ls | grep -P "filtered_.*\\.fastq$" | xargs cat > filtered_data.fastq

# remove intermediate results to conserve space
rm to_ngsfilter*

# Split reads by locus.
filtered=$( ls | grep -P "^filtered_"  )
parallel -j$N --result to_split_{} obisplit -p MS.PCR_ -t experiment {} ::: $filtered

ls | grep -P "\\.err$|\\.seq$" | xargs rm -r
ls | grep -P "^to_split_" | xargs rm -r

# Wasn't able to parallelize this, so it's just pushing tasks into the background.
obiuniq -m sample MS.PCR_UA_MxRout1_03.fastq > MS.PCR_UA_MxRout1_03.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_06.fastq > MS.PCR_UA_MxRout1_06.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_14.fastq > MS.PCR_UA_MxRout1_14.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_16.fastq > MS.PCR_UA_MxRout1_16.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_17.fastq > MS.PCR_UA_MxRout1_17.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_25.fastq > MS.PCR_UA_MxRout1_25.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_51.fastq > MS.PCR_UA_MxRout1_51.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_57.fastq > MS.PCR_UA_MxRout1_57.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_63.fastq > MS.PCR_UA_MxRout1_63.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_64.fastq > MS.PCR_UA_MxRout1_64.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_65.fastq > MS.PCR_UA_MxRout1_65.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_67.fastq > MS.PCR_UA_MxRout1_67.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_68.fastq > MS.PCR_UA_MxRout1_68.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_ZF.fastq > MS.PCR_UA_MxRout1_ZF.uniq.fasta &
wait

uniqfiles=$( ls | grep -P ".{2}+\\.uniq\\.fasta$" )

for i in $uniqfiles
do
    # this section removes the file suffix
    # https://stackoverflow.com/questions/125281/how-do-i-remove-the-file-suffix-and-path-portion-from-a-path-string-in-bash
    tab=${i%.fasta}
    tab=${tab##*/}
    obigrep -p 'count>1' $i | obiannotate -k merged_sample -k count | obiannotate --length | obisort -r -k seq_length | obitab --no-definition --output-seq > $tab.tab
done

ls | grep -P "\\.uniq\\.fasta" | xargs rm -r

echo "Done. I think."
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.