How do I compare reads between two bam file?

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:

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

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

thanks

samtools reads bam compare

• 8,366 views • link

updated 4.8 years ago by Pierre Lindenbaum

• written 4.8 years ago by Assa Yeroslaviz

1 answer

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

• 0 views • link

written 4.8 years ago by Pierre Lindenbaum

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?

• 0 views • link

written 4.8 years ago by Assa Yeroslaviz

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)

• 0 views • link

written 4.8 years ago by Pierre Lindenbaum

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.

• 0 views • link

written 4.8 years ago by Assa Yeroslaviz

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

• 0 views • link

written 4.8 years ago by Pierre Lindenbaum