This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Split SAM file into smaller SAM files by line count

So, I am trying to split a SAM file into many smaller SAM files -- by line number, where the line number is a multiple of 8.

I tried the following series of commands:

samtools view -H Mapped_Pila/LBG07/LBG07_pe.sam > Mapped_Pila/LBG07/LBG07_pe.header.sam && samtools view Mapped_Pila/LBG07/LBG07_pe.sam | split -l 8000000 - Mapped_Pila_Split/LBG07/LBG07_pe.split. && find -L . | grep Mapped_Pila_Split/LBG07/LBG07_pe.split. | parallel --gnu -j4 "echo Mapped_Pila/LBG07/LBG07_pe.header.sam {} >> {}.tmp && rm -f {}" && rename "s/\.tmp$/\.sam/" Mapped_Pila_Split/LBG07/LBG07_pe.split.*

Unpacking this:

1) get the header in a separate file

2) get the headerless SAM alignments and split them by a multiple of 8 (8000000)

3) find the split files from the original files, and for each of them, cat the header, then the headerless records into a new file

4) Remove the original split file

5) rename the tmp files to sam

This results in a SAM file that is parsed as truncated by samtools. So, I think I am missing something.

Generally, I want to split sam files that are very large into more tractable parts and then later recombine them with samtools merge.

sam
• 7,152 views • link
updated by simon.burgermeister • simon.burgermeister • written by cacampbell • cacampbell

3 answers

You should add the header back to each of the splitted files.

something line cat Mapped_Pila/LBG07/LBG07_pe.header.sam split_file > split_file.sam

It should be 

cat Mapped_Pila/LBG07/LBG07_pe.header.sam {} > {}.tmp && rm -f {}

in your command.

• 1 views • link
written by GouthamAtla GouthamAtla

YEP... I definitely did use echo instead of cat. Oops.

• 1 views • link
written by cacampbell • cacampbell

I actually end up doing this with some regularity, so here's a smooth all-in-one method:

samtools view -H yourfile.bam > header

then

samtools view yourfile.bam | split - yourprefix -l 8000000 --filter='cat header - > $FILE.sam'

or

samtools view yourfile.bam | split - yourprefix -l 8000000 --filter='cat header - | samtools view -b - > $FILE.bam'
• 1 views • link
written by george.ry george.ry

Thanks, this does the trick -- I think my problem was that I stupidly used "echo" instead of "cat" -- so of course the file was invalid.

• 1 views • link
written by cacampbell • cacampbell

With python using pysam:

import pysam base_path_save=path_2_save samfile = pysam.AlignmentFile("sequences5/aligned.sam", "r") header1=samfile.header.to_dict() for read in samfile.fetch(): file_path=base_path_save+read.query_name+'.sam' with pysam.AlignmentFile(file_path, "w", header=header1) as outf: outf.write(read)

• 0 views • link
written by simon.burgermeister • simon.burgermeister

Log in to answer this question.

More posts like this