Use bbmap reformat.sh to convert from paired fq files to a bam file
1
0
Entering edit mode
10 months ago
berndmann ▴ 10

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?

reformat.sh bbmap • 577 views
ADD COMMENT
3
Entering edit mode
10 months ago
GenoMax 142k

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.

ADD COMMENT

Login before adding your answer.

Traffic: 2042 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6