Hello,

I have 50 bam files that I'd like to go through.  I have a list of mutations in 21 genes in those bamfiles, and I'd like to write some code to go through each gene, locate the base that is mutated, and get a readout of what's there in the tumor sample and the normal sample.

In the past I've done this manually in IGV, but I've gotten better at coding and would like to give this a try.  I've used several bioconductor packages before, but never these.  I've used samtools sparingly on the command line, so am somewhat familiar with it, but am no expert.

What I need help with is using rsamtools or rbamtools to extract a base from the genome.  For example, I would give it chr10: 387,878,500-387,878,500, and it would give me some information.  Information I could use is: 1) all the bases read at that position, 2) information on how many of each base are read at that position (IGV shows the total reads and # of reads with each base pair when you hover.  This would work).

Is this possible?  Thanks for your help!

Hi,

Some functions in Rsamtools and GenomicAlignments package can fullfil this task.

Following were from the help document of GenomicAlignments package

1 stackStringsFromBam function for stacking the read sequences (or their quality strings) stored in a BAM file on a region of interest.

2 pileLettersAt extracts the letters/nucleotides of a set of reads that align to a set of genomic positions of interest. The extracted letters are returned as "piles of letters" (one per genomic position of interest) stored in an XStringSet (typically DNAStringSet) object.

3 The pileup and applyPileups functions defined in the Rsamtools package, as well as the SAMtools mpileup command (available at http://samtools.sourceforge.net/ as part of the SAMtools project), for more powerful flexible alternatives.

---

Gang