Calculate coverage from SAM file
I have a really old sam format file that does not have any header or information on what build it has been aligned to. It looks something like this
NS500377:28:H16A7BGXX:3:21410:1181:5368 99 chr11 3000035 255 40M = 3000039 44 GTTTTTCACTGTTTCTCCCCATATTCCAGGTCTTACAGTG AAAAAFA<AFFFFFFFFAFFFFFFFFFFFFFFFFFFFFFF XA:i:1 MD:Z:17G10A11 NM:i:2
NS500377:28:H16A7BGXX:3:21410:1181:5368 147 chr11 3000039 255 40M = 3000035 -44 TTCACTGTTTCTCCCCATATTCCAGGTCTTACAGTGTGTG FFF<FFFFFFFFFFFFFAFFFFAFFFFFFFFFFFFAAAAA XA:i:2 MD:Z:13G10A15 NM:i:2
NS500377:28:H16A7BGXX:2:22105:14414:16397 163 chr11 3000074 255 39M = 3000199 164 GTGTGTTTCTTATTTTGTAGGTTTATTCAGTGTTTGTGG AAAAAFF<FFAFFFFFFFAFFFFF.FFFFFAFFFFFFFA XA:i:1 MD:Z:24T14 NM:i:1
NS500377:28:H16A7BGXX:1:23206:13035:2291 163 chr11 3000187 255 40M = 3000325 178 ACCATAGGTGTGTCTCATTTTCACATTTTATCAGTTTACA AA<AAFF7<FFAFFAFA)7.)FFFAFF)<)F.<F)AF.F. XA:i:0 MD:Z:29T9T0 NM:i:2
I have been asked to check for coverage, ie, calculating how many times part of a read maps to each coordinate in the genome.
But cannot think of any way of doing so other than using the tools such as picard, bedtools etc. But those will not work as this SAM file is incomplete/invalid without header/build info.
Does anyone have a custom script or know of a way to calculate coverage for such files ?
Thank you very much.
• 2,736 views
•
link
updated
by
onestop_data
•
written
by
NB
•
written
by
NB
0 answers
No answers yet.
Log in to answer this question.
Do you know the reference genome? If so you can simply take the header of any BAM file that was produced against the same reference. Extract the @HD and @SQ part of the header and then do something like:
From there on use any standard tool of choice.
unfortunately not. these are really old clinical samples that have been sent across.
See if you can infer the genome version from the first/last coordinates of mapped reads on each chromosome. The chromosome lengths and how much of them are hard masked has changed between genome versions.
If that doesn't work, just choose what you hope to be the correct genome version, since the results won't be off by much.
Samtools recently added a coverage sub-command, didn't it? Try
samtools coverageand make sure you have the latest version of the tool (as I said it was added recently.The problem here is that OP has no information about the genome version, it is not about tools.
I see. In this case, I would re-create the original FASTA from the SAM he has (it contains sequence id and sequence), and would properly realign to the reference (which seems to be the human genome). And then use samtools coverage.
Does it make sense @Nandini?