I have a .bam, as ouput from STAR aligner, from which I need to extract some info using RSeQC while using all the computational resources available to increase speed. To accomplish that I am using GNU parallel. Currently I am doing it like this:

  1. Split the main .bam per chromosome, generating files with the suffix Aligned.sortedByCoord.out.REF

    bamtools split -in Aligned.sortedByCoord.out.bam  -reference
  2. Calculate the clipping profile

    ls Aligned.sortedByCoord.out.REF_* | parallel --joblog clipping_profile.log --eta --bar --max-procs 20 --jobs 20 'clipping_profile.py -i {} -s "SE" -o RSeQC.clipping.profile/clipping_profile.{.}'

This will generate a .r, .pdfand a .xlxs for each chromosome, thus ending up with a lot of files. I am forced to split the .bam per chromosome as I am not able to find a --recstart nor a --recend in the compressed .bam.


Is there a better way to split a .bam file using GNU parallel and RSeQC' programs (e.g. split_bam.py), such that one gets the outputs from a single input file?

Inspired by ole.tange's thread, I have also tried this, though unsuccessfully:

cat star_outputAligned.sortedByCoord.out.bam | parallel --block 1G -k --joblog splitbam3.log --max-procs 20 --jobs 20 --pipe 'split_bam.py -i {} -r ../../My_genomes/GCA.GRCh38/GRCh38_wgEncodeGencodeAttrsV28_all_rRNA.bed12 -o splitbam3'

Your Answer

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

Browse other questions tagged or ask your own question.