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?
1 answer
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.
Log in to answer this question.
-
How to set fq.gz as input for bbduk (BBmap)?
written by marongiu.luigi
-
Low SNP Overlap with Michigan 1KG and TopMed reference panel
written by berndmann •
-
Clumpify on multiple files
written by liorglic
-
merge individual runs after bcftools mpileup | bcftools call
written by Leonor •
-
samtools mpileup fail to create bcf
written by mthm •
-
Extracting reads from fastq files with BBMAP?
written by simplitia •
-
Split a FastQ file into smaller files
written by Matthias Zepper
-
extract unmapped pair end read from sam files
written by Bioinfonext
-
Pipe BBMap reformat output into SAMtools
written by bruce.moran
-
How to present "N" in gap region of consensus in sequence assembly?
written by kuosheng0629 •
