(How to) Fix a badly formatted BAM

Geraldine_VdAuweraGeraldine_VdAuwera Posts: 9,977Administrator, Dev admin
edited December 2015 in Tutorials

Fix a BAM that is not indexed or not sorted, has not had duplicates marked, or is lacking read group information. The options on this page are listed in order of decreasing complexity.

You may ask, is all of this really necessary? The GATK imposes strict formatting guidelines, including requiring certain read group information, that other software packages do not require. Although this represents a small additional processing burden upfront, the downstream benefits are numerous, including the ability to process library data individually, and significant gains in speed and parallelization options.

Prerequisites

  • Installed Picard tools
  • If indexing or marking duplicates, then coordinate sorted reads
  • If coordinate sorting, then reference aligned reads
  • For each read group ID, a single BAM file. If you have a multiplexed file, separate to individual files per sample.

Jump to a section on this page

  1. Add read groups, coordinate sort and index using AddOrReplaceReadGroups
  2. Coordinate sort and index using SortSam
  3. Index an already coordinate-sorted BAM using BuildBamIndex
  4. Mark duplicates using MarkDuplicates

Tools involved

Related resources

  • Our dictionary entry on read groups discusses the importance of assigning appropriate read group fields that differentiate samples and factors that contribute to batch effects, e.g. flow cell lane. Be sure that your read group fields conform to the recommendations.
  • Picard's standard options includes adding CREATE_INDEX to the commands of any of its tools that produce coordinate sorted BAMs.


1. Add read groups, coordinate sort and index using AddOrReplaceReadGroups

Use Picard's AddOrReplaceReadGroups to appropriately label read group (@RG) fields, coordinate sort and index a BAM file. Only the five required @RG fields are included in the command shown. Consider the other optional @RG fields for better record keeping.

java -jar picard.jar AddOrReplaceReadGroups \ 
    INPUT=reads.bam \ 
    OUTPUT=reads_addRG.bam \ 
    RGID=H0164.2 \ #be sure to change from default of 1
    RGLB= library1 \ 
    RGPL=illumina \ 
    RGPU=H0164ALXX140820.2 \ 
    RGSM=sample1 \ 

This creates a file called reads_addRG.bam with the same content and sorting as the input file, except the SAM record header's @RG line will be updated with the new information for the specified fields and each read will now have an RG tag filled with the @RG ID field information. Because of this repetition, the length of the @RG ID field contributes to file size.

To additionally coordinate sort by genomic location and create a .bai index, add the following options to the command.

    SORT_ORDER=coordinate \ 
    CREATE_INDEX=true

To process large files, also designate a temporary directory.

    TMP_DIR=/path/shlee #sets environmental variable for temporary directory


2. Coordinate sort and index using SortSam

Picard's SortSam both sorts and indexes and converts between SAM and BAM formats. For coordinate sorting, reads must be aligned to a reference genome.

java -jar picard.jar SortSam \ 
    INPUT=reads.bam \ 
    OUTPUT=reads_sorted.bam \ 
    SORT_ORDER=coordinate \

Concurrently index by tacking on the following option.

    CREATE_INDEX=true

This creates a file called reads_sorted.bam containing reads sorted by genomic location, aka coordinate, and a .bai index file with the same prefix as the output, e.g. reads_sorted.bai, within the same directory.

To process large files, also designate a temporary directory.

    TMP_DIR=/path/shlee #sets environmental variable for temporary directory


3. Index an already coordinate-sorted BAM using BuildBamIndex

Picard's BuildBamIndex allows you to index a BAM that is already coordinate sorted.

java -jar picard.jar BuildBamIndex \ 
    INPUT=reads_sorted.bam 

This creates a .bai index file with the same prefix as the input file, e.g. reads_sorted.bai, within the same directory as the input file. You want to keep this default behavior as many tools require the same prefix and directory location for the pair of files. Note that Picard tools do not systematically create an index file when they output a new BAM file, whereas GATK tools will always output indexed files.

To process large files, also designate a temporary directory.

    TMP_DIR=/path/shlee #sets environmental variable for temporary directory


4. Mark duplicates using MarkDuplicates

Picard's MarkDuplicates flags both PCR and optical duplicate reads with a 1024 (0x400) SAM flag. The input BAM must be coordinate sorted.

java -jar picard.jar MarkDuplicates \ 
    INPUT=reads_sorted.bam \ 
    OUTPUT=reads_markdup.bam \
    METRICS_FILE=metrics.txt \
    CREATE_INDEX=true

This creates a file called reads_markdup.bam with duplicate reads marked. It also creates a file called metrics.txt containing duplicate read data metrics and a .bai index file.

To process large files, also designate a temporary directory.

    TMP_DIR=/path/shlee #sets environmental variable for temporary directory
  • During sequencing, which involves PCR amplification within the sequencing machine, by a stochastic process we end up sequencing a proportion of DNA molecules that arose from the same parent insert. To be stringent in our variant discovery, GATK tools discount the duplicate reads as evidence for or against a putative variant.
  • Marking duplicates is less relevant to targeted amplicon sequencing and RNA-Seq analysis.
  • Optical duplicates arise from a read being sequenced twice as neighboring clusters.

back to top


Post edited by shlee on

Geraldine Van der Auwera, PhD

Comments

  • HeidihuangHeidihuang ChinaPosts: 18Member

    If I want to add more header information for example, HD, VN, GO,SO, SQ,SN,LN,RG,ID,LB,SM,CL how to do that

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 9,977Administrator, Dev admin

    To add or modify read group information, use Picard's AddOrReplaceReadGroup. For other header information you'll have to look at other tools in Picard or samtools. GATK does not provide functions to do this.

    Geraldine Van der Auwera, PhD

  • roger_droger_d Redwood CityPosts: 3Member

    This was very helpful to me, but I needed to make a few changes to get the examples to work:

    In 2) the required input METRICS_FILE should be added, like this:
    java -jar MarkDuplicates.jar \
    INPUT=sorted_reads.bam \
    OUTPUT=dedup_reads.bam \
    METRICS_FILE=sorted_reads.dup

    in 4) BuildBamIndex should be BuildBamIndex.jar

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 9,977Administrator, Dev admin

    Thanks for pointing this out, @roger_d. I've made the necessary corrections.

    Note for others: if you are using a very recent version of Picard, you may need to adapt the command lines to the new Picard syntax, which uses one monolithic jar to which you pass tool names (e.g. -jar picard.jar BuildBamIndex etc.), like we have in GATK.

    Geraldine Van der Auwera, PhD

  • shleeshlee CambridgePosts: 84Member, Broadie, Moderator, Dev admin

    I've recently revamped this page and it is pending review.

  • AmandaAmanda North CarolinaPosts: 8Member

    Has there been a change in the options for adding read groups? I am getting an error where it says it does not recognize the option RGID. Thank you!

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 9,977Administrator, Dev admin

    I can't think of anything major like that. Can you please post your full command and console output? And version too.

    Geraldine Van der Auwera, PhD

  • AmandaAmanda North CarolinaPosts: 8Member

    Hi Geraldine, Please disregard, I found a typo that had been integrated somehow into a saved script that was causing the error. Thank you!

Sign In or Register to comment.