Use bbmap reformat.sh to convert from paired fq files to a bam file

Question

[berndmann]

As outlined here I was able to create paired-end fastq files with the help of GenoMax . Now I wonder how I can use reformat.sh from bbmap to convert this files to a valid bam file given some reference.fasta.

If I try to use

call = glue("bbmap/reformat.sh in=data/1kg_hg37/out/embryo1_reads_R1.fq.gz in2=data/1kg_hg37/out/embryo1_reads_R2.fq.gz out=embryo1_reads.bam samplerate=0.004 ref=hg19.fa") 
system(call)

This is generating a bam file but if I use samtools mpileup to extract the reads for each position the output is just empty. What am I doing wrong here?

Answer 1

[GenoMax]

One can have a BAM file that contains unaligned sequence data. This is what you are getting from reformat.sh step outlined above. I am not sure what you mean by "valid" BAM file.

If you need to use the BAM file for downstream analysis with other tools they likely expect the BAM file to contain actual alignments against a reference. You will need to use bbmap.sh (the aligner from BBTools or any other suitable other aligner) to create that aligned BAM file. You can't use reformat.sh to create a real aligned data containing BAM file.