Heads up:
We’re moving the GATK website, docs and forum to a new platform. Read the full story and breakdown of key changes on this blog.
Notice:
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.

Depth of Coverage: Only first gene summary output

KurtKurt Member ✭✭✭
edited January 2013 in Ask the GATK team

Hello,

I've been trying to run DOC on a sorted RefSeq Gene list downloaded from the UCSC Table Browser as described in the DOC guide doc. After sorting/formatting, only the first record in the ref seq gene list gets summarize in the foo.sample_gene_summary file. I've tried multiple one gene files and they all work individually, but when there is more than one record only the first record gets summarized. I'm using 2.2-8 (tried 2.0-25 as well) and I'm totally flabbergasted. Below is an example command line.

java -jar ./GenomeAnalysisTK-2.2-8-gec077cd/GenomeAnalysisTK.jar \
    -T DepthOfCoverage \
    -R human_g1k_v37_decoy.fasta \
    -I in.bam \
    -geneList:REFSEQ 2genes.txt \
    -L 22:1615633-18714498 \
    -mmq 1 \
    --outputFormat table \
    -o foo \
    -l DEBUG

also attached is the 2genes.txt file (put it at the end as well in case there is problem with the attachment).

Below is a little bit from the logging when DEBUG is switched on (can add more if needed).

INFO  15:13:12,293 RMDTrackBuilder - Loading Tribble index from disk for file /isilon/sequencing/Kurt/Genes/2genes.txt 
DEBUG 15:13:12,364 DepthOfCoverage - Refseq init done. 
DEBUG 15:13:12,364 DepthOfCoverage - Examining 22:1615633-18714498 
DEBUG 15:13:12,365 DepthOfCoverage - Annotation list is anonymous 
DEBUG 15:13:12,366 DepthOfCoverage - We do overlap 000  NM_001136213    22      -       16256331        16287937        16258185        16287885        11      16256331,16258184,16266928,16268136,16269872,16275206,16277747,16279194,16282144,16282477,16287253, 16256677,16258303,16267095,16268181,16269943,16275277,16277885,16279301,16282318,16282592,16287937,      0       POTEH   cmpl    cmpl    -1,2,0,0,1,2,2,0,0,2,0, 
INFO  15:13:12,462 DepthOfCoverage - Printing summary info 
INFO  15:13:12,479 DepthOfCoverage - Printing locus summary 
INFO  15:13:12,772 ProgressMeter -            done        1.71e+07    2.3 m        8.1 s    100.0%         2.3 m     0.0 s 
INFO  15:13:12,773 ProgressMeter - Total runtime 138.10 secs, 2.30 min, 0.04 hours 
INFO  15:13:12,928 MicroScheduler - 808 reads were filtered out during traversal out of 49775 total (1.62%) 
INFO  15:13:12,928 MicroScheduler -   -> 607 reads (1.22% of total) failing DuplicateReadFilter 
INFO  15:13:12,928 MicroScheduler -   -> 201 reads (0.40% of total) failing UnmappedReadFilter 
INFO  15:13:12,929 NSRuntimeProfile - Input   time:   22.5 s (18.59%) 
INFO  15:13:12,929 NSRuntimeProfile - Map     time:   36.2 s (29.93%) 
INFO  15:13:12,929 NSRuntimeProfile - Reduce  time:   59.8 s (49.47%) 
INFO  15:13:12,929 NSRuntimeProfile - Outside time:    2.4 s ( 2.01%) 
DEBUG 15:13:12,931 GATKRunReport - Aggregating data for run report 
DEBUG 15:13:13,112 GATKRunReport - Posting report of type STANDARD 

709     NM_001136213    22      -       16256331        16287937        16258185        16287885        11      16256331,16258184,16266928,16268136,16269872,16275206,16277747,16279194,16282144,16282477,16287253,     16256677,16258303,16267095,16268181,16269943,
16275277,16277885,16279301,16282318,16282592,16287937,  0       POTEH   cmpl    cmpl    -1,2,0,0,1,2,2,0,0,2,0,
90      NM_018943       22      +       18593452        18614498        18593631        18613903        5       18593452,18604245,18606922,18609120,18613609,   18593634,18604468,18607071,18609801,18614498,   0       TUBA8   cmpl    cmpl    0,0,1,0,0,

Is this just user error?

Post edited by Geraldine_VdAuwera on

Best Answer

Answers

  • KurtKurt Member ✭✭✭

    Ah, Thank you very much Geraldine. That's definitely not how I read it (but then again my brain was fried). So I made 2 bed files, one for the transcription start/end intervals for the two genes and then another that had just the exon start/end positions for both genes (one row for each exon), and got the expected results. Thanks again!

  • mariagrmariagr ZMBHMember

    Hello ,
    I ve been trying to write a script for calculating coverage per gene,unsuccessfully(!) ,and I found now that is nicely done by GATK !
    I would very much need to use this calculation of depthOfCoverage for each gene but I cannot find the geneList needed in the format explained here. I have a RefSeq gene list downloaded from UCSC table which contains RefSeq name ,cds_start & end and "chr" information. Is this acceptable? I want to do it for exons falling inside the genes, which I have downloaded from :
    ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/exome_pull_down_targets/ (phase3),
    so this would be my Intervals List. It contains only chr info and start-end.
    Moreover I cannot find in UCSC a table which combines RefSeq name with used gene Name. Is it combined in this genesList you provide from GATK?
    Can you guide me where I could find the exact url for /humgen/.../geneList.txt or if mine could work ,and if exons table is ok with only these 3 columns ?
    I m a registered member, as for writing in the forum. Is there any extra procedure needed to access your database ?
    Thanks in advance !

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @mariagr Unfortunately the gene list referenced above is not currently provided to the external community outside of Broad. Have a look at this FAQ for producing an acceptable RefSeq file: https://www.broadinstitute.org/gatk/guide/article?id=1329

  • seruseru BergenMember ✭✭

    Hi,

    I have discovered that coverage of overlapping genes/exons is not calculated properly by the DoC tool. Coverage stats generated for an interval are attributed only to the first gene on the list that overlaps the interval (similarly to @Kurt 's example). Subsequent genes overlapping the interval are ignored. For reference, please see getGeneName() function in DepthOfCoverage.java. I have made a pull request with a fix here: https://github.com/broadgsa/gatk/pull/16

    I have also a question/comment regarding the way the coverage stats are computed if -geneList is used. As far as I understand the comments above and the code, the walker iterates over intervals provided by the -L argument. If an interval overlaps a gene (i.e. one of its exons), the depth stats of that interval are taken into account when calculating per-gene stats. Is then only a fraction that actually overlaps taken into account, or entire interval? If the second, then the BED file (-L) should contain exact exon coordinates, otherwise coverage of some non-coding regions gets also included, correct?

    A related concern is that, it looks like overlapping intervals in the BED file are merged together before/during traversal. In effect, coverage stats of any overlapping exons (of possibly different genes) are lumped together. With the fix I made, they will be attributed to all overlapping genes, but in exactly the same way. This means that two 1-exon genes: A (900bp long), and B (101bp long) which share 1bp, will have identical coverage stats (being truly stats for the 1kb interval covering A&B). Only workaround I have come up with so far is to chop the exon intervals into pieces, such that each piece covers exactly the same set of exons/genes (the extreme case being 1bp long intervals:)

    Best wishes,
    Paweł

  • seruseru BergenMember ✭✭

    Hi,

    I was just wondering if you had any comments on this? I tested the chopping-idea, and to make it work I had to introduce 1bp breaks between exon/interval fragments in the BED file. Otherwise, adjacent ends were merged together. I guess, the coverage in breaks is not taken into account in the calculation, so the stats are not 100% accurate. On the positive side, there is only 8.9k affected intervals among over 300k (one interval per each exon/CCDS pair), and most of them get introduced only one such 1bp break, so the inaccuracy is marginal.

    Paweł

  • seruseru BergenMember ✭✭

    Ahh, I found --interval_merging OVERLAPPING_ONLY option now:) That removes the problem with merging adjacent intervals, so after all, it should be possible to compute the stats correctly.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @seru
    Hi Pawel,

    Yes, I am happy you found that option! Sorry for the late response. I was thinking about your question, and I thought about the interval_merging option, but then I did not think it solved the pull request you made. I will have someone look at the pull request. Thank you for it :smile:

    -Sheila

  • seruseru BergenMember ✭✭

    Hi Sheila,

    I should have been more clear with my last statement. Sorry. You're right, it does not solve the pull request. The overlapping genes will still need the patch I submitted.

    The --interval_merging setting solved the issue of 1bp breaks I had to introduce in the BED file to stop GATK from merging fragmented exons back again. It is explain in more detail in the script I made for fragmenting the exons:
    https://github.com/seru71/ExomeSeqGenotypingPipeline/blob/diagnostic/reference/b37/chop_exons.py

    Pawel

    Issue · Github
    by Sheila

    Issue Number
    558
    State
    closed
    Last Updated
    Assignee
    Array
    Closed By
    vdauwera
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @seru, we're having someone look at your pull request, thanks for contributing this. Would you have a test case available that we can use to add an integration test?

  • seruseru BergenMember ✭✭

    Hi and sorry for the delay. I am attaching an archive with a short gene file with 3 transcripts for genes MUTYH, TOE1, and TESK2 (test_exons.txt). The file is modified so that the coverage should be reported for every exon (not every gene) separately. The last exon of MUTYH overlaps with the first exon of TOE1. When you run default DoC (test.sh) the first exon of TOE1 should be missing from the gene_summary output.

    Best regards,
    Paweł

  • SheilaSheila Broad InstituteMember, Broadie admin

    @seru
    Hi Paweł,

    Thanks! I will get this to our developer working on this immediately.

    -Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @seru, your patch has been merged into the development repo. The functionality will be available in the next nightly build. Thanks for your contribution!

  • seruseru BergenMember ✭✭

    Great! Thanks.

  • hns04hns04 Member

    Hi @Sheila @Geraldine_VdAuwera I had a question regarding the genelist option. I want the coverage per gene. I created the refseq gene list as described by the GATK team and it works great! But, the intervals file that I have have +-5bp buffer to each exon and so I wanted to be sure if the gene level coverage that DOC will produce will include all the regions from the intervals file?. For example ( interval file has 0 195 but the refseq has 5 190) Will the gene level coverage include bases 0 - 195?.
    Thank you for your help
    Himanshu

  • SheilaSheila Broad InstituteMember, Broadie admin

    @hns04
    Hi Himanshu,

    I just tried this out on a small interval. For example, I tried -L 20:1000000-1000005 -L 20:1000001-1000006. My output file had information for positions 1000000-1000006.

    It looks like overlapping positions are all output as one interval.

    -Sheila

  • hns04hns04 Member

    Thanks for your response @Sheila .. Actually I have the question when you supply the refseq genelist and the intervals, then the output at Gene level includes these additional bases or no?.
    Example ( Refseq file has GENE A chr1 10 - 200 and the -L or interval file has the position chr1 8-205) Will the bases 8-10 and 200-206 be reported in the gene level output?.
    Basically, if our intervals have 5 bp buffer to the coding exons, will we get the information for them in the gene level file?.

    Thanks,
    Himanshu

  • SheilaSheila Broad InstituteMember, Broadie admin

    @hns04
    Hi Himanshu,

    Sorry for the delay. I am about to test this out, so I will get back to you soon.

    -Sheila

  • SheilaSheila Broad InstituteMember, Broadie admin

    @hns04
    Hi Himanshu,

    I just tested this, and it looks like the -L intervals take precedence over the -geneList intervals. So, in your example, yes, the extra padding in the -L intervals will be included in the output. However, I should note, that if you have intervals in your -geneList interval that are not in the -L interval, those will not be included in the output.

    -Sheila

  • lalirlalir CanadaMember

    Hello @Geraldine_VdAuwera and @Sheila I just have a follow up question for this thread. If I am using the -genelist argument, will the per gene statistics be based on the intervals specified in the -L argument? For instance, will the per gene statistics be based on the intervals supplied in the -L argument?

    Also, if there are multiple transcripts per gene in the refgene.sorted file, for which transcript are the gene level statistics calculated?

    Thank you for your assistance.

    Issue · Github
    by Sheila

    Issue Number
    1334
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    vdauwera
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hey @lalir, I'm pretty sure the answer to your first question is the post that is highlighted in green at the top of this thread.

    Regarding transcripts, I believe only the first is used.

Sign In or Register to comment.