Convert BAM to a GenomeDataObject

There must be a better way to do this, but I can't figure it out.

For ChIP SEQ analysis I usually start with Eland Export files that I read into R with ShortReads, and convert to GenomeData objects for use with the chipseq package. However, lately I've been stuck with BAM files, which are still mysterious and awkward to me.

In order to get the files into R in a form I can start analyzing, I use Rsamtools, which isn't officially out yet, but it's installed in my environment and can read BAM files. The code below basically reads in the BAM file, extracts out the chromosome, strand, and position, for every read, and builds a GenomeData object.

A genome data object is basically a list of chromosomes, whereby each chromosome has s sublist of "+" and "-" vectors of integers representing alignment positions for each strand.

library(Rsamtools)
library(chipseq)

# create a dataframe for chr name substitutions
chrnames <- data.frame(sc=c("chrI", "chrII", "chrIII", "chrIV", "chrV", "chrVI", "chrVII", "chrVIII", "chrIX", "chrX", "chrXI", "chrXII", "chrXIII", "chrXIV", "chrXV", "chrXVI", "mito"), bam=c("Scchr01", "Scchr02", "Scchr03", "Scchr04", "Scchr05", "Scchr06", "Scchr07", "Scchr08", "Scchr09", "Scchr10", "Scchr11", "Scchr12", "Scchr13", "Scchr14", "Scchr15", "Scchr16", "Scmito"))

# since I run this as a script
# capture command line args
Args <- commandArgs()

bamfile <- Args[4]
sname <- Args[5]


# read in bam data
bd <- scanBam(bamfile)

# get chr names
bamchrs <- levels(bd[[1]][["rname"]])

# create empty list
mygd <- list()
for( i in 1:length(bamchrs) ){
  # create index vectors for positions and strands
  chriv <- bd[[1]][["rname"]] == bamchrs[i]
  chriv[is.na(chriv)] <- FALSE
  plusStriv <- bd[[1]][["strand"]] == "+"
  minusStriv <- bd[[1]][["strand"]] == "-"
  # create genome data-like list
  mygd[[i]] <- list("+"=bd[[1]][["pos"]][chriv & plusStriv], "-"=bd[[1]][["pos"]][chriv & minusStriv])
}

# substitute my names for the ones in the BAM object
names(mygd) <- chrnames[match(bamchrs, chrnames[,"bam"]),"sc"]

# convert to a GenomeData object
mygd <- GenomeData(listData = mygd)

# re-assign the mygd object
# to a name I handed the script
assign(sname, mygd)

# save the GD object to disk
filename <- paste(sname, ".RData", sep="")

save(list=sname, file=filename)

The syntax of the save call might seem a little bit odd (list=sname), but that's what's required to save the object with the new name I assigned it using the "assign" function.

With the script above, I can convert a bunch of BAM files to GD objects using a shell script as follows:

#!/usr/local/bin/bash

cat bam2GD.r | R --vanilla --args lane1.bam mutant_wce
cat bam2GD.r | R --vanilla --args lane2.bam mutant_ip
cat bam2GD.r | R --vanilla --args lane3.bam wt_wce
cat bam2GD.r | R --vanilla --args lane4.bam wt_ip

I save the above code ina file called: bam2gd.sh, and it creates a bunch of GenomeData objects for further analysis.