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

vanish007vanish007 Cleveland, OHMember

Hi everyone,

Thank you in advance for reading this. Has anyone else had the following problem when running CollectHSmetrics:

[Tue Apr 12 12:17:52 EDT 2016] picard.analysis.directed.CollectHsMetrics done. Elapsed time: 14.54 minutes.
Runtime.totalMemory()=920649728
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
Exception in thread "main" java.lang.ArrayIndexOutOfBoundsException: -1
at picard.analysis.directed.TargetMetricsCollector$PerUnitTargetMetricCollector.calculateTargetCoverageMetrics(TargetMetricsCollector.java:637)
at picard.analysis.directed.TargetMetricsCollector$PerUnitTargetMetricCollector.finish(TargetMetricsCollector.java:577)
at picard.metrics.MultiLevelCollector$AllReadsDistributor.finish(MultiLevelCollector.java:208)
at picard.metrics.MultiLevelCollector.finish(MultiLevelCollector.java:324)
at picard.analysis.directed.CollectTargetedMetrics.doWork(CollectTargetedMetrics.java:152)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:209)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:95)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:105)

The Bait and Target I used was just the header of my BAM file pulled by Samtools view header.

I saw something in the FAQ saying that, "The BAM index format imposes a limit of 512MB for a single reference sequence. If this limit is exceeded, various errors may occur depending on what steps have been taken."

Any explanation on how to address this? Thank you again!

Best Answer

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @vanish007
    Hi,

    What kind of reference are you working with? How many contigs does it have?

    Are you using the latest version of Picard?

    Thanks,
    Sheila

  • vanish007vanish007 Cleveland, OHMember

    HI Sheila,

    Thank you for answering! I'm using BAM files downloaded from the TCGA data matrix. I'm not sure how many contigs these files have nor do I know how to get that info.
    I am using the Picardtools 2.1.1 so it's fairly new.
    Currently I am trying to generate a BED file from the BAM file and then change that BED to an intervals file, though I don't know if this will solve the problem or not... Sorry I'm new at this.

    Best,

    Vai

  • vanish007vanish007 Cleveland, OHMember

    Actually never mind...I guess I don't really have a sequence dictionary since the files from the TCGA are just Whole Exome sequencing BAM files and not FastQ files...

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @vanish007
    Hi Vai,

    What is the exact command you are running? Are you trying to use this tool to get a bed/intervals file? That is not what the tool does. You may look into the TCGA site to find a pre-made intervals file.

    -Sheila

  • vanish007vanish007 Cleveland, OHMember

    Hi Sheila,

    Yes, I am trying to run the HSmetrics and had no idea where or how to obtain interval files for the BAIT and TARGET properties.

    I essentially just used the Samtools "view header" command to get a header from the BAM file and used that as my BAIT and TARGET option, but that resulted in the error above. Therefore I was seeking out a way to create an Interval file of my own. Is trying to find a pre-made interval the only way to get the collecthsmetrics tool to work? Thanks again for your patience in all this!

    Best,

    Vai

  • vanish007vanish007 Cleveland, OHMember

    Hi Sheila,

    Awesome, thank you so much for answering! That helps clear things up quite a bit. Have a great rest of your week! =)

    Best,

    Vai

  • RashRash Member

    @Sheila Dear Sheila, I also have the same error.
    I did try to build the bait and target interval (using one time the headers from .bam file and one time with headers from mm10.dict) like this:
    m=0
    for i in 'cat mm10_list.interval_list' ; \
    do echo "$i"; echo "MPMPMP$m"; m=$(($m + 1)) ; \
    done > mm10_list.interval_list.tmp2
    perl -0777 -i -pe "s/\nMPMPMP/\ttarget_/g" mm10_list.interval_list.tmp2
    perl -0777 -i -pe "s/\r\t/\t/g" mm10_list.interval_list.tmp2
    perl -0777 -i -pe "s/:/\t/g" mm10_list.interval_list.tmp2
    perl -0777 -i -pe "s/-/\t/g" mm10_list.interval_list.tmp2
    cat header.txt mm10_list.interval_list.tmp2 > mm10_list.interval_list.picard

    but I have still error while running calculateHSmetrics. Would you please let me know if the way I am building the bait and target files is correct? and I am using the same file for bait_interval and also target_interval, is this also correct?
    Many thanks in advance,
    Rahel

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Rash
    Hi Rahel,

    As mentioned in my previous post, the bam file header is not a valid intervals list. Where did you get your BAM file from? You should ask the sequencing company/sequence provider to give you the intervals list.

    -Sheila

  • RashRash Member
    edited November 2016

    Dear sheila @Sheila ,
    Many thanks for your response. My BAM file is mapped.realigned.recal.bam that I get from "PrintReads" arguments. Now I could do HsMetrics with exon.mm10.interval.list from UCSC. To do so, I got first exon_mm10_target.intervals.bed from UCSC and then I used "BedToIntervalList" to convert it to interval list and then I used for CollectHsMetrics. Is this correct? Many thanks in advance for your time,
    Rahel

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Rash
    Hi Rahel,

    Yes, that is fine. Did everything work alright?

    -Sheila

  • RashRash Member

    @Sheila Dear Sheila,
    Yes, it works well without error. but in the the Hsmetrics ourput, I have 0 for all HS_PENALTY (*X) and in the histogram part, after some numbers, I just have 0 (in baseq_count column).
    I add it here, is that correct or normal output? Or something went wrong? Many many thanks for your help, Rahel

    METRICS CLASS picard.analysis.directed.HsMetrics

    BAIT_SET GENOME_SIZE BAIT_TERRITORY TARGET_TERRITORY BAIT_DESIGN_EFFICIENCY TOTAL_READS PF_READS PF_UNIQUE_READS PCT_PF_READS PCT_PF_UQ_READS PF_UQ_READS_ALIGNED PCT_PF_UQ_READS_ALIGNED PF_BASES_ALIGNED PF_UQ_BASES_ALIGNED ON_BAIT_BASES NEAR_BAIT_BASES OFF_BAIT_BASES ON_TARGET_BASES PCT_SELECTED_BASES PCT_OFF_BAIT ON_BAIT_VS_SELECTED MEAN_BAIT_COVERAGE MEAN_TARGET_COVERAGE MEDIAN_TARGET_COVERAGE PCT_USABLE_BASES_ON_BAIT PCT_USABLE_BASES_ON_TARGET FOLD_ENRICHMENT ZERO_CVG_TARGETS_PCT PCT_EXC_DUPE PCT_EXC_MAPQ PCT_EXC_BASEQ PCT_EXC_OVERLAP PCT_EXC_OFF_TARGET FOLD_80_BASE_PENALTY PCT_TARGET_BASES_1X PCT_TARGET_BASES_2X PCT_TARGET_BASES_10X PCT_TARGET_BASES_20X PCT_TARGET_BASES_30X PCT_TARGET_BASES_40X PCT_TARGET_BASES_50X PCT_TARGET_BASES_100X HS_LIBRARY_SIZE HS_PENALTY_10X HS_PENALTY_20X HS_PENALTY_30X HS_PENALTY_40X HS_PENALTY_50X HS_PENALTY_100X AT_DROPOUT GC_DROPOUT HET_SNP_SENSITIVITY HET_SNP_Q SAMPLE LIBRARY READ_GROUP
    hg38_list 3099750718 83532717 83532717 1 105776794 105776794 105776794 1 1 105776794 1 15599449101 15599449101 11347371082 4252078019 0 11347371082 1 0 0.727421 135.843433 135.843433 20 0.727391 0.727391 26.993311 0.034427 0 0 0 0 0.272579 67.921717 0.885454 0.823709 0.53052 0.500379 0.48641 0.474407 0.463272 0.409135 0 0 0 0 0 0 7.25602 2.402312 0.726945 6

    HISTOGRAM java.lang.Integer

    coverage count baseq_count
    0 11514921 0
    1 5157722 80325
    2 7038330 481773
    ..
    ..
    37 97755 1998330137
    38 97301 1411706604
    39 95360 0
    40 95007 0
    41 93775 0
    42 94352 0
    ..
    ..
    197 97736 0
    198 97066 0
    199 97768 0
    200 24539388 0

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Rash
    Hi Rahel,

    Have you solved your issue? I think perhaps some of the other posts may have helped?

    -Sheila

  • RashRash Member
    edited November 2016

    @Sheila Dear Sheila, I have still problem to get the right interval list. I used mm10 (GRCm38) and GRCh38 from Ensemble for mapping. But I get the exome.intervals.bed with 10bp before and after exome area from UCSC table browser (using GRCm38 and GRCh38). Then, I did convert it to interval list using Picard. But it seems that the intervals are not correct. I also have many "./." in the final combined VCF output. Would you please help me with finding the correct interval list for the genomes I am using in this study (Considering that I do not have access to the sequencing center!). Many thanks in advance, Rahel

    Issue · Github
    by Sheila

    Issue Number
    1468
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    chandrans
  • RashRash Member

    @Sheila Dear sheila, I finally could get the interval file from illumina which it is for hg19 in bed format. I did first convert it to bam and after that I used liftoverintervals from picard to make it suitable for mm10 and hg38 but I have this errors, would you please help me with this problem. Many thanks in advance, Rahel

    java -Xmx4g -jar picard.jar LiftOverIntervalList \
    I=hg19_illumina_nochr.interval_list \
    O=mm10_illumina_list.interval_list \
    SD=Mus_musculus.GRCm38.dna.primary_assembly.dict \
    CHAIN=hg19ToMm10.over.chain
    picard.util.LiftOverIntervalList INPUT=hg19_illumina_nochr.interval_list OUTPUT=mm10_illumina_list.interval_list SEQUENCE_DICTIONARY=Mus_musculus.GRCm38.dna.primary_assembly.dict CHAIN=hg19ToMm10.over_nochr.chain MIN_LIFTOVER_PCT=0.95 VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=false CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json
    [Thu Nov 24 18:10:54 CET 2016] Executing as [email protected] on Linux 2.6.32-504.12.2.el6.x86_64 amd64; Java HotSpot(TM) 64-Bit Server VM 1.8.0_111-b14; Picard version: 2.5.0(2c370988aefe41f579920c8a6a678a201c5261c1_1466708365)
    [Thu Nov 24 18:11:07 CET 2016] picard.util.LiftOverIntervalList done. Elapsed time: 0.22 minutes.
    Runtime.totalMemory()=2620391424
    To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
    Exception in thread "main" htsjdk.samtools.SAMException: Reached end of chain without seeing terminal block in chain file hg19ToMm10.over_nochr.chain at line 34330636
    at htsjdk.samtools.liftover.Chain.throwChainFileParseException(Chain.java:412)
    at htsjdk.samtools.liftover.Chain.loadChain(Chain.java:386)
    at htsjdk.samtools.liftover.Chain.loadChains(Chain.java:318)
    at htsjdk.samtools.liftover.LiftOver.(LiftOver.java:62)
    at picard.util.LiftOverIntervalList.doWork(LiftOverIntervalList.java:101)
    at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:208)
    at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:95)
    at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:105)

  • RashRash Member

    @Sheila I solved the problem by removing the chr and extra contigs from .chain file but I have the same result for Hsmetrics!!! and for Depth of coverage, I just have these four output (attached). Do you have any idea where is the problem? Can I use any of the output from DepthofCoverage to run XHMM? Many thanks in advance,

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Rash
    Hi Rahel,

    Okay. So, you have an input interval file that works with one Picard tool but not another? And, that same interval file works with GATK? Can you post the exact command and error you get with CollectHsMetrics? What do you mean by "run XHMM"?

    Thanks,
    Sheila

  • RashRash Member

    @Sheila Hi Sheila, Many thanks for your reply. I got an interval list from illumina for hg19 and then by using "LiftOverIntervalList" I made it suitable for hg38 and also mm10 (I have data from mouse and also human). I have solved the error of this step and There is no error in this step.
    After that I used this 2nd interval list for collecting HsMetrics and I did not get any error here too. But in the output, I have many zero in the part for HISTOGRAM. I just want to be sure if there is a problem or it is normal to have such an output?

    HISTOGRAM java.lang.Integer
    coverage count baseq_count
    0 11514921 0
    1 5157722 80325
    2 7038330 481773
    ..
    ..
    37 97755 1998330137
    38 97301 1411706604
    39 95360 0
    40 95007 0
    41 93775 0
    42 94352 0
    ..
    ..
    197 97736 0
    198 97066 0
    199 97768 0
    200 24539388 0

    And by XHMM (http://atgu.mgh.harvard.edu/xhmm/tutorial.shtml), I wanted to know which output from DepthofCovergae is suitable for it and I found out that the interval.summary is the correct one.

    Many thanks for your help, Rahel

    Issue · Github
    by Sheila

    Issue Number
    1512
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    chandrans
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
    edited December 2016

    @Rash
    Hi Rahel,

    It turns out baseq_count refers to the actual number of bases that have the base quality in the first column. The second and third column are completely independent of one another.

    In your case, you have 0's in the baseq_count for base qualities above 38 because the sequencer you used does not emit base qualities higher than 38. This is all perfectly acceptable once you understand what the columns mean :smiley:

    -Sheila

    Post edited by Geraldine_VdAuwera on
  • Hello,

    I don't have bait interval files as the company no longer provides it. So I used the target interval file for both BAIT_INTERVALS and TARGET_INTERVALS. But the values of on_bait_bases and on_target_bases are different. I tried to set NEAR_DISTANCE=0 but the results are similar to that using NEAR_DISTANCE=250. Any explanation?

    -Ming

    Issue · Github
    by Sheila

    Issue Number
    2144
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    sooheelee
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @ct6464983
    Hi Ming,

    I am checking with the team and will get back to you.

    -Sheila

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭
    edited June 2017

    Hi @ct6464983,

    ON_BAIT_BASES The number of PF_BASES_ALIGNED that are mapped to the baited regions of the genome.
    ON_TARGET_BASES The number of PF_BASES_ALIGNED that are mapped to a targeted region of the genome.
    NEAR_DISTANCE The maximum distance between a read and the nearest probe/bait/amplicon for the read to be considered 'near probe' and included in percent selected. Default value: 250.

    How different are your values of on_bait_bases and on_target_bases? Can you give us some examples?

    When you say your results are similar when changing the NEAR_DISTANCE parameter, what do you mean? This option is about reads that align a certain distance away from the regions you provide. The reads outside this distance are not counted towards the metrics. The expectation is that off-target alignments will be minimal, no? Is there any reason for you to expect your data to have large off-target effects?

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @ct6464983,
    I followed up on your question with a developer. At first I thought perhaps you may be observing minor differences in rounding (hence my request for example numbers). However, our developer has confirmed that with NEAR_DISTANCE set to other than 0, for your case where you use the same intervals for on-bait and on-target, you could expect different metrics. This is because the NEAR_DISTANCE parameter applies only to the on-bait target intervals. On-target intervals do not use the argument. I will clarify this in the tool parameter description. Thanks for bringing this to our attention.

    Issue · Github
    by shlee

    Issue Number
    2170
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    sooheelee
  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭
    edited June 2017

    Here is some more information from our developer @ct6464983. Other parameters, e.g. MINIMUM_BASE_QUALITY, can explain the differences you are observing when NEAR_DISTANCE is set to zero. They apply to one of the interval sets. Sorry I cannot be more detailed--this is as much as I can do for you. I think the way to test exactly which of these are involved, you can set each parameter to zero and see your two sets of metrics converge.

Sign In or Register to comment.