How do I compare reads between two bam file?

Question

[Assa Yeroslaviz]

I would like to find reads in one bam files which are not in a second bam file.

Is there an easy way to do this?

My way for now is:

1. sort both bam files by read name
2. extract read headers of both files (using cut -f1) for easier handling
3. run diff on both files
4. manually search for differences

But this seems to me very tedious. Is there a better way?

thanks

Answer 1

[Pierre Lindenbaum]

I wrote CompareBams https://lindenb.github.io/jvarkit/CompareBams.html (I haven't used this for years...)

[Assa Yeroslaviz]

Thanks, my goal is to identify reads which are in bam1 but are not in bam2. I would go about it by extracting the reads in bam1 using samtools view, then list the reads of both bam files using samgrep. Next, i can comapre them. But how do I see if a read is in one and not in the other? Can it show me that?

[Pierre Lindenbaum]

Thanks, my goal is to identify reads which are in bam1 but are not in bam2.

how, then you just need comm

 comm -3 \
      <(samtools view input1.bam | cut -f 1 | sort | uniq) \
      <(samtools view input2.bam | cut -f 1 | sort | uniq)

[Assa Yeroslaviz]

Thanks again, i have found your tool, but I get this strange error:

ava -jar jvarkit/dist/commbams.jar --samtools \
>  ../bamUnmapped/L75442_Track-118107.Unmapped.sorted.bam ../bamFiles/L75442_Track-118107.Aligned.sortedByCoord.out.bam 
[INFO][CommBams]Expected two and only two bams please, but got 3
[INFO][Launcher]commbams Exited with failure (-1)

any ideas what causes it? I do have only two bam files in the input.

[Pierre Lindenbaum]

this option --samtools was removed , the doc hasn't been updated.