Convert a BAM file from one reference to another?

Question

I have a set of BAM files that are aligned using the NCBI GRCh37 human genome reference (with the chromosome names as NC_000001.10) but I want to analyze it using a BED file that has the UCSC hg19 chromosome names (e.g. chr1). I want to use bedtools to pull out all the on-target and off-target reads.

1. Are NCBI and UCSC directly comparable? Or do I need to re-align the BAM/lift-over the BED to the UCSC reference?
2. Should I convert the BED file or the BAM file? Everyone here uses the UCSC chromosome names/positions so I'll need to convert the eventual files to UCSC anyway.

Answer 1

You're the second person I have ever seen using NCBI "chromosome names" (they're more like supercontig IDs). Normally I would point you to a resource providing mappings between chromosome names, but since no one has added NCBI names (yet, maybe I'll add them now) you're currently out of luck there.

Anyway, the quickest way to do what you want is to samtools view -H foo.bam > header to get the BAM header and then change each NCBI "chromosome name" to its corresponding UCSC chromosome name. DO NOT REORDER THE LINES! You can then use samtools reheader and be done.

Why, you might ask, would this work? The answer is that chromosome/contig names in BAM files aren't stored in each alignment. Rather, the names are stored in a list in the header and each alignment just contains the integer index into that list (read group IDs are similar, for what it's worth). This also leads to the warning above against reordering entries, since that's a VERY convenient way to start swapping alignments between chromosomes.

As an aside, you'd be well served switching to Gencode or Ensembl chromosome names, they're rather more coherent than the something_random mess that's present in hg19 from UCSC.

Update: Because I'm nice, here is the conversion between NCBI and UCSC. Note that if you have any alignments to patches that there is simply no UCSC equivalent. One of the many reasons not to use UCSC (avoid their annotations too).

Answer 2

The "right" solution would be realignment, but that's expensive and most of us would not go that route. My preferred solution would be to convert the bed file, as opposed to the bam. Here's why:

1) Reheadering the bam means that you may have reads aligned to contigs without a corresponding entry in UCSC (see Devon's list for the mappings). This is a problem because:

• Some of those reads would likely have been mapped elsewhere if a reference without those contigs was used.
• I'm not even sure what happens to those reads after reheadering - I guess they would need to be marked as unmapped? Lots of potential for screwiness there.

2) It seems cleaner to convert the bed file from UCSC->NCBI, where you are guaranteed that every entry has a "home". Then, after you pull your info from the bam, you can always convert chromosome names back if you need to.

Answer 3

Just to point out that if you want to follow @Devon Ryan's answer for a different organism/assembly, that is not in his very useful linked resource, you can download NCBI to UCSC contig to chromosome number mappings from https://www.ncbi.nlm.nih.gov/assembly.

To the site and search for your assembly name. At the bottom of the page is a box called "Global assembly definition" containing a link titled "Download full sequence report".

The downloaded file contains a table with:

• Chromosome numbers in "Sequence-Name"/"Assigned-molecule"
• NCBI names in Refseq-Accn
• UCSC contig name in "UCSC-style-name"