CompareBams4

Last commit

Compare two query-name sorted BAM files. Print a tab-delimited report

Usage

This program is now part of the main jvarkit tool. See jvarkit for compiling.

Usage: java -jar dist/jvarkit.jar cmpbams4  [options] Files

Usage: cmpbams4 [options] Files
  Options:
    -c, --chain
      LiftOver chain file. Can be a local chain file, a URL 'https://hgdownload.soe.ucsc.edu/goldenpath/hg19/liftOver/hg19ToCriGri1.over.chain.gz', 
      or a chain identifier like 'hg19ToHg38'.. Lift Over file from bam1 to 
      bam2. Optional
    -h, --help
      print help and exit
    --helpFormat
      What kind of help. One of [usage,markdown,xml].
    -m, --mismatch
      Default Lift Over mismatch. negative=use default
      Default: 0.95
    -novalidchain, --novalidchain
      Disable Lift Over chain validation
      Default: false
    -o, --output
      Output file. Optional . Default: stdout
    -R, --reference
      For Reading CRAM. Reference for first BAMIndexed fasta Reference file. 
      This file must be indexed with samtools faidx and with picard/gatk 
      CreateSequenceDictionary or samtools dict
    -R2, --reference2
      For Reading CRAM. Reference for second BAM, if different from 1st bam. 
      Indexed fasta Reference file. This file must be indexed with samtools 
      faidx and with picard/gatk CreateSequenceDictionary or samtools dict
    -sortmethod, --sortmethod
      [20171110]Method used to sort the read on query name. (samtools != 
      picard) see https://github.com/samtools/hts-specs/issues/5
      Default: picard
      Possible Values: [samtools, picard]
    --version
      print version and exit

Keywords

  • sam
  • bam
  • compare

Creation Date

20161206

Source code

https://github.com/lindenb/jvarkit/tree/master/src/main/java/com/github/lindenb/jvarkit/tools/cmpbams/CompareBams4.java

Contribute

License

The project is licensed under the MIT license.

Citing

Should you cite cmpbams4 ? https://github.com/mr-c/shouldacite/blob/master/should-I-cite-this-software.md

The current reference is:

http://dx.doi.org/10.6084/m9.figshare.1425030

Lindenbaum, Pierre (2015): JVarkit: java-based utilities for Bioinformatics. figshare. http://dx.doi.org/10.6084/m9.figshare.1425030

Example

The following Makefile compare the bam for hg19 and hg38 on chr22 and 21


include ../../config/config.mk

CHROMS=21 22
OUTDIR=tmp

define run
${OUTDIR}/$(1).bam : ${OUTDIR}/$(1).fa.bwt R1.fastq.gz R1.fastq.gz
    ${bwa.exe} mem -M   -R '@RG\tID:SAMPLE\tLB:SAMPLE\tSM:SAMPLE\tPL:illumina\tCN:Nantes' ${OUTDIR}/$(1).fa $$(word 2,$$^) $$(word 3,$$^) |  ${samtools.exe} view -b -u -S -F4 - | ${samtools.exe} sort -n -o $$@ -T ${OUTDIR}/$(1)_tmp -

${OUTDIR}/$(1).fa.bwt : ${OUTDIR}/$(1).fa
    ${bwa.exe} index $$<

${OUTDIR}/$(1).dict : ${OUTDIR}/$(1).fa
    ${java.exe} -jar $(picard.jar) CreateSequenceDictionary R=$$< O=$$@

${OUTDIR}/$(1).fa.fai : ${OUTDIR}/$(1).fa
    ${samtools.exe} faidx $$<

${OUTDIR}/$(1).fa : 
    mkdir -p $$(dir $$@) && rm -f $$@
    $$(foreach C,${CHROMS}, curl "http://hgdownload.cse.ucsc.edu/goldenPath/$(1)/chromosomes/chr$${C}.fa.gz" | gunzip -c >> $$@;)

endef

all:  ${OUTDIR}/diff.txt 

${OUTDIR}/diff.txt : ${OUTDIR}/hg19.bam ${OUTDIR}/hg38.bam ${OUTDIR}/hg19ToHg38.over.chain ${OUTDIR}/hg19.dict ${OUTDIR}/hg38.dict
    mkdir -p $(dir $@) &&  $(call run_jvarkit,cmpbams4) --novalidchain -st -c $(word 3,$^) $(word 1,$^)  $(word 2,$^) > $@

${OUTDIR}/hg19ToHg38.over.chain :
    mkdir -p $(dir $@) && curl "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.gz" | gunzip -c > $@

$(eval $(call run,hg19))
$(eval $(call run,hg38))

Output

onlyIn  liftover    compareContig   shift   diffCigarOperations diffNM  diffFlags   diffChroms  Count
BOTH    SameChrom   DiscordantContig    .   -1  0   147/163 chr22/chr21 2
BOTH    SameChrom   SameContig  Gt100   3   15  83/83   chr22/chr22 1
BOTH    SameChrom   DiscordantContig    .   3   5   147/129 chr21/chr22 1
BOTH    SameChrom   SameContig  Gt100   -1  1   163/163 chr22/chr22 22
BOTH    SameChrom   DiscordantContig    .   0   1   83/99   chr21/chr22 32
BOTH    SameChrom   SameContig  Gt100   0   2   99/99   chr21/chr21 22
BOTH    SameChrom   DiscordantContig    .   2   6   81/65   chr22/chr21 1
BOTH    SameChrom   SameContig  Gt100   0   0   185/137 chr22/chr22 20
BOTH    SameChrom   SameContig  Zero    0   0   177/177 chr22/chr22 1417
(...)