The current GATK version is 3.7-0
Examples: Monday, today, last week, Mar 26, 3/26/04

Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

Get notifications!


You can opt in to receive email notifications, for example when your questions get answered or when there are new announcements, by following the instructions given here.

Did you remember to?


1. Search using the upper-right search box, e.g. using the error message.
2. Try the latest version of tools.
3. Include tool and Java versions.
4. Tell us whether you are following GATK Best Practices.
5. Include relevant details, e.g. platform, DNA- or RNA-Seq, WES (+capture kit) or WGS (PCR-free or PCR+), paired- or single-end, read length, expected average coverage, somatic data, etc.
6. For tool errors, include the error stacktrace as well as the exact command.
7. For format issues, include the result of running ValidateSamFile for BAMs or ValidateVariants for VCFs.
8. For weird results, include an illustrative example, e.g. attach IGV screenshots according to Article#5484.
9. For a seeming variant that is uncalled, include results of following Article#1235.

Did we ask for a bug report?


Then follow instructions in Article#1894.

Formatting tip!


Wrap blocks of code, error messages and BAM/VCF snippets--especially content with hashes (#)--with lines with three backticks ( ``` ) each to make a code block as demonstrated here.

Jump to another community
Picard 2.9.4 is now available. Download and read release notes here.
GATK 3.7 is here! Be sure to read the Version Highlights and optionally the full Release Notes.

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.

Thanks,

Stefania

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    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.
    Thanks!!

  • shleeshlee CambridgeMember, Broadie, Moderator

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

    SM:13-2250_L001
    SM:13-2250_L002
    SM:13-2250_L003
    SM:13-2250_L004
    

    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,

    Stefania

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    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
    INPUT=16-1377/16-1377_L001.RG.bam
    INPUT=16-1377_L002.RG.bam
    INPUT=16-1377/16-1377_L003.RG.bam
    INPUT=16-1377_L004.RG.bam
    OUTPUT=16-1377_RGdedup.bam
    METRICS_FILE=16-1377_metrics
    REMOVE_DUPLICATES=true
    ASSUME_SORTED=true

    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

    Stefania

    Issue · Github
    by Sheila

    Issue Number
    1763
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    sooheelee
  • shleeshlee CambridgeMember, Broadie, Moderator

    Hi @merella.stefania,

    16% can be reasonable depending on the technology you are sequencing with. Take a look at http://gatkforums.broadinstitute.org/gatk/discussion/6747/how-to-mark-duplicates-with-markduplicates-or-markduplicateswithmatecigar for more discussion on factors to consider when marking duplicates.

Sign In or Register to comment.