This is a test version of Biostars. For the public version, visit https://www.biostars.org.
samtools sort generates over millions of files

Hi, I am using bwa mem to map illumina pair-end reads to the rat genome, and then using samtools sort to sort the bam file. The bam file is about 15Gb. However, it appeared that over millions of temp files have been generated during samtools sort step and it cannot be completed because it reached the max directory limit of the school computer. I checked the temp files, and it looked like only one read was recorded in each file. That doesn't look right. Does anyone know what could be the problem? Thanks a lot!

alignment genome sequencing
• 5,676 views • link
written by sherry760302 • sherry760302

you should start by adding the command that you used

• 0 views • link
written by Chris Miller Chris Miller

Can you provide the command you used to sort with samtools? samtools sort does produce multiple temp files but millions is not normal.

• 0 views • link
written by GenoMax GenoMax

Here is the command: REF, FASTQ1, FASTQ2, OUTDIR were defined elsewhere:

SAMPLE=NP2

RG="@RG\tID:$SAMPLE\tSM:$SAMPLE\tPL:illumina\tLB:$SAMPLE\tPU:$SAMPLE"

NT=8

MEM="8G"

module load bwa
module load java
module load samtools

bwa mem -t $NT -M -R "$RG" $REF  $FASTQ1 $FASTQ2 | samtools view -bS - > $OUTDIR/$SAMPLE.bam

samtools sort -m $MEM $OUTDIR/$SAMPLE.bam $OUTDIR/$SAMPLE.sorted
• 0 views • link
updated by Alex Reynolds Alex Reynolds • written by sherry760302 • sherry760302

What version of samtools are you using?

• 0 views • link
written by Tonor Tonor

Eeek! that is an ancient version of samtools. Please consider upgrading to the latest. Not sure why samtools put one read per file in temp files.

• 0 views • link
written by GenoMax GenoMax

OK. I will try on a newer version. I just found out we do have samtools/1.3.1, but 0.1.18 is the default setting.... Thanks for the suggestion.

• 0 views • link
written by sherry760302 • sherry760302

Also - what operating system is your school computer?

• 0 views • link
written by Tonor Tonor

I would guess that you forgot to set the -m parameter appropriately. It defines the amount of memory to be used at max before spilling the data to a tmp file. Check if you defined the unit correctly. Simply typing -m 1 would be 1 Byte of memory I think (thus explaining the abundant number tmp files as basically every read gets its own file), but you need something like 1G for appropriate performance.

• 0 views • link
written by ATpoint ATpoint

Actually $MEM is set to 8G in the script. Assigning the memory via a variable must be causing some strange problem here (as linked by @tonor in one of the threads above).

• 0 views • link
written by GenoMax GenoMax

Yeah, it is probably the variable. I had the same issue on our CentOS server. Setting -m in a variable always caused trouble there, but the same script on OS X Mavericks worked fine.

• 0 views • link
written by ATpoint ATpoint

1 answer

Thank you all for your comments. I used the 1.3.1 version of samtools and manually defined -m, and now it works! Here is the command I used:

samtools sort -m 8G -T $OUTDIR/$SAMPLE.sorted -o $OUTDIR/$SAMPLE.sorted.bam $OUTDIR/$SAMPLE.bam
• 0 views • link
updated by Alex Reynolds Alex Reynolds • written by sherry760302 • sherry760302

Log in to answer this question.

More posts like this