Question

Convert BAM to FASTQ

1

Entering edit mode

ioannis.vardaxis ▴ 30

@ioannisvardaxis-11763

Last seen 11 months ago

Norway/Oslo

Hey,

In my package I was using rbamtools from CRAN for converting some BAM files to FASTQ. However rbamtools is now depressed, so I was wondering if there is any function in Bioconductor which can do that. I tried to find it myself with no luck.

Best,

BAM FASTQ • 3.2k views

ADD COMMENT • link updated 5.1 years ago by Martin Morgan 25k • written 5.1 years ago by ioannis.vardaxis ▴ 30

1

Entering edit mode

Hi ioannis.vardaxis, I never done this in R environment but I can assure that is very easy using samtools in unix. this is the code:

samtools fastq input.bam > output.fastq

if your BAM is a paired end you have to split the fastq generated

cat output.fastq | grep '^@.*/1$' -A 3 --no-group-separator > sample_mate1.fastq 
cat output.fastq | grep '^@.*/2$' -A 3 --no-group-separator > sample_mate2.fastq

ADD REPLY • link 5.1 years ago Merlin ▴ 10

0

Entering edit mode

Thanks, but I need it in R since it is part of my R-package :/

ADD REPLY • link 5.1 years ago ioannis.vardaxis ▴ 30

score 5 · Accepted Answer · 2019-03-29

Load Rsamtools and ShortRead

library(Rsamtools)
library(ShortRead)

Create a BamFile that references your file, with a yieldSize (number of records available at one time). Open it

bf = BamFile("<path/to/bam">, yieldSize = 1000000)
open(bf)
to = "<path/to/fastq>"

Write a helper function to read and write one chunk

fun <- function(bf, to) {
    chunk <- scanBam(bf, param = ScanBamParam(what = c("seq", "qual")))
    fq <- ShortReadQ(chunk[[1]]$seq, chunk[[1]]$qual)
    writeFastq(fq, to, "a")
    length(fq)
}

and apply this to each chunk

repeat {
    ## maybe indicate progress?
    len <- fun(bf, to)
    if (len == 0L)
        break
}