Run Samtools on multiple files

Question: Run Samtools on multiple files

3.6 years ago by

V • 230

UK/London

V • 230 wrote:

I'm trying to run samtools sort in a folder containing about 380 bam files.

I've tried the following command and nothing happens. No errors, it just returns blank. Can anyone help me with what I've done wrong? Thank you :)

samtools sort -@ 7 *.bam > *.sorted.bam

samtools rnaseq • 8.2k views

ADD COMMENT • link •

modified 3.6 years ago by WouterDeCoster ♦ 45k • written 3.6 years ago by V • 230

3.6 years ago by

Medhat ♦ 8.8k

Texas

Medhat ♦ 8.8k wrote:

Try this;

for f in *.bam; do filename="${f%%.*}"; samtools sort -@ 7 $f ${filename}.sorted.bam; done

for f in *.bam; do  samtools sort -@ 7 $f   ${f/.bam/sorted.bam}; done

ADD COMMENT • link modified 3.6 years ago • written 3.6 years ago by Medhat ♦ 8.8k

Seems to be working, hasn't thrown out an error! thanks alot!!

ADD REPLY • link written 3.6 years ago by V • 230

I modified a bit based on Medhat's script. Seems better. for f in *.bam; do samtools sort -@ 7 $f ${f/bam/sorted}; done

ADD REPLY • link written 2.0 years ago by yolandaly • 0

3.6 years ago by

WouterDeCoster ♦ 45k

Belgium

WouterDeCoster ♦ 45k wrote:

Using gnu parallel in its simplest form:

ls *.bam | parallel -j 8 'samtools sort -@ 7  {} > {.}_sorted.bam'

To show progress and expected time of completion

ls *.bam | parallel --progress --eta -j 8 'samtools sort -@ 7  {} > {.}_sorted.bam'

ADD COMMENT • link written 3.6 years ago by WouterDeCoster ♦ 45k

--eta implies --progress, so --progress is not needed. Also try the newer --bar. If you want one per cpu core, leave out -j8.

ADD REPLY • link modified 3.6 years ago • written 3.6 years ago by ole.tange • 4.0k

Thanks! I am huge fan of parallel, but so many options :o

ADD REPLY • link written 3.6 years ago by WouterDeCoster ♦ 45k

3.6 years ago by

Pierre Lindenbaum ♦ 133k

France/Nantes/Institut du Thorax - INSERM UMR1087

Pierre Lindenbaum ♦ 133k wrote:

create the following 'Makefile'

SHELL=/bin/bash
BAMS=$(shell ls -1 *.bam)
define run

$$(addsuffix .sorted.bam,$$(basename $(1))): $(1)
    samtools sort $$< $$(basename $$@)
    samtools index $$@
endef

all: $(addsuffix .sorted.bam,$(basename ${BAMS}))

$(eval $(foreach B,${BAMS},$(call run,${B})))

and run it in parallel for you 380 bam using the option -j for the number of parallel operations

make -j 16

ADD COMMENT • link written 3.6 years ago by Pierre Lindenbaum ♦ 133k

On a separate note, this is really cool, suggestions to learn how to use "make" (tutorial links) would be helpful for all.

ADD REPLY • link written 3.6 years ago by lakhujanivijay ♦ 5.4k

May be I am asking silly question, does that "16" means the no. of threads?

ADD REPLY • link written 3.6 years ago by lakhujanivijay ♦ 5.4k

does that "16" means the no. of threads?

yes

ADD REPLY • link written 3.6 years ago by Pierre Lindenbaum ♦ 133k

3.6 years ago by

Charles Plessy • 2.7k

Japan

Charles Plessy • 2.7k wrote:

The * wildcard character can only expand to already existing files (*.bam in your example). The command-line interpreter can not guess your intent and provide contents for *.sorted.bam. Moreover, samtools sort is not designed to run on multiple files in parallel.

However, there are generic tools for doing what you want. See for instance the two following links.

ADD COMMENT • link written 3.6 years ago by Charles Plessy • 2.7k

Hi Charles,

Sorry My bad, I seem to have explained what I want to do wrong. I didn't mean I want to run it on all of them in parallel. I meant I want samtools to sort them out one by one. Usually what I would do would be:

samtools sort -@ 7 File1.bam File1.sorted.bam
samtools sort -@ 7 File2.bam File2.sorted.bam

and then as soon as it finishes with file1 it would go automatically to file2. But having 380 or so files, you can see the issue with having to write all of them out manually :)

ADD REPLY • link written 3.6 years ago by V • 230

Indeed, a for loop better answers to your question. Nevertheless, if you run parallel and limit the number of concurrent jobs to 1 (--jobs 1) you will also be able to run the commands sequentially, basically like in a for loop. And you will have learned a new tool :)

ADD REPLY • link written 3.6 years ago by Charles Plessy • 2.7k

3.6 years ago by

lakhujanivijay ♦ 5.4k

India/Ahmedabad

lakhujanivijay ♦ 5.4k wrote:

try this:

for in in *.bam; do sample_name=`echo $i | awk -F "." '{print $1}'`; samtools sort -@ 7 $i >  ${sample_name}.sorted.bam ; done

ADD COMMENT • link modified 3.6 years ago • written 3.6 years ago by lakhujanivijay ♦ 5.4k

Just a comment, when you are defining output filename, it's always better to use basename in order to remove just the extension and add output format. Something like outFile=$(basename "$BAM" .bam | awk '{print $1".sorted.bam"}').

Why? Because sometimes the filenames include several . which explain a lot about the file. E.g. ICGC_File1.cohortA.rmdup.bam. With your approach you'll get an output which looses some information ICGC_File1.sorted.bam.

ADD REPLY • link modified 3.6 years ago • written 3.6 years ago by venu ♦ 6.7k

+1 Thanks venu for the suggestion, though I am clueless why the loop I suggested does not work for V

ADD REPLY • link written 3.6 years ago by lakhujanivijay ♦ 5.4k

Hello,

Thanks. I've tried it but keep getting this error appearing multiple times:

awk: read error (Is a directory)

Usage:   samtools sort [options] <in.bam> <out.prefix>

Options: -n        sort by read name
         -f        use <out.prefix> as full file name instead of prefix
         -o        final output to stdout
         -l INT    compression level, from 0 to 9 [-1]
         -@ INT    number of sorting and compression threads [1]
         -m INT    max memory per thread; suffix K/M/G recognized [768M]

ADD REPLY • link written 3.6 years ago by V • 230

try this: Removed ">" and its for i in not for in in ; i edited

for i in *.bam; do sample_name=`echo $i | awk -F "." '{print $1}'`; samtools sort -@ 7 $i ${sample_name}.sorted; done

I assume you run it for a directory having files:

sample1.bam
sample2.bam
sample3.bam

and so on

ADD REPLY • link modified 3.6 years ago • written 3.6 years ago by lakhujanivijay ♦ 5.4k

Same Error coming up :/

ADD REPLY • link written 3.6 years ago by V • 230

yes, thats correct. They aren't named "sample" they have a serial number but yes all bam

ADD REPLY • link written 3.6 years ago by V • 230

Please log in to add an answer.

Similar posts • Search »