Select concordantly mapped reads from a sam file without header

Question

[valerie]

Hi everyone.

I have paired-end sam file without header. Unfortunately I removed it earlier as I didn't suppose I will need it. Now I am trying to select only concordantly mapped reads using the following command:

samtools view -f 0x2 myfile.sam > myfile_concordant.sam

But I get an error:

[main_samview] truncated file.[E::sam_parse1] missing SAM header
[W::sam_read1] parse error at line 1
[main_samview] truncated file.

Is there any possibility to select only concordant reads? Maybe somehow filter sam file without samtools view?

Thanks!

Answer 1

[valerie]

I guess I found the solution here: http://www.acgt.me/blog/2015/4/15/filtering-a-sam-file-generated-by-tophat-to-find-uniquely-mapped-concordant-read-pairs-awk-vs-samtools

So if you have a TopHat SAM file, and you wanted to filter it to only keep uniqueley mapped, concordant read pairs, you could use two of the options that the samtools view command provides:

-q 50 — This filters on the MAPQ field (5th column of SAM file). TopHat uses a value of 50 in this field to denote unique mappings (this important piece of information is not in the TopHat manual). -f 0x2 — This option filters on the bitwise FLAG score (column 2 of SAM file), and will extract only those mappings where the 2nd bit is set. From the SAM documentation, this is described as each segment properly aligned according to the aligner. In practice this means looking for values that are either 83, 89, 147, or 163 (see this helpful Biobeat blog post for more information about this).