YOKOFAKUN: April 2013

26 April 2013

Java JNI bindings for BWA(mem-lite)

Motivation

BWA 7.4(http://bio-bwa.sourceforge.net/) contains a small C example(https://github.com/lh3/bwa/blob/master/example.c) for running bwa-mem as a library (bwamem-lite). I created some JNI bindings to see if I can bind the C bwa library to java and get the same output than bwamem-lite. I put the code on github at https://github.com/lindenb/jbwa.

Example

(compare to https://github.com/lh3/bwa/blob/master/example.c )

System.loadLibrary("bwajni");
BwaIndex index=new BwaIndex(new File("hg19.fa"));
BwaMem mem=new BwaMem(index);
KSeq kseq=new KSeq(new File("input.fastq.gz");
ShortRead read=null;
while((read=kseq.next())!=null)
        {
        for(AlnRgn a: mem.align(read))
                {
                if(a.getSecondary()>=0) continue;
                System.out.println(  read.getName()+"\t"+  a.getStrand()+"\t"+  a.getChrom()+"\t"+
                        a.getPos()+"\t"+ a.getMQual()+"\t"+ a.getCigar()+"\t"+  a.getNm() );
                }
        }
kseq.dispose();
index.close();
mem.dispose();

Testing

Here is the ouput of the JAVA version:

gunzip -c input.fastq.gz | head -n 4000 |\
java  -Djava.library.path=src/main/native -cp src/main/java \
   com.github.lindenb.jbwa.jni.Example human_g1k_v37.fasta -| tail 


HWI-1KL149:20:C1CU7ACXX:4:1101:3077:33410       +       3       38647538        60      89M11S  1
HWI-1KL149:20:C1CU7ACXX:4:1101:3396:33445       +       8       52567289        60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:10013:33288      -       1       156104115       60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:10390:33496      -       6       123824853       60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:13537:33483      +       2       157367092       60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:14139:33390      +       20      31413797        60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:14514:33458      +       2       179401813       60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:15292:33282      +       15      63335820        60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:16960:33276      -       12      110782784       60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:17355:33322      +       6       126077895       60      100M    1

And the ouput of the Native C version:

gunzip -c input.fastq.gz | head -n 4000 |\
bwa-0.7.4/bwamem-lite human_g1k_v37.fasta - | tail 

HWI-1KL149:20:C1CU7ACXX:4:1101:3077:33410       +       3       38647538        60      89M11S  1
HWI-1KL149:20:C1CU7ACXX:4:1101:3396:33445       +       8       52567289        60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:10013:33288      -       1       156104115       60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:10390:33496      -       6       123824853       60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:13537:33483      +       2       157367092       60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:14139:33390      +       20      31413797        60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:14514:33458      +       2       179401813       60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:15292:33282      +       15      63335820        60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:16960:33276      -       12      110782784       60      100M    1
HWI-1KL149:20:C1CU7ACXX:4:1101:17355:33322      +       6       126077895       60      100M    1

GUI

As a test I also created a swing-Based interface for BWA:

java  -Djava.library.path=src/main/native  -cp src/main/java \
    com.github.lindenb.jbwa.jni.BwaFrame human_g1k_v37.fasta

That's it,

Pierre

23 April 2013

playing with BWA-MEM : my notebook

BWA-MEM, is a new tool that is part of the latest version of BWA. As said Heng Li on Biostar bwa-mem can be used to identify the viral integration sites in the human genome. Here I've used various short reads to explore how bwa maps the pairs. The sequences I used below are:

NOTCH2a and NOTCH2b two sequences on the chromosomes 1, on the same strand
NOTCH2del : same as NOTCH2a but with a deletion
ROXAN1a and ROXAN1b: two sequences on the chromosome 22, on the same strand
NOTCH2ins is NOTCH2a with and insertion of ROXAN1a

The following results were generated with the makefile below. Each pair of fastq was submitted twice (1_* and 2_*). I also added some pairs of reads from another source of FASTQs to let bwa infer the size of the fragments (see biostar: http://www.biostars.org/p/69694).

NOTCH2a vs NOTCH2b BWA-MEM options: none

NOTCH2a vs NOTCH2b
BWA-MEM options: none
NOT Properly mapped reads because they're on the same strand.
Read	Flag	Chrom	Pos	qual	Cigar	Extra
1_Notch2	p1=65	1	120572000	9	50M	NM:i:0 AS:i:50 XS:i:45
1_Notch2	p2=129	1	120572200	9	50M	NM:i:0 AS:i:50 XS:i:45
2_Notch2	p1	1	120572000	9	50M	NM:i:0 AS:i:50 XS:i:45
2_Notch2	p2	1	120572200	9	50M	NM:i:0 AS:i:50 XS:i:45

NOTCH2a vs revcomp(NOTCH2b) BWA-MEM options: none

NOTCH2a vs revcomp(NOTCH2b)
BWA-MEM options: none
Properly mapped reads.
Read	Flag	Chrom	Pos	qual	Cigar	Extra
1_Notch2	pPR1=99	1	120572000	39	50M	NM:i:0 AS:i:50 XS:i:45
1_Notch2	pPr2=147	1	120572200	39	50M	NM:i:0 AS:i:50 XS:i:45
2_Notch2	pPR1	1	120572000	39	50M	NM:i:0 AS:i:50 XS:i:45
2_Notch2	pPr2	1	120572200	39	50M	NM:i:0 AS:i:50 XS:i:45

NOTCH2del vs revcomp(NOTCH2b) BWA-MEM options: none

NOTCH2del vs revcomp(NOTCH2b)
BWA-MEM options: none
Deletion in Read1
Read	Flag	Chrom	Pos	qual	Cigar	Extra
1_Notch2del	pPR1=99	1	120572000	39	37M20D31M	NM:i:20 AS:i:42 XS:i:37
1_Notch2del	pPr2=147	1	120572200	39	50M	NM:i:0 AS:i:50 XS:i:45
2_Notch2del	pPR1	1	120572000	39	37M20D31M	NM:i:20 AS:i:42 XS:i:37
2_Notch2del	pPr2	1	120572200	39	50M	NM:i:0 AS:i:50 XS:i:45

NOTCH2ins vs revcomp(NOTCH2b) BWA-MEM options:

NOTCH2ins vs revcomp(NOTCH2b)
BWA-MEM options: non
Soft clipping of chr1: two hits on different chromosomes.
Read	Flag	Chrom	Pos	qual	Cigar	Extra
1_NOTCH2insRoxan	pR1=97	22	41721566	60	23S50M27S	NM:i:0 AS:i:50 XS:i:0
1_NOTCH2insRoxan	pr2=145	1	120572200	9	50M	NM:i:0 AS:i:50 XS:i:45
2_NOTCH2insRoxan	pR1	22	41721566	60	23S50M27S	NM:i:0 AS:i:50 XS:i:0
2_NOTCH2insRoxan	pr2	1	120572200	9	50M	NM:i:0 AS:i:50 XS:i:45

NOTCH2a+ROXAN1a vs revcomp(NOTCH2b) BWA-MEM options: none

NOTCH2a+ROXAN1a vs revcomp(NOTCH2b)
BWA-MEM options: none
two hits on chromosome chr1 and chr22 for the 5' seq
nevertheless, properly paired was set.
Read	Flag	Chrom	Pos	qual	Cigar	Extra
1_Notch2Roxan	pPR1=99	1	120572000	9	50M50S	NM:i:0 AS:i:50 XS:i:45 XP:Z:22,+41721566,50S50M,9,0;
1_Notch2Roxan	pPR1=99	22	41721566	9	50S50M	NM:i:0 AS:i:50 XS:i:0 XP:Z:1,+120572000,50M50S,9,0;
1_Notch2Roxan	pPr2=147	1	120572200	9	50M	NM:i:0 AS:i:50 XS:i:45
2_Notch2Roxan	pPR1	1	120572000	9	50M50S	NM:i:0 AS:i:50 XS:i:45 XP:Z:22,+41721566,50S50M,9,0;
2_Notch2Roxan	pPR1	22	41721566	9	50S50M	NM:i:0 AS:i:50 XS:i:0 XP:Z:1,+120572000,50M50S,9,0;
2_Notch2Roxan	pPr2	1	120572200	9	50M	NM:i:0 AS:i:50 XS:i:45

NOTCH2a+ROXAN1a vs revcomp(NOTCH2b+ROXAN1b) BWA-MEM options: none

NOTCH2a+ROXAN1a vs revcomp(NOTCH2b+ROXAN1b)
BWA-MEM options: none
All fragments on chr1/chr22.
Reads are NOT properly paired.
Read	Flag	Chrom	Pos	qual	Cigar	Extra
1_Notch2Roxan	pR1=97	1	120572000	9	50M50S	NM:i:0 AS:i:50 XS:i:45 XP:Z:22,+41721566,50S50M,9,0;
1_Notch2Roxan	pR1=97	22	41721566	9	50S50M	NM:i:0 AS:i:50 XS:i:0 XP:Z:1,+120572000,50M50S,9,0;
1_Notch2Roxan	pr2=145	22	41721766	60	50S50M	NM:i:0 AS:i:50 XS:i:0 XP:Z:1,-120572200,50M50S,9,0;
1_Notch2Roxan	pr2=145	1	120572200	9	50M50S	NM:i:0 AS:i:50 XS:i:45 XP:Z:22,-41721766,50S50M,60,0;
2_Notch2Roxan	pR1	1	120572000	9	50M50S	NM:i:0 AS:i:50 XS:i:45 XP:Z:22,+41721566,50S50M,9,0;
2_Notch2Roxan	pR1	22	41721566	9	50S50M	NM:i:0 AS:i:50 XS:i:0 XP:Z:1,+120572000,50M50S,9,0;
2_Notch2Roxan	pr2	22	41721766	60	50S50M	NM:i:0 AS:i:50 XS:i:0 XP:Z:1,-120572200,50M50S,9,0;
2_Notch2Roxan	pr2	1	120572200	9	50M50S	NM:i:0 AS:i:50 XS:i:45 XP:Z:22,-41721766,50S50M,60,0;

NOTCH2a+ROXAN1a vs revcomp(NOTCH2b) BWA-MEM options: -C

NOTCH2a+ROXAN1a vs revcomp(NOTCH2b)
BWA-MEM options: -C
I cannot see the effect of the option '-C':
"append FASTA/FASTQ comment to SAM output"
Read	Flag	Chrom	Pos	qual	Cigar	Extra
1_Notch2Roxan	pPR1=99	1	120572000	9	50M50S	NM:i:0 AS:i:50 XS:i:45 XP:Z:22,+41721566,50S50M,9,0;
1_Notch2Roxan	pPR1=99	22	41721566	9	50S50M	NM:i:0 AS:i:50 XS:i:0 XP:Z:1,+120572000,50M50S,9,0;
1_Notch2Roxan	pPr2=147	1	120572200	9	50M	NM:i:0 AS:i:50 XS:i:45
2_Notch2Roxan	pPR1	1	120572000	9	50M50S	NM:i:0 AS:i:50 XS:i:45 XP:Z:22,+41721566,50S50M,9,0;
2_Notch2Roxan	pPR1	22	41721566	9	50S50M	NM:i:0 AS:i:50 XS:i:0 XP:Z:1,+120572000,50M50S,9,0;
2_Notch2Roxan	pPr2	1	120572200	9	50M	NM:i:0 AS:i:50 XS:i:45

NOTCH2a+ROXAN1a vs revcomp(NOTCH2b) BWA-MEM options: -M

NOTCH2a+ROXAN1a vs revcomp(NOTCH2b)
BWA-MEM options: -M
Set the duplicate flags for the multiple hits.
Read	Flag	Chrom	Pos	qual	Cigar	Extra
1_Notch2Roxan	pPR1=99	1	120572000	9	50M50S	NM:i:0 AS:i:50 XS:i:45 XP:Z:22,+41721566,50S50M,9,0;
1_Notch2Roxan	pPR1s=355	22	41721566	9	50S50M	NM:i:0 AS:i:50 XS:i:0 XP:Z:1,+120572000,50M50S,9,0;
1_Notch2Roxan	pPr2	1	120572200	9	50M	NM:i:0 AS:i:50 XS:i:45
2_Notch2Roxan	pPR1	1	120572000	9	50M50S	NM:i:0 AS:i:50 XS:i:45 XP:Z:22,+41721566,50S50M,9,0;
2_Notch2Roxan	pPR1s	22	41721566	9	50S50M	NM:i:0 AS:i:50 XS:i:0 XP:Z:1,+120572000,50M50S,9,0;
2_Notch2Roxan	pPr2	1	120572200	9	50M	NM:i:0 AS:i:50 XS:i:45

NOTCH2ins vs revcomp(ROXAN1b) BWA-MEM options: none

NOTCH2ins vs revcomp(ROXAN1b)
BWA-MEM options: none
Reads are mapped on chr22.
Read	Flag	Chrom	Pos	qual	Cigar	Extra
1_NOTCH2insRoxan	pPR1	22	41721566	60	23S50M27S	NM:i:0 AS:i:50 XS:i:0
1_NOTCH2insRoxan	pPr2	22	41721766	60	50M	NM:i:0 AS:i:50 XS:i:0
2_NOTCH2insRoxan	pPR1	22	41721566	60	23S50M27S	NM:i:0 AS:i:50 XS:i:0
2_NOTCH2insRoxan	pPr2	22	41721766	60	50M	NM:i:0 AS:i:50 XS:i:0

That's it,

Pierre

YOKOFAKUN