Playing with #magicblast, the #NCBI Short read mapper. My notebook

NCBI MAGIC Blast was recently mentioned by BioMickWatson on twitter.

Looks pretty cool. Perhaps once again the answer to all bfx questions will be BLAST
RE https://t.co/4D5e9QQnrb pic.twitter.com/bwW3y0yl2n

- Mick Watson (@BioMickWatson) September 9, 2016

Here, I'll be playing with magicblast and I'll compare its output with bwa (Makefile below).

First, here is an extract of the manual for magicblast.

USAGE
DESCRIPTION
   Short read mapper

 *** Input query options
 -query <File_In>
   Input file name
   Default = `-'
    * Incompatible with:  sra
 -infmt <String, `asn1', `asn1b', `fasta', `fastc', `fastq'>
   Input format for sequences
   Default = `fasta'
    * Incompatible with:  sra
 -paired
   Input query sequences are paired
 -query_mate <File_In>
   FASTA file with mates for query sequences (if given in another file)
    * Requires:  query
 -sra <String>
   Comma-separated SRA accessions
    * Incompatible with:  query, infmt

 *** Formatting options
 -outfmt <String, Permissible values: 'asn' 'sam' 'tabular' >
   alignment view options:
   sam = SAM format,
   tabular = Tabular format,
   text ASN.1
   Default = `sam'

Indexing the reference for magicblast

I've indexed the human chr22 using the good old NCBI makeblastdb:

$ wget -O "chr22.fa.gz" "http://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/chr22.fa.gz"
$ gunzip -f chr22.fa.gz
$ makeblastdb -in chr22.fa -dbtype nucl -title chr22

Mapping Paired-End reads with magicblast

First I've removed the read names' suffixes '/1' and '/2' because they still appear in the final bam.

gunzip -c R1.aa.fq.gz | sed 's/\/1$$//' | gzip > R1.fq.gz
gunzip -c R2.aa.fq.gz | sed 's/\/2$$//' | gzip > R2.fq.gz

The reads are then mapped with magicblast using the following command:

magicblast -db chr22 \
 -infmt fastq -paired \
 -query R1.fq.gz \
 -query_mate R2.fq.gz |\
 sed 's/gnl\|BL_ORD_ID\|0/chr22/g' |\
 samtools sort -o child.magic.bam -T child.magic --reference chr22.fa -O BAM

As far as I could see, there was no option to specify the read group (RG).

Here,I've transformed the name of the contig because it looks like a blast-id identifier:
I've also mapped the reads using bwa:

bwa mem  chr22.fa R1.fq.gz R2.fq.gz |\
 samtools sort -o child.bwa.bam -T child.bwa --reference chr22.fa -O BAM

BAM headers

Here is the BAM header for magicblast: Magic blast uses the spec v1.2 and has a Group-Order flag.

@HD VN:1.2 GO:query SO:coordinate
@SQ SN:chr22 LN:50818468
@PG ID:0 PN:magicblast magicblast -db chr22.fa -infmt fastq -paired -query R1.fq.gz -query_mate R2.fq.gz

And the BAM header for bwa:

@HD VN:1.3 SO:coordinate
@SQ SN:chr22 LN:50818468
@PG ID:bwa PN:bwa VN:0.7.15-r1140 bwa mem chr22.fa R1.fq.gz R2.fq.gz

Samtools samflags

for magicblast:

24387 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
24387 + 0 mapped (100.00% : N/A)
24387 + 0 paired in sequencing
12222 + 0 read1
12222 + 0 read2
24330 + 0 properly paired (99.77% : N/A)
24387 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
57 + 0 with mate mapped to a different chr
57 + 0 with mate mapped to a different chr (mapQ>=5)

for bwa:

24001 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
1 + 0 supplementary
0 + 0 duplicates
23976 + 0 mapped (99.90% : N/A)
24000 + 0 paired in sequencing
12000 + 0 read1
12000 + 0 read2
23840 + 0 properly paired (99.33% : N/A)
23952 + 0 with itself and mate mapped
23 + 0 singletons (0.10% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

Validationg with picard

Validation of child.magic.bam with picard generates many errors:

$ java -jar picard.jar  ValidateSamFile  I=child.magic.bam IGNORE=RECORD_MISSING_READ_GROUP
ERROR: Read groups is empty
ERROR: Header version: 1.2 does not match any of the acceptable versions: 1.0, 1.3, 1.4, 1.5
ERROR: Record 4, Read name HWI-1KL149:97:C6Y6VACXX:4:2306:7588:23627, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 6, Read name HWI-1KL149:97:C6Y6VACXX:5:2303:6772:28953, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 15, Read name HWI-1KL149:97:C6Y6VACXX:4:1207:16508:83772, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 17, Read name HWI-1KL149:97:C6Y6VACXX:4:1216:17305:32405, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 19, Read name HWI-1KL149:97:C6Y6VACXX:5:1315:20109:45547, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 29, Read name HWI-1KL149:97:C6Y6VACXX:4:2313:12478:85202, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 36, Read name HWI-1KL149:97:C6Y6VACXX:4:1202:11437:96159, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 40, Read name HWI-1KL149:97:C6Y6VACXX:4:1213:16611:87818, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 45, Read name HWI-1KL149:97:C6Y6VACXX:4:1208:7944:83271, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 49, Read name HWI-1KL149:97:C6Y6VACXX:4:1216:17305:32405, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 22, Read name HWI-1KL149:97:C6Y6VACXX:5:1105:17083:36022, Mate alignment does not match alignment start of mate
ERROR: Record 51, Read name HWI-1KL149:97:C6Y6VACXX:5:1105:17083:36022, Mate alignment does not match alignment start of mate
ERROR: Record 51, Read name HWI-1KL149:97:C6Y6VACXX:5:1105:17083:36022, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 25, Read name HWI-1KL149:97:C6Y6VACXX:4:2209:8319:82198, Mate negative strand flag does not match read negative strand flag of mate
(...)

while there is ~no error for child.bwa.bam:

$ java -jar picard.jar  ValidateSamFile  I=child.bwa.bam IGNORE=RECORD_MISSING_READ_GROUP
ERROR: Read groups is empty
[Fri Sep 09 16:09:28 CEST 2016] picard.sam.ValidateSamFile done. Elapsed time: 0.01 minutes.
Runtime.totalMemory()=78643200

Variant Calling

Both BAMs are mapped with samtools/vcftools (min DP=10)

samtools mpileup  --uncompressed --output-tags DP --reference chr22.fa child.bwa.bam |\
 bcftools call --ploidy GRCh38  --multiallelic-caller --variants-only --format-fields GQ,GP - |\
 bcftools filter -e 'DP<10'  --output-type v -o child.bwa.vcf -
samtools mpileup  --uncompressed --output-tags DP --reference chr22.fa child.magic.bam |\
 bcftools call --ploidy GRCh38  --multiallelic-caller --variants-only --format-fields GQ,GP - |\
 bcftools filter -e 'DP<10'  --output-type v -o child.magic.vcf -

Number of variants:

$ grep -v "#" child.magic.vcf  | wc -l
111
$ grep -v "#" child.bwa.vcf  | wc -l
166

Here is the diff for CHROM/POS/REF/ALT:

• Uniq to BWA: 'comm -13 <(grep -v "#" child.magic.vcf|cut -f 1,2,4,5 | sort | uniq) <(grep -v "#" child.bwa.vcf|cut -f 1,2,4,5 | sort | uniq) | wc -l' : 64 variants
• Uniq to Magic: 'comm -23 <(grep -v "#" child.magic.vcf|cut -f 1,2,4,5 | sort | uniq) <(grep -v "#" child.bwa.vcf|cut -f 1,2,4,5 | sort | uniq) | wc -l' : 9 variants
• Common variants: 'comm -12 <(grep -v "#" child.magic.vcf|cut -f 1,2,4,5 | sort | uniq) <(grep -v "#" child.bwa.vcf|cut -f 1,2,4,5 | sort | uniq) | wc -l' : 102 variants

The Makefile

SHELL=/bin/bash
data.dir=${HOME}/src/DATA
ncbi.bin=${HOME}/packages/magicblast/bin
REF=chr22.fa
samtools.exe=${HOME}/packages/samtools/samtools 
bwa.exe=${HOME}/packages/bwa-0.7.15/bwa 

all: child.magic.bam child.bwa.bam

R1.fq.gz : ${HOME}/src/DATA/child.R1.aa.fq.gz 
 gunzip -c $< | sed 's/\/1$$//' | gzip > $@
R2.fq.gz : ${HOME}/src/DATA/child.R2.aa.fq.gz 
 gunzip -c $< | sed 's/\/2$$//' | gzip > $@

child.bwa.bam : ${REF}.bwt R1.fq.gz R2.fq.gz 
 ${bwa.exe} mem  ${REF} $(word 2,$^) $(word 3,$^) |\
 ${samtools.exe} sort -o $@ -T $(basename $@) --reference ${REF} -O BAM 


child.magic.bam : ${REF}.nin R1.fq.gz R2.fq.gz 
 ${ncbi.bin}/magicblast -db ${REF} -infmt fastq -paired -query $(word 2,$^) -query_mate $(word 3,$^) |\
 sed 's/gnl\|BL_ORD_ID\|0/$(basename ${REF})/g' |\
 ${samtools.exe} sort -o $@ -T $(basename $@) --reference ${REF} -O BAM 

${REF}.bwt : ${REF}
 ${bwa.exe} index $<

${REF}.nin : ${REF}
 ${ncbi.bin}/makeblastdb -in $< -dbtype nucl -title $(basename ${REF})

${REF} :
 wget -O "$@.gz" "http://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/$@.gz"
 gunzip -f $@.gz

That's it,
Pierre

It's a kind of magic(5)

The following question was recently asked on biostars.org : Tool that detects data types:"More specifically I am interested for detecting between SAM, BAM, FASTA, FASTQ (if possible descriminating between these types), BED, BED12, BED15, GFF, GFF2, GFF3). The detection should be performed without just looking at the extension of the file....":
I suggested to create a new magic file for the linux file command. As an example, a BAM file starts (position=0) with the magic bytes BAM\1. A magic file definition 'bam' for BAM would be:

0 string BAM\1 BAM file v1.0

Compile the new magic file:

file -C -m bam

Now use this magic:

file -z -m bam.mgc file.bam
file.bam: BAM file v1.0 (data)

Now, I've started a github repo containing some 'magic' patterns for bioinformatics (Fastq, blast, bam, bigwig, etc... ):

https://github.com/lindenb/magic

.
(My current problem is to prioritize some results like differentiating a 'Nucleotide' and a 'Protein' Fasta sequence ( http://unix.stackexchange.com/questions/154096/file1-and-magic5-prioritizing-a-result).)

That's it,
Pierre

Parsing a BAM file with javascript, yes we can. (Node.js and V8)

Node.js is an event-driven I/O server-side JavaScript environment based on V8, Google's open source JavaScript engine. In the current post I will describe how I've used Node/V8 to parse a BAM file. Here I've used node v0.5.5 and my code is hosted in a new git repository bionode.

Designing a C++ native extension for Node.js wrapping the bgzf format

BAM files are stored using the bgzf format. We must first create a C++ extension wrapping the methods related to bgzf. This process is nicely described in "Writing Node.js Native Extensions". Here, a BGZF* pointer is wrapped into a class BGZFSupport that extends v8::ObjectWrapp:

class BGZFSupport: public ObjectWrap
 {
 private:
   BGZF* file;
 public:
    (...)
   int close()
    {
    int ret=0;
    if(file!=NULL) ret=::bgzf_close(file);
    file=NULL;
    return ret;
    }
   ~BGZFSupport()
    {
   if(file!=NULL) ::bgzf_close(file);
    }
 (...)

The javascript constructor for BGZFSupport opens the bgzfile and is implemented on the C++ side as:

  static Handle<Value> New(const Arguments& args)
    {
    HandleScope scope;
    if (args.Length() < 2)
      {
      RETURN_THROW("Expected two parameters for bgfz");
      }
    if(!args[0]->IsString())
     {
     RETURN_THROW("1st argument is not a string");
     }
    if(!args[1]->IsString())
     {
     RETURN_THROW("2nd argument is not a string");
     }
    
    v8::String::Utf8Value filename(args[0]);
    v8::String::Utf8Value mode(args[1]);
    BGZF* file= ::bgzf_open(ToCString(filename),ToCString(mode));
    if(file==NULL)
     {
     RETURN_THROW("Cannot open \"" << ToCString(filename) <<  "\"");
     }
    BGZFSupport* instance = new BGZFSupport(file);
    instance->Wrap(args.This());
    return args.This();
    }

... and so on for the other functions...

Implementing the javascript-based BAM-Reader

Next, we can embbed this BGZFSupport in a javascript file that will read a BAM file:

var bgzf=require("bgzf");

and we create a javascript class/function BamReader that will open the file as bgzf and will read the BAM header:

var bgzf=require("bgzf");
var Buffer = require('buffer').Buffer;


function BamReader(path)
 {
 this.fd= new bgzf.bgzf(path,"r");
 var b=new Buffer(4);
 var n = this.fd.read(b,0,4);
 if(n!=4) throw new Error("Cannot read 4 bytes");
 if(b[0]!=66)  throw new Error("Error MAGIC[0]");
 if(b[1]!=65)  throw new Error("Error MAGIC[1] got"+b[1]);
 if(b[2]!="M".charCodeAt(0))  throw new Error("Error MAGIC[2]");
 if(b[3]!="\1".charCodeAt(0))  throw new Error("Error MAGIC[3]");
 
 /* l_text */
 n = this.fd.read(b,0,4);
 if(n!=4) throw new Error("Cannot read 4 bytes");
 var l_text=b.readInt32LE(0);
 b=new Buffer(l_text);
 n = this.fd.read(b,0,l_text);
 if(n!=l_text) throw new Error("Cannot read "+l_text+" bytes (l_text)");
 this.text=b.toString('utf-8', 0, l_text);
 
 /* n_seq */
 b=new Buffer(4);
 n = this.fd.read(b,0,4);
 if(n!=4) throw new Error("Cannot read 4 bytes");
 var n_ref=b.readInt32LE(0);
 this.references=[];
 this.name2seq={};
 for(var i=0;i< n_ref;++i)
  {
  var refseq={};
  /* l_name */
  b=new Buffer(4);
  n = this.fd.read(b,0,4);
  if(n!=4) throw new Error("Cannot read 4 bytes");
  var l_name=b.readInt32LE(0);
  /* name */
  b=new Buffer(l_name);
  n = this.fd.read(b,0,l_name);
  if(n!=l_name) throw new Error("Cannot read "+l_name+" bytes (name)");
  refseq.name=b.toString('utf-8', 0,l_name-1);//\0 terminated
  /* l_ref */
  b=new Buffer(4);
  n = this.fd.read(b,0,4);
  if(n!=4) throw new Error("Cannot read 4 bytes");
  refseq.l_ref=b.readInt32LE(0);
  this.references.push(refseq);
  this.name2seq[refseq.name]=refseq;
  }
 //console.log(this.name2seq);
 }

Another function next() reads the next alignment or returns null ( see the code ).

Testing

$ export NODE_PATH=/path/to/bionode/build

the script reads a simple BAM file and prints the positions of the reads:

(...)

var r= new BamReader("/path/to/samtools-0.1.17/examples/toy.bam");
var align;
while((align=r.next())!=null)
 {
 console.log(
  r.references[align.refID].name+"\t"+
  align.read_name+"\t"+
  align.pos
  );
 }
r.close();

Result

$ node bgzf.js
ref r001 6
ref r002 8
ref r003 8
ref r004 15
ref r003 28
ref r001 36
ref2 x1 0
ref2 x2 1
ref2 x3 5
ref2 x4 9
ref2 x5 11
ref2 x6 13

Remaining questions:

At the moment, I don't know how to correctly package the C++ and javascript files for node.js, how to correctly include the files, how to group the different files under a common 'namespace', etc...

That's It,
Pierre

Memory-Mapping the Human Genome with 'mmap': my notebook

In this post, I've explored how to use a memory-mapped file to handle the 3Go of the Human Genome sequence as if it was entirely loaded in memory.
Via wikipedia: "A memory-mapped file is a segment of virtual memory which has been assigned a direct byte-for-byte correlation with some portion of a file or file-like resource. This resource is typically a file that is physically present on-disk... In computing, mmap is a POSIX-compliant Unix system call that maps files or devices into memory. It is a method of memory-mapped file I/O. It naturally implements demand paging, because initially file contents are not entirely read from disk and do not use physical RAM at all. The actual reads from disk are performed in "lazy" manner, after a specific location is accessed."
Using a C++ program, I'm going to memory-map the fasta sequence of the human genome (indexed with samtools faidx) and search the position of some short-reads using a BoyerMoore algorithm:

Required members:

/* maps a chromosome to its samtools faidx index */
map<string,faidx1_t> name2index;
/* used to get the size of the file */
struct stat buf;
/* genome fasta file file descriptor */
int fd;
/* the mmap (memory mapped) pointer */
char *mapptr;

Opening the mmap

string faidx(fastaFile);
string line;
faidx.append(".fai");
/* open *.fai file */
ifstream in(faidx.c_str(),ios::in);
/* read indexes in the .fai file that was created with samtools */
while(getline(in,line,'\n'))
 {
 faidx1_t index;
 //parse the faidx line...
 //...
 name2index.insert(make_pair(chrom,index));
 }
/* close index file */
in.close();

/* get the whole size of the fasta file */
stat(fasta, &buf);
/* open the fasta file */
fd = open(fastaFile, O_RDONLY);
/* open a memory mapped file associated to this fasta file descriptor */
mapptr = (char*)mmap(0, buf.st_size, PROT_READ, MAP_SHARED, fd, 0);

It's a kind of MAGIC: Getting the base at index 'position-th' of chromosome 'chrom'

std::map<std::string,faidx1_t>::iterator r=name2index.find(chrom);
faidx1_t& index=r->second;
char base=at(&index,position);
(...)
/* returns the base at position 'index' for the chromosome indexed by faidx */
char at(const FaidxPtr faidx,int64_t index)
 {
 long pos= faidx->offset +
  index / faidx->line_blen * faidx->line_len +
  index % faidx->line_blen
  ;
 /* here is the magic: no need to fseek/fread/ftell the file */
 return toupper(mapptr[pos]);
 }

Mapping the short reads

I've hacked a simple Boyer-Moore-Horspool algorithm from ttp://en.wikipedia.org/wiki/Boyer-Moore-Horspool_algorithm. Of course, you wouldn't use this algorithm to map your short reads for real :-) .

Disposing the mmap

/* close memory mapped map */
if(mapptr!=NULL) munmap(mapptr,buf.st_size);
 /* dispose fasta file descriptor */
if(fd!=-1) close(fd);

Compile & run

$  g++ -Wall testmmap.cpp -lz

$ ./a.out -g /path/tp/hg19.fa /path/to/my.fastq.gz

ATCATTTTCCTCCTAACAGATTAAAAATCAAGAAATATAAACCAGATGTAGCAG chr11 93065362
(...)

Source code

That's it,

Pierre