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.

gender and/or R problem ?

Will_GilksWill_Gilks University of Sussex, UKMember ✭✭
edited December 2015 in GenomeSTRiP

Hi,

My deletion calling pipeline runs to completion but there are error messages (see below). The first error message points to SVDiscovery-1.out
which contains the error message on .R: Error in chisq.test(m) : all entries of 'x' must be nonnegative and finite. In this case, it looks like the problem is negative values in the first cluster... so why could these be generated..?

The gender report lists dosage of X and Y as zero for all my 16 samples. I made a custom ploidy info file for D.melanogaster so maybe this file is wrong. In the reference sequence, there are unmapped scaffolds which have X or Y assignment, so in the ploidy info file, I indicated that they had the same dosage as X and Y for male and female.

The first possible oddity in the output is that the individual depth files all have Frags_per_base as Infinity. No discovery.vcf or genotypes.vcf are generated. Just unfiltered.vcf, although this contains no variants (just headers).I'm considering removing all unmapped scaffold/chromosome data from the bams, and from the reference, then running genomestrip but it seems like a bit of a drag. I don't know if the two error messages are related to the sex chromosome dosage issue. Also we had an issue with R compatibility for plotting regional graphs - I guess it's not critical, or is it for other calculations?
Any advice would be appreciated,

Cheers,

Will

The first error message in the log is:
INFO 19:59:32,028 FunctionEdge - Output written to /lustre/top/department/me/project/genotyping/CNVs/test3/logs/SVDiscovery-1.out INFO 19:59:32,132 DrmaaJobRunner - Submitted job id: 6979402 INFO 19:59:32,133 QGraph - 1 Pend, 1 Run, 0 Fail, 0 Done ERROR 20:06:02,035 FunctionEdge - Error: 'java' '-Xmx4096m' '-XX:+UseParallelOldGC' '-XX:ParallelGCThreads=4' '-XX:GCTimeLimit=50' '-XX:GCHeapFreeLimit=10' '-Djava.io.tmpdir=/lustre/scratch/tmp' '-cp' '/cm/shared/apps/svtoolkit/2.0.1602/lib/SVToolkit.jar: ERROR 20:06:02,044 FunctionEdge - Contents of /lustre/top/department/me/project/genotyping/CNVs/test3/logs/SVDiscovery-1.out: INFO 19:59:39,508 HelpFormatter - --------------------------------------------------------------------------- INFO 19:59:39,512 HelpFormatter - The Genome Analysis Toolkit (GATK) v<unknown>, Compiled 2015/05/15 09:12:56

The second error message in the main log is:
INFO 20:05:29,637 03-Dec-2015 ReadCountDiskCache - Initialized read count disk cache with 1 file. INFO 20:05:29,640 03-Dec-2015 SVDiscovery - No hapmap snp genotype directory specified INFO 20:05:29,641 03-Dec-2015 SVDiscovery - No array intensity data specified INFO 20:05:30,141 03-Dec-2015 SVDiscovery - Clustering: Generating clusters for 3 read pairs. INFO 20:05:30,144 03-Dec-2015 SVDiscovery - Clustering: Generating clusters for 8 read pairs. INFO 20:05:30,144 03-Dec-2015 SVDiscovery - Processing cluster chr2L:16452-16555 chr2L:19389-19502 LR 3 ERROR: ClusterDepthModule.computePValue 883 -885 5448 -5462 Error: Exception processing cluster: Error running script /cm/shared/apps/svtoolkit/2.0.1602/R/discovery/compute_depth_pvalue.R: Error in chisq.test(m) : all entries of 'x' must be nonnegative and finite Calls: main -> cat -> sprintf -> depth.chisq.1 Execution halted Cluster: chr2L:16452-16555 chr2L:19389-19502 LR 3 INFO 20:05:33,906 03-Dec-2015 GATKRunReport - Uploaded run statistics report to AWS S3

My commands are:
. /etc/profile.d/modules.sh module load genomestrip/2.0 module load sge
#export SV_DIR=cd .. && pwd/
#SV_TMPDIR=./tmpdir
echo ${RIGHT_NOW} ".... WORKFLOW PREPARATION .... " runDir=test3 bam=lhm_test_bams.list disco_genos=test3/test3.discovery.vcf final_genos=test3/test3.genotypes.vcf ref_seq=../../reference_sequences/dmel/v6.0/dm6.fa mkdir -p ${runDir}/logs || exit 1 mkdir -p ${runDir}/metadata || exit 1 which java > /dev/null || exit 1 which Rscript > /dev/null || exit 1 which samtools > /dev/null || exit 1 export PATH=${SV_DIR}/bwa:${PATH} export LD_LIBRARY_PATH=${SV_DIR}/bwa:${LD_LIBRARY_PATH} mx="-Xmx4g" classpath="${SV_DIR}/lib/SVToolkit.jar:${SV_DIR}/lib/gatk/GenomeAnalysisTK.jar:${SV_DIR}/lib/gatk/Queue.jar" java -cp ${classpath} ${mx} -jar ${SV_DIR}/lib/SVToolkit.jar

echo ${RIGHT_NOW} " .... PREPROCESSING DELETIONS .... " java -cp ${classpath} ${mx} org.broadinstitute.gatk.queue.QCommandLine \ -S ${SV_DIR}/qscript/SVPreprocess.q \ -S ${SV_DIR}/qscript/SVQScript.q \ -gatk ${SV_DIR}/lib/gatk/GenomeAnalysisTK.jar \ -jobRunner Drmaa \ -gatkJobRunner Drmaa \ -jobNative '-V -pe openmp 20 -q bioinf.q' \ -cp ${classpath} \ -configFile conf/genstrip_test3_parameters.txt \ -tempDir ${SV_TMPDIR} \ -R ${ref_seq} \ -runDirectory ${runDir} \ -md ${runDir}/metadata \ -disableGATKTraversal \ -ploidyMapFile conf/ploidy_dm6.map \ -genderMapFile data/test3_gender.map \ -useMultiStep \ -computeGCProfiles true \ -reduceInsertSizeDistributions true \ -computeReadCounts true \ -jobLogDir ${runDir}/logs \ -I ${bam} \ -run || exit 1

echo ${RIGHT_NOW} " .... DISCOVERING DELETIONS .... " java -cp ${classpath} ${mx} org.broadinstitute.gatk.queue.QCommandLine \ -S ${SV_DIR}/qscript/SVDiscovery.q \ -S ${SV_DIR}/qscript/SVQScript.q \ -gatk ${SV_DIR}/lib/gatk/GenomeAnalysisTK.jar \ -jobRunner Drmaa \ -gatkJobRunner Drmaa \ -jobNative '-V -pe openmp 20 -q bioinf.q' \ -cp ${classpath} \ -configFile conf/genstrip_test3_parameters.txt \ -tempDir ${SV_TMPDIR} \ -R ${ref_seq} \ -runDirectory ${runDir} \ -md ${runDir}/metadata \ -disableGATKTraversal \ -jobLogDir ${runDir}/logs \ -ploidyMapFile conf/ploidy_dm6.map \ -genderMapFile data/test3_gender.map \ -L chr2L \ -minimumSize 100 \ -maximumSize 20000 \ -I ${bam} \ -O ${disco_genos} \ -run || exit 1

echo ${RIGHT_NOW} " .... GENOTYPING DELETIONS .... " java -cp ${classpath} ${mx} org.broadinstitute.gatk.queue.QCommandLine \ -S ${SV_DIR}/qscript/SVGenotyper.q \ -S ${SV_DIR}/qscript/SVQScript.q \ -gatk ${SV_DIR}/lib/gatk/GenomeAnalysisTK.jar \ -jobRunner Drmaa \ -gatkJobRunner Drmaa \ -jobNative '-V -pe openmp 20 -q bioinf.q' \ -cp ${classpath} \ -configFile conf/genstrip_test3_parameters.txt \ -tempDir ${SV_TMPDIR} \ -R ${ref_seq} \ -runDirectory ${runDir} \ -md ${runDir}/metadata \ -disableGATKTraversal \ -jobLogDir ${runDir}/logs \ -ploidyMapFile conf/ploidy_dm6.map \ -genderMapFile data/test3_gender.map \ -parallelJobs 10 \ -L 2L \ -I ${bam} \ -vcf ${disco_genos} \ -O ${final_genos} \ -run || exit 1

Post edited by Will_Gilks on

Best Answers

Answers

  • bhandsakerbhandsaker Member, Broadie, Moderator admin
    edited December 2015

    There is probably something wrong with your reference metadata.

    If you can post or email me the following, that would be helpful:
    metadata/genome_sizes.txt
    metadata/depth.dat
    conf/ploidy_dm6.map
    dm6.fa.fai

    Also, did you create any other reference metadata files that are stored parallel to dm6.fa ?
    For example, dm6.svmask.fasta or dm6.gcmask.fasta ?
    Note that currently the automatic detection of the svmask and gcmask files requires the reference extension to be .fasta (not .fa).

  • Will_GilksWill_Gilks University of Sussex, UKMember ✭✭

    Thanks for getting back to me.

    I've attached the suggested metadata.

    There are other files stored in the same directory as dm6.fa. I'll delete these out as they're not for GS, and let you know how things go.

  • Will_GilksWill_Gilks University of Sussex, UKMember ✭✭

    Deleted other files in the reference assembly directory but still similar error messages. Starts with an R error in preprocessing. Not sure if this is critical though:
    INFO 17:00:47,638 QJobsReporter - Plotting JobLogging GATKReport to file /lustre/scratch/bioenv/wg39/LHm_analysis/genotyping/CNVs/SVPreprocess.jobre WARN 17:00:49,613 RScriptExecutor - RScript exited with 1. Run with -l DEBUG for more info. INFO 17:00:49,623 QCommandLine - Script completed successfully with 7615 total jobs

    Then in deletion discovery:
    ERROR 17:08:06,580 FunctionEdge - Error: 'java' '-Xmx4096m' '-XX:+UseParallelOldGC' '-XX:ParallelGCThreads=4' '-XX:GCTimeLimit=50' '-XX:GCHeapF ERROR 17:08:06,587 FunctionEdge - Contents of /lustre/scratch/bioenv/wg39/LHm_analysis/genotyping/CNVs/test3/logs/SVDiscovery-1.out:

    ERROR: ClusterDepthModule.computePValue 883 -885 5448 -5462 Error: Exception processing cluster: Error running script /cm/shared/apps/svtoolkit/2.0.1602/R/discovery/compute_depth_pvalue.R: Error in chisq.test( all entries of 'x' must be nonnegative and finite Calls: main -> cat -> sprintf -> depth.chisq.1

  • Will_GilksWill_Gilks University of Sussex, UKMember ✭✭

    The other thing I was wondering about was Samtools. I read somewhere that GS requires it but I don't see it specified in any command line. I'm currently using v1.2.

  • bhandsakerbhandsaker Member, Broadie, Moderator admin

    In rare cases, GS uses samtools just for samtools index to index bam files, so pretty much any version should work.

  • Will_GilksWill_Gilks University of Sussex, UKMember ✭✭
    edited December 2015

    Thanks @bhandsaker Have corrected ploidy file and made rdmask. To make the svmask do you still recommend computing the mask for each chromosome separately, and then concatenate? (Each of the four major autosomes are about 30MB.)

    http://gatkforums.broadinstitute.org/discussion/1499/computegenomemask

    Post edited by Will_Gilks on
  • Will_GilksWill_Gilks University of Sussex, UKMember ✭✭

    Hallelujah. Made ploidy file as you advised. svmask required an indexed fasta but easy to do.

    Below is some code for generating the svmask on a linux SGE without the hashed headers in case anyone else needs it.

    module load sge module load genomestrip/2.0 which java > /dev/null || exit 1 which Rscript > /dev/null || exit 1 which samtools > /dev/null || exit 1 which bwa > /dev/null || exit 1 export PATH=${SV_DIR}/bwa:${PATH} export LD_LIBRARY_PATH=${SV_DIR}/bwa:${LD_LIBRARY_PATH} ${SV_DIR}/bwa/bwa index -a bwtsw dm6.fa mx="-Xmx4g" classpath="${SV_DIR}/lib/SVToolkit.jar:${SV_DIR}/lib/gatk/GenomeAnalysisTK.jar:${SV_DIR}/lib/gatk/Queue.jar" java -cp ${classpath} ${mx} -jar ${SV_DIR}/lib/SVToolkit.jar java -cp ${classpath} ${mx} \ org.broadinstitute.sv.apps.ComputeGenomeMask \ -R dm6.fa \ -O dm6.svmask.fasta \ -readLength 100

    ${SV_DIR}/bwa/bwa index -a bwtsw dm6.svmask.fasta module load samtools/1.0 samtools faidx dm6.svmask.fasta module unload samtools/1.0 module unload genomestrip/2.0 module unload sge

Sign In or Register to comment.