Compare mapped reads in SAM/BAM files

[Brian Bushnell]

It's possible to do so with BBTools like this:

reformat.sh in=file1.sam out=mapped1.sam mappedonly
reformat.sh in=file2.sam out=mapped2.sam mappedonly

That gets you the mapped reads only. Then:

filterbyname.sh in=mapped1.sam names=mapped2.sam out=shared.sam include=t

...which gets you the set intersection;

filterbyname.sh in=mapped1.sam names=mapped2.sam out=only1.sam include=f
filterbyname.sh in=mapped2.sam names=mapped1.sam out=only2.sam include=f

...which get you the set subtractions.

[cmbetts]

If you can get the read names into two R character vectors, you could use the various set operations like setdiff.
There's probably a much more elegant (and faster) way to do this, but quick and dirty:

library("GenomicAlignments");
reads_alnmt1 <- names(readGAlignments(bam_file1, use.names=TRUE));
reads_alnmt2 <- names(readGAlignments(bam_file2, use.names=TRUE));

unique1 <- setdiff(reads_alnmt1, reads_alnmt2);
unique2 <- setdiff(reads_alnmt2, reads_alnmt1);

[Marisa_Miller]

Hello,

I have two different SAM/BAM files for a sequenced library (Illumina library sequenced paired-end) which were generated by mapping the reads (with Bowtie2) to two different reference genomes. I would like to see how many mapped reads (can be mapped concordant or discordant, doesn't matter) are shared between the two alignment files and how many mapped reads are unique to one file or the other.

Is there a tool available to do this, possibly in samtools? If not, does anyone have some suggestions to get me started on the right track?

Thank you!