No mutect output for certain intervals - Any idea why?

AnushiAnushi SydneyMember

Hi,

I ran mutect on certain specific intervals e.g :- chr1:40505653-40505653, chr1:139897149-139897149, chr1:198760696-198760696, etc.
However, it does not print any output for such intervals (whether it is REJECT or KEEP).
I also checked my tumor bam file to see if any genotype exists using "samtools mpileup". And it does show the genotype as well as mutation for those intervals.
Any idea why could the reason for this? I would really appreciate some information on this.

Thanks,
Anushi

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @Anushi,

    Have you checked the reference sequence at those positions for N bases? If a reference base is N, the program will skip it.

  • AnushiAnushi SydneyMember

    Thanks Geraldine for your reply. I checked and the reference base is not N for those positions. Do you think there could be some other reason for this?

    Anushi

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    It's difficult to say without seeing the data. Could you maybe post a screenshot of the region in IGV?

  • AnushiAnushi SydneyMember

    Hi Geraldine,

    I trying to paste the screen shot from IGV using image button above. (I give the path to the image file located on my computer).
    But it does not show up. Is there some other way to do it?

    Thanks a lot for your help,
    Anushi

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Ah, the image button only works if your image is online somewhere. You can use a sharing service like dropbox to share image links easily. In the main post it's possible to attach an image directly, but not in comments, unfortunately.

  • AnushiAnushi SydneyMember
    edited November 2014

    Oh, I see. thanks for the information. I just shared the image from dropbox. Here's the link, just in case you cannot view the pasted image below. https://www.dropbox.com/s/a7vlmvb3uufh4gl/igv_snapshot.png?dl=0
    The image shows the tumor and normal bam files in IGV around region 40505653 in chr1.

    image

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @Anushi,

    I see, thanks for posting the picture. This does look like a reasonable site, I'm not sure why it's not getting output. To be clear, do some of the sites get output and not others, or is none of them output? And can you please post your full command line and log output please?

  • AnushiAnushi SydneyMember

    Yes, some of the sites get output and some don't.

    The command that I use is :

    java -Xmx5g -jar /home/ays561/softwares/muTect-1.1.4-bin/muTect-1.1.4.jar \
    --analysis_type MuTect \
    --reference_sequence /home/ays561/sequence_alignment/data/ucsc.hg19.fasta \
    --cosmic /home/ays561/sequence_alignment/data/hg19_cosmic_v54_120711.vcf \
    --dbsnp /mnt/rdsi/tcga/TCGA/cghub_data/anushi/data/dbsnp_137.hg19.vcf \
    --intervals "chr1:40505653-40505653" \
    --input_file:normal /mnt/rdsi/tcga/TCGA/cghub_data/mutect_test1/normal_bw2_uniq_0.bam \
    --input_file:tumor /mnt/rdsi/tcga/TCGA/cghub_data/mutect_test1/tumour_bw2_uniq_0.bam \
    --out /mnt/rdsi/tcga/TCGA/cghub_data/mutect_test1/call_stats.out \
    --coverage_file /mnt/rdsi/tcga/TCGA/cghub_data/mutect_test1/coverage.wig.txt \
    -nt 4 \
    --vcf /mnt/rdsi/tcga/TCGA/cghub_data/mutect_test1/out.vcf

    Output log :

    INFO 11:29:38,071 HelpFormatter - ---------------------------------------------------------------------------------
    INFO 11:29:38,073 HelpFormatter - The Genome Analysis Toolkit (GATK) v2.2-25-g2a68eab, Compiled 2012/11/08 10:30:02
    INFO 11:29:38,073 HelpFormatter - Copyright (c) 2010 The Broad Institute
    INFO 11:29:38,073 HelpFormatter - For support and documentation go to http://www.broadinstitute.org/gatk
    INFO 11:29:38,077 HelpFormatter - Program Args: --analysis_type MuTect --reference_sequence /home/ays561/sequence_alignment/data/ucsc.hg19.fasta --cosmic /home/ays561/sequence_alignment/data/hg19_cosmic_v54_120711.vcf --dbsnp /mnt/rdsi/tcga/TCGA/cghub_data/anushi/data/dbsnp_137.hg19.vcf --intervals chr1:40505653-40505653 --input_file:normal /mnt/rdsi/tcga/TCGA/cghub_data/mutect_test1/normal_bw2_uniq_0.bam --input_file:tumor /mnt/rdsi/tcga/TCGA/cghub_data/mutect_test1/tumour_bw2_uniq_0.bam --out /mnt/rdsi/tcga/TCGA/cghub_data/mutect_test1/call_stats.out --coverage_file /mnt/rdsi/tcga/TCGA/cghub_data/mutect_test1/coverage.wig.txt -nt 4 --vcf /mnt/rdsi/tcga/TCGA/cghub_data/mutect_test1/out.vcf
    INFO 11:29:38,077 HelpFormatter - Date/Time: 2014/11/10 11:29:38
    INFO 11:29:38,077 HelpFormatter - ---------------------------------------------------------------------------------
    INFO 11:29:38,077 HelpFormatter - ---------------------------------------------------------------------------------
    INFO 11:29:38,105 ArgumentTypeDescriptor - Dynamically determined type of /mnt/rdsi/tcga/TCGA/cghub_data/anushi/data/dbsnp_137.hg19.vcf to be VCF
    INFO 11:29:38,111 ArgumentTypeDescriptor - Dynamically determined type of /home/ays561/sequence_alignment/data/hg19_cosmic_v54_120711.vcf to be VCF
    INFO 11:29:38,119 GenomeAnalysisEngine - Strictness is SILENT
    INFO 11:29:38,235 GenomeAnalysisEngine - Downsampling Settings: Method: BY_SAMPLE Target Coverage: 1000
    INFO 11:29:38,251 SAMDataSource$SAMReaders - Initializing SAMRecords in serial
    INFO 11:29:38,600 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 0.35
    INFO 11:29:38,625 RMDTrackBuilder - Loading Tribble index from disk for file /mnt/rdsi/tcga/TCGA/cghub_data/anushi/data/dbsnp_137.hg19.vcf
    INFO 11:29:39,040 RMDTrackBuilder - Loading Tribble index from disk for file /home/ays561/sequence_alignment/data/hg19_cosmic_v54_120711.vcf
    INFO 11:29:39,114 GenomeAnalysisEngine - Processing 1 bp from intervals
    INFO 11:29:39,119 MicroScheduler - Running the GATK in parallel mode with 4 total threads, 1 CPU thread(s) for each of 4 data thread(s), of 16 processors available on this machine
    INFO 11:29:39,127 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
    INFO 11:29:39,127 ProgressMeter - Location processed.sites runtime per.1M.sites completed total.runtime remaining
    INFO 11:29:39,945 Walker - [REDUCE RESULT] Traversal result is: 0
    INFO 11:29:39,958 ProgressMeter - done 1.00e+00 0.8 s 9.6 d 100.0% 0.8 s 0.0 s
    INFO 11:29:39,959 ProgressMeter - Total runtime 0.83 secs, 0.01 min, 0.00 hours
    INFO 11:29:40,117 MicroScheduler - 0 reads were filtered out during traversal out of 29 total (0.00%)
    INFO 11:29:40,118 NSRuntimeProfile - Input time: 0.0 s ( 0.07%)
    INFO 11:29:40,118 NSRuntimeProfile - Map time: 0.1 s ( 2.57%)
    INFO 11:29:40,119 NSRuntimeProfile - Reduce time: 0.0 s ( 0.00%)
    INFO 11:29:40,120 NSRuntimeProfile - Outside time: 3.9 s (97.36%)
    WARN 11:29:41,953 RestStorageService - Error Response: PUT '/GATK_Run_Reports/Mx2GXXCC07gn1tbYbFOPxUSspry08jfU.report.xml.gz' -- ResponseCode: 403, ResponseStatus: Forbidden, Request Headers: [Content-Length: 398, Content-MD5: 1ebwYx4Cng01cT16ozc+hg==, Content-Type: application/octet-stream, x-amz-meta-md5-hash: d5e6f0631e029e0d35713d7aa3373e86, Date: Mon, 10 Nov 2014 00:29:40 GMT, Authorization: AWS AKIAJXU7VIHBPDW4TDSQ:rkQGqGOxvy3VJDtvu562FBuYip8=, User-Agent: JetS3t/0.8.1 (Linux/3.0.26-0.7-default; amd64; en; JVM 1.7.0_51), Host: s3.amazonaws.com, Expect: 100-continue], Response Headers: [x-amz-request-id: 07B4D25F9291DBD6, x-amz-id-2: 7qVsEHuAxLG4rjJxv29gSBNuayUsGd7zPB6D0vZme4//t4xI//DJRiD70kzSQwGe, Content-Type: application/xml, Transfer-Encoding: chunked, Date: Mon, 10 Nov 2014 00:29:41 GMT, Connection: close, Server: AmazonS3]

    Thanks,

    Anushi

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Does the same thing happen if you remove -nt 4 from your command"?

  • AnushiAnushi SydneyMember

    Thanks Geraldine for your suggestion. I tried by removing "-nt 4" option. Unfortunately there's no output yet. In another experiment, I also tried removing --cosmic and/or --dbsnp option. Still no output.

    Anushi

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @Anushi, sorry for the response delay. Can you please try running the GATK HaplotypeCaller on this site, with -allsites in your command line?

Sign In or Register to comment.