Issue with Sam-Bam conversion samtools - how to remove last line of Sam file?

[TabeaK]

Hello guys and gals!

I have a bit of a problem with samtools Sam-bam conversion. I have 8 files containing single read illumina data that I have aligned to my reference genome using bwa samse. No problems. Upon conversion, however, 4 out of 8 Sam files create a problem. The other 4 convert without a hitch.

I am writing this from my mobile so cannot copy in the error message, but it is the same that has been discussed here before - samtools complains about the sequence and quality information not matching and aborts.

Funnily, this only happens in the very last line of each Sam file! So I figured I just remove it and try again.

Any ideas as to how to accomplish this? I have tried opening the file in a simple text editor, but they are too massive (15GB)...

[dpryan]

If all you want to do is trim off the last line in a file, I can think of a whole slew of ways. You could count the number of lines in the files (wc -l) and then "head -n X" into a new SAM files, where X is 1 less than the number of lines in the file. Alternatively, you could just grep for everything put the last read name and send that to a new file. Or, you could just tell samtools how many of the reads to convert to a BAM file (the -c switch), giving it all but the last one.

Of course, this assumes that removing the last line will actually fix things. It'd be helpful to see the last read that seems to be causing the problem and the error message.

[xied75]

Agree that one line might only be a tip of iceberg.

[TabeaK]

Thanks for your input!

I have, however, since sorted and merged all my bam files. Samtools complains about the missing EOF in the 4 bams in question and warns of possibly truncated data, but creates a merged file, which contains all my reads ( I can see that from samtools flagstat) but the 4 last ones from the respective files. The alignment looks absolutely fine when visualized in IGV.

Anyway, it is of course better to find out where this pesky error is coming from, I'll figure out how to look at those lines tomorrow...read about the bamseek tool, maybe I'll try that.