Problem running GenomeSTRiP CNVDiscovery on X chromosome
Hi there,
We have run the GenomeSTRiP CNVDiscovery pipeline on all of the autosomes as well as the Y chromosome of ~500 whole genome sequences successfully (we ran each chromosome separately). This is the script that was used: Just to note that we excluded some nodes because the bam files can only be accessed from certain nodes.
#!/bin/bash
#SBATCH --exclude=n01,n03,n04,n05,n06,n07,n09,n11,n12,n13,n14,n15,n16,n17,n18,n19,n20,n21,n22,n23,n25,n26,n27,n28
#SBATCH --mem=1gb
BAMFILE=/spaces/lcottino/H3A_CNVs/NEW_input.list
REF=/spaces/lcottino/H3A_CNVs//Homo_sapiens_assembly19/Homo_sapiens_assembly19.fasta
MD=/spaces/lcottino/H3A_CNVs/GSRun1/output_metadata_directory
GENDER=/spaces/lcottino/H3A_CNVs/gender2.map
CHROM=$1
export PATH=$PATH:/home/lcottino/samtools-1.9
export PATH=$PATH:/home/lcottino/attempt1/tabix-0.2.6
export LD_LIBRARY_PATH=/opt/exp_soft/slurm/lib:/opt/exp_soft/slurm/lib64
export LD_RUN_PATH=/opt/exp_soft/slurm/lib
export SV_DIR=/opt/exp_soft/bioinf/svtoolkit
classpath="${SV_DIR}/lib/SVToolkit.jar:${SV_DIR}/lib/gatk/GenomeAnalysisTK.jar:${SV_DIR}/lib/gatk/Queue.jar"
java -Xmx2g -cp ${classpath} \
org.broadinstitute.gatk.queue.QCommandLine \
-S ${SV_DIR}/qscript/discovery/cnv/CNVDiscoveryPipeline.q \
-S ${SV_DIR}/qscript/SVQScript.q \
-cp ${classpath} \
-gatk ${SV_DIR}/lib/gatk/GenomeAnalysisTK.jar \
-configFile ${SV_DIR}/conf/genstrip_parameters.txt \
-R ${REF} \
-I ${BAMFILE} \
-md ${MD} \
-genderMapFile ${GENDER} \
-runDirectory Run_Chr${CHROM} \
-jobLogDir Run_Chr${CHROM} \
-jobRunner Drmaa \
-gatkJobRunner Drmaa \
-jobNative "-x n01,n03,n04,n05,n06,n07,n09,n11,n12,n13,n14,n15,n16,n17,n18,n19,n20,n21,n22,n23,n25,n26,n27,n28"\
-intervalList ${CHROM} \
-tilingWindowSize 1000 \
-tilingWindowOverlap 500 \
-maximumReferenceGapLength 1000 \
-boundaryPrecision 100 \
-minimumRefinedLength 500 \
-maxConcurrentRun 20 \
-run
However, when we try to run the X chromosome it keeps failing at stage 7. When we look at some of the error files, there seems to be no actual error (see example below). They all end with “Done. There were no warn messages.” We have also tried limiting the number of jobs to 60, 40 and 20 and it still fails.
INFO 15:18:37,050 HelpFormatter - ------------------------------------------------------------------
INFO 15:18:37,054 HelpFormatter - Program Name: org.broadinstitute.sv.genotyping.RefineCNVBoundaries
INFO 15:18:37,061 HelpFormatter - Program Args: -I /spaces/lcottino/H3A_CNVs/GSRun2/Run_ChrX/cnv_stage6/seq_X/seq_X.merged_headers.bam -O /spaces/lcottino/H3A_CNVs/GSRun2/Run_ChrX/cnv_stage7/seq_X/P6063/
seq_X.merged.brig.vcf -R /spaces/lcottino/H3A_CNVs/Homo_sapiens_assembly19/Homo_sapiens_assembly19.fasta -md /spaces/lcottino/H3A_CNVs/GSRun1/output_metadata_directory -configFile /opt/exp_soft/bioinf/svt
oolkit/conf/genstrip_parameters.txt -P depth.readCountCacheIgnoreGenomeMask:true -genomeMaskFile /spaces/lcottino/H3A_CNVs/Homo_sapiens_assembly19/Homo_sapiens_assembly19.svmask.fasta -genomeMaskFile /spa
ces/lcottino/H3A_CNVs/Homo_sapiens_assembly19/Homo_sapiens_assembly19.lcmask.fasta -genderMapFile /spaces/lcottino/H3A_CNVs/gender2.map -ploidyMapFile /spaces/lcottino/H3A_CNVs/Homo_sapiens_assembly19/Hom
o_sapiens_assembly19.ploidymap.txt -vcf Run_ChrX/cnv_stage4/seq_X/seq_X.merged.genotypes.vcf.gz -site Run_ChrX/cnv_stage7/seq_X/P6063.sites.list -boundaryPrecision 100 -minimumRefinedLength 500 -maximumRe
ferenceGapLength 1000
INFO 15:18:37,066 HelpFormatter - Executing as lcottino@n10.core.wits.ac.za on Linux 3.10.0-957.27.2.el7.x86_64 amd64; OpenJDK 64-Bit Server VM 1.8.0_242-b08.
INFO 15:18:37,066 HelpFormatter - Date/Time: 2020/02/13 15:18:37
INFO 15:18:37,066 HelpFormatter - ------------------------------------------------------------------
INFO 15:18:37,067 HelpFormatter - ------------------------------------------------------------------
INFO 15:18:37,076 RefineCNVBoundaries - Opening reference sequence ...
INFO 15:18:37,077 RefineCNVBoundaries - Opened reference sequence.
INFO 15:18:37,090 RefineCNVBoundaries - Opening genome mask ...
INFO 15:18:37,092 RefineCNVBoundaries - Opened genome mask.
INFO 15:18:37,095 MetaData - Opening metadata ...
INFO 15:18:37,095 MetaData - Adding metadata location /spaces/lcottino/H3A_CNVs/GSRun1/output_metadata_directory ...
INFO 15:18:37,100 MetaData - Opened metadata.
INFO 15:18:37,112 RefineCNVBoundaries - Initializing input data set ...
INFO 15:18:37,479 RefineCNVBoundaries - Initialized data set: 1 file, 548 read groups, 548 samples.
INFO 15:18:37,630 ReadCountCache - Initializing read count cache with 1 file.
INFO 15:43:21,207 CommandLineProgram - Program completed.
------------------------------------------------------------------------------------------
Done. There were no warn messages.
------------------------------------------------------------------------------------------
We ran the X chromosome on a smaller set of genomes (n=20) and it ran to completion so there doesn’t seem to be a problem with the program itself but rather the volume.
Do you have any suggestions as to how to proceed?
Thank you in advance,
Laura
-
Thank you for your question! I am tagging Bob Handsaker in this thread and he should be able to help out with this. Thanks in advance Bob!
-
As you probably know, you can run Queue multiple times with the same command line and it will retry any jobs that previously failed. This allows you to work around transient problems (disk timeout, hung node, etc.).
When you rerun, does Queue retry the jobs? Does Queue report that some jobs are still failing? Do some jobs succeed when you retry?
-
We have resumed it and it take some time (over 12 hours) to fail again. When we look at the output files it seems as if it is the same jobs that are failing.
-
You'll have to provide more information for me to be able to help. You said it seemed as if the same jobs were failing. Presumably k jobs were retried - how many succeeded and how many failed? Were there error message in the log files from the failed jobs?
If you are getting no error messages in the log files and yet Queue thinks the jobs failed, that is probably some interaction with the job scheduler (e.g. maybe the job scheduler has lost the job status).
Something you may be able to use to help debug is the -jobWrapperScript path/to/script.sh option. This script is "wrapped" around each command that is run. So you could use a wrapper like this for example:
#!/bin/bash
echo "Wrapped command: ${@}"
"${@}" 2>&1
status=$?
if [ $status != 0 ]; then echo "Command failed: $status"; exit $status; fi
echo "Command completed successfully." -
If you can be sure that the commands succeeded, it is also possible to manually touch the .*.done files Queue uses for job tracking and then Queue will believe these jobs ran successfully on the next retry.
Please sign in to leave a comment.
5 comments