If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

Test-drive the GATK tools and Best Practices pipelines on Terra

Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

CollectHsMetrics for merged bams

Dear all,
I am using PICARD version 2.1.1 to generate metrics for my samples.
The samples were sequenced with llumina NextSeq, we loaded a single library that is a pool of 12 samples.
The flow cell has 4 lanes so at the end of de-multiplexing we have 4 paired fastq files per sample.
I aligned them separately with BWA and then I merged the 4 bams with MarkDuplicates in one step as suggested in your guidelines.
Read groups assigned:

@RG ID:13-2250_L001 LB:L001 PL:illumina SM:13-2250_L001 PU:170129
@RG ID:13-2250_L002 LB:L002 PL:illumina SM:13-2250_L002 PU:170129
@RG ID:13-2250_L003 LB:L003 PL:illumina SM:13-2250_L003 PU:170129
@RG ID:13-2250_L004 LB:L004 PL:illumina SM:13-2250_L004 PU:170129

After the CollectHsMetrics I obtained a MEAN_TARGET_COVERAGE of 165 per sample that is far away from the results expected. We used a high output flow cell so we expected to double (at least) that value.
And in fact, looking at the illumina report generated by base space it is indicated a mean coverage of 400x.
Am I doing something wrong? Is there some step that I am missing?

Any help is really appreciated.




  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Stefania,

    I can't speak to what is reported in BaseSpace (we don't know what metrics are reported or how they're calculated) but let's see what you have here. I'm going to assume that you provided the target intervals list that correspond to your experiment (if not, you should). Have you looked at the rest of the metrics to see if there might be a problem with reads failing filters, distribution of the coverage generated, or some other sequencing problem? Did you look at your duplication metrics?

    Make sure also that you're using the most recent version of Picard to do this. There has been a lot of improvement done in the metrics tools over the past year.

  • Hi @Geraldine_VdAuwera
    thanks for the answer.
    I provided the target intervals to the tool.
    Looking better at my results I found a problem in the bam file using samtools flagstat:

    19529552 + 0 in total (QC-passed reads + QC-failed reads)
    47468 + 0 secondary
    0 + 0 supplementary
    0 + 0 duplicates
    19470160 + 0 mapped (99.70% : N/A)
    19482084 + 0 paired in sequencing
    9741042 + 0 read1
    9741042 + 0 read2
    0 + 0 properly paired (0.00% : N/A)
    19422692 + 0 with itself and mate mapped
    0 + 0 singletons (0.00% : N/A)
    145424 + 0 with mate mapped to a different chr
    113 + 0 with mate mapped to a different chr (mapQ>=5)

    I suspect that the problem is related with "0 + 0 properly paired (0.00% : N/A)".
    I am now performing other tests and controls to verify where is the problem.
    I will let you know.

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Are the @RG groups you list all for the same sample?


    If so, then the sample names should be identical.

  • Hi @shlee and @Geraldine_VdAuwera,

    I aligned bam again and it seems the results are corrected now.
    As procedure I aligned per-lane bam with bwa and then I used AddOrReplaceReadGroup. Then I merged bam files using MarkDuplicates.
    The RGs are:

    @RG ID:16-1377_L001 LB:13-2250 PL:illumina SM:16-1377 PU:170129
    @RG ID:16-1377_L002 LB:13-2250 PL:illumina SM:16-1377 PU:170129
    @RG ID:16-1377_L003 LB:13-2250 PL:illumina SM:16-1377 PU:170129
    @RG ID:16-1377_L004 LB:13-2250 PL:illumina SM:16-1377 PU:170129

    I checked the number of reads with samtools flagstat and in the merged file I have less reads than in the separated bam, about 15 million reads less. In your experience is this expected or could be related with the merge?

    Thanks in advance,


  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin
    It seems odd to me that you would lose reads in the merge, unless you're instructing MarkDuplicates to remove duplicates in your command. Can you post the command line you used for the merge?

    Also, have you looked at whether you're losing reads from one of the readgroups or from several, or from a particular config or chromosome?
  • Hi @Geraldine_VdAuwera
    the command I used is:

    java -jar /usr/local/cluster/bin/picard.jar MarkDuplicates

    As you supposed, I am instructing MarkDuplicates to remove duplicated reads.
    The duplicated reads are about 16% of the total reads.
    Here the schema:

    RG Merge file Single BAM
    16-1377_L001 20650653 24560568
    16-1377_L002 20412747 24086065
    16-1377_L003 20639487 24674281
    16-1377_L004 20425025 24035773

    Is 16% too high or is it an expected result?
    Is there any threshold of acceptability for duplicate reads?
    Thanks for the kind support


    Issue · Github
    by Sheila

    Issue Number
    Last Updated
    Closed By
  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @merella.stefania,

    16% can be reasonable depending on the technology you are sequencing with. Take a look at for more discussion on factors to consider when marking duplicates.

Sign In or Register to comment.