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!

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

#### ☞ 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.

GATK 3.7 is here! Be sure to read the Version Highlights and optionally the full Release Notes.

# CollectHsMetrics for merged bams

ItalyMember Posts: 28

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.

@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

Tagged:

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.

Geraldine Van der Auwera, PhD

• ItalyMember Posts: 28

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

47468 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
19470160 + 0 mapped (99.70% : N/A)
19482084 + 0 paired in sequencing
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!!

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.

• ItalyMember Posts: 28

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?

Stefania

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?

Geraldine Van der Auwera, PhD

• ItalyMember Posts: 28

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.
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 Number
1763
State
closed
Last Updated
Assignee
Array
Milestone
Array
Closed By
sooheelee