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.

GenomeSTRiP query: Empty seq_1.merged.sites.vcf.gz file in Stage 3 of CNVDiscovery

Dear all,

While running the CNVDiscovery pipeline, we have found that the seq_1.merged.sites.vcf.gz file which is outputted by Stage 3 is empty. The CNVDiscovery pipeline does continue to run, however, we would like to know whether this is normal or whether it may cause a problem later on.

Below is the CNVDiscoveryStage3-1.out log for chromosome 1: we noticed that 0 output sites were generated
(the same is true for all chromosomes)

```
INFO 11:40:24,538 HelpFormatter - --------------------------------------------------------------------
INFO 11:40:24,549 HelpFormatter - Program Name: org.broadinstitute.sv.discovery.MergeDepthScannerSites
INFO 11:40:24,553 HelpFormatter - Program Args: -O /home/lcottino/attempt5/run1/cnv_stage3/seq_1/seq_1.merged.sites.vcf.gz -vcf run1/cnv_stage2/seq_1/seq_1.genotypes.vcf.gz -R /home/lcottino/attempt3/Hom
o_sapiens_assembly19/Homo_sapiens_assembly19.fasta -genomeMaskFile /home/lcottino/attempt3/Homo_sapiens_assembly19/Homo_sapiens_assembly19.svmask.fasta -genomeMaskFile /home/lcottino/attempt3/Homo_sapiens
_assembly19/Homo_sapiens_assembly19.lcmask.fasta -reportFile run1/cnv_stage3/seq_1/MergeScannerSites.report.dat -outputVariants ALL -duplicateScoreThresholdMax 0.0
INFO 11:40:24,557 HelpFormatter - Executing as [email protected] on Linux 3.10.0-957.1.3.el7.x86_64 amd64; OpenJDK 64-Bit Server VM 1.8.0_191-b12.
INFO 11:40:24,557 HelpFormatter - Date/Time: 2019/02/04 11:40:24
INFO 11:40:24,557 HelpFormatter - --------------------------------------------------------------------
INFO 11:40:24,557 HelpFormatter - --------------------------------------------------------------------
INFO 11:40:24,558 MergeDepthScannerSites - Opening reference sequence ...
INFO 11:40:24,569 MergeDepthScannerSites - Opened reference sequence.
INFO 11:40:24,570 MergeDepthScannerSites - Opening genome mask ...
INFO 11:40:24,572 MergeDepthScannerSites - Opened genome mask.
INFO 11:40:30,789 MergeDepthScannerSites - Processed 427869 input sites.
INFO 11:40:30,791 MergeDepthScannerSites - Generated 0 output sites.
INFO 11:40:30,791 CommandLineProgram - Program completed.
------------------------------------------------------------------------------------------
Done. There were no warn messages.
------------------------------------------------------------------------------------------
```

Many thanks for any help.

Answers

  • bhandsakerbhandsaker Member, Broadie, Moderator admin

    This definitely indicates a problem. Probably the problem is happening earlier, in stage2.
    I would start by reviewing the command line carefully to see if there is some problem (you could post it).
    I would also suggest looking in cnv_stage2/seq_1/eval. For example at CopyNumberClass.report.dat, which is a text file, or some of the other report files there, which may give some hints, and possibly the log files from stage2.

  • @bhandsaker

    Thank you for your response.

    Here is the code that we ran:
    ```
    BAMFILE=/home/lcottino/attempt5/sahgp.list

    export PATH=$PATH:/home/lcottino/samtools-1.9
    export PATH=$PATH:/home/lcottino/attempt1/tabix-0.2.6
    export PATH=$PATH:/home/lcottino/attempt1/samtools-1.9/htslib-1.9/tabix
    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/SVPreprocess.q \
    -S ${SV_DIR}/qscript/SVQScript.q \
    -cp ${classpath} \
    -gatk ${SV_DIR}/lib/gatk/GenomeAnalysisTK.jar \
    -configFile ${SV_DIR}/conf/genstrip_parameters.txt \
    -R /home/lcottino/attempt3/Homo_sapiens_assembly19/Homo_sapiens_assembly19.fasta \
    -I ${BAMFILE} \
    -md /home/lcottino/attempt5/output_metadata_directory_1 \
    -jobRunner Drmaa \
    -bamFilesAreDisjoint true \
    -jobLogDir logDir \
    -run

    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 /home/lcottino/attempt3/Homo_sapiens_assembly19/Homo_sapiens_assembly19.fasta \
    -I ${BAMFILE}\
    -md /home/lcottino/attempt5/output_metadata_directory_1 \
    -genderMapFile /home/lcottino/attempt5/gender.map \
    -runDirectory run1 \
    -jobLogDir run1 \
    -jobRunner Drmaa \
    -gatkJobRunner Drmaa \
    -tilingWindowSize 1000 \
    -tilingWindowOverlap 500 \
    -maximumReferenceGapLength 1000 \
    -boundaryPrecision 100 \
    -minimumRefinedLength 500 \
    -run
    ```

    The CopyNumberClass.report.dat (from Stage 2) appears to be empty: below are the first few lines
    ```
    ID CALLRATE CNMIN CNMAX CNALLELES NNONREF NVARIANT CNCATEGORY CNDIST
    CNV_1_13028_16230 0.000 NA NA 0 0 0 NA NA
    CNV_1_14941_18893 0.000 NA NA 0 0 0 NA NA
    CNV_1_16230_29415 0.000 NA NA 0 0 0 NA NA
    CNV_1_18893_40236 0.000 NA NA 0 0 0 NA NA
    CNV_1_29415_47299 0.000 NA NA 0 0 0 NA NA
    CNV_1_40236_51628 0.000 NA NA 0 0 0 NA NA
    CNV_1_47299_52709 0.000 NA NA 0 0 0 NA NA
    CNV_1_51628_54502 0.000 NA NA 0 0 0 NA NA
    CNV_1_52709_55412 0.000 NA NA 0 0 0 NA NA
    ```
    The VariantsPerSample.report.dat, NonVariant.report.dat, SelectedVariants.list, ClusterSeparation.report.dat, GenotypeLikelihoodStats.report.dat and NonVariant.report.dat from Stage 2 all appear to be empty as well.



    We had to reheader the BAM files as they did not have a library tag. Below is one of the BAM headers - could we possibly have made an error here?
    ```
    @HD VN:1.0 GO:none SO:coordinate
    @SQ SN:chrM LN:16571 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:d2ed829b8a1628d16cbeee88e88e39eb
    @SQ SN:chr1 LN:249250621 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:1b22b98cdeb4a9304cb5d48026a85128
    @SQ SN:chr2 LN:243199373 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:a0d9851da00400dec1098a9255ac712e
    @SQ SN:chr3 LN:198022430 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:641e4338fa8d52a5b781bd2a2c08d3c3
    @SQ SN:chr4 LN:191154276 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:23dccd106897542ad87d2765d28a19a1
    @SQ SN:chr5 LN:180915260 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:0740173db9ffd264d728f32784845cd7
    @SQ SN:chr6 LN:171115067 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:1d3a93a248d92a729ee764823acbbc6b
    @SQ SN:chr7 LN:159138663 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:618366e953d6aaad97dbe4777c29375e
    @SQ SN:chr8 LN:146364022 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:96f514a9929e410c6651697bded59aec
    @SQ SN:chr9 LN:141213431 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:3e273117f15e0a400f01055d9f393768
    @SQ SN:chr10 LN:135534747 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:988c28e000e84c26d552359af1ea2e1d
    @SQ SN:chr11 LN:135006516 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:98c59049a2df285c76ffb1c6db8f8b96
    @SQ SN:chr12 LN:133851895 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:51851ac0e1a115847ad36449b0015864
    @SQ SN:chr13 LN:115169878 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:283f8d7892baa81b510a015719ca7b0b
    @SQ SN:chr14 LN:107349540 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:98f3cae32b2a2e9524bc19813927542e
    @SQ SN:chr15 LN:102531392 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:e5645a794a8238215b2cd77acb95a078
    @SQ SN:chr16 LN:90354753 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:fc9b1a7b42b97a864f56b348b06095e6
    @SQ SN:chr17 LN:81195210 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:351f64d4f4f9ddd45b35336ad97aa6de
    @SQ SN:chr18 LN:78077248 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:b15d4b2d29dde9d3e4f93d1d0f2cbc9c
    @SQ SN:chr19 LN:59128983 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:1aacd71f30db8e561810913e0b72636d
    @SQ SN:chr20 LN:63025520 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:0dec9660ec1efaaf33281c0d5ea2560f
    @SQ SN:chr21 LN:48129895 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:2979a6085bfe28e3ad6f552f361ed74d
    @SQ SN:chr22 LN:51304566 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:a718acaa6135fdca8357d5bfe94211dd
    @SQ SN:chrX LN:155270560 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:7e0e2e580297b7764e31dbc80c2540dd
    @SQ SN:chrY LN:59373566 UR:/illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa M5:1fa3474750af0948bdf97d5a0ee52e51
    @RG ID:0 LB:LP6005857-DNA_A01 PL:ILLUMINA PU:C3Y3VACXX_0:5:none SM:LP6005857-DNA_A01
    @RG ID:1 LB:LP6005857-DNA_A01 PL:ILLUMINA PU:C3Y3VACXX_0:6:none SM:LP6005857-DNA_A01
    @RG ID:2 LB:LP6005857-DNA_A01 PL:ILLUMINA PU:C3Y3VACXX_0:7:none SM:LP6005857-DNA_A01
    @PG ID:iSAAC PN:iSAAC VN:iSAAC-01.13.04.29 CL:/illumina/development/Isis/2.3.23.1/Workflows/ResequencingWorker/isaac-align --qscore-bin 1 --stop-at Bam --reference-genome /illumina/development/Isis/Genomes/Homo_sapiens/UCSC/hg19/Sequence/IsaacIndex/sorted-reference.xml --memory-limit 122 --output-directory /scratch/LP6005857-DNA_A01_WF40184/Alignment/IsaacAlignment --temp-directory /scratch/LP6005857-DNA_A01_WF40184/Alignment/IsaacTemp --base-quality-cutoff 15 --pf-only 1 --gap-scoring bwa --keep-duplicates 1 --keep-unaligned back --stats-image-format none --variable-fastq-read-length yes --ignore-missing-bcls 1 --ignore-missing-filters 1 --realign-gaps no --base-calls-directory /scratch/LP6005857-DNA_A01_WF40184/RunData/140404_SN439_0263_AC3Y3VACXX_sd-seq/Data/Intensities/BaseCalls --barcode-mismatches 1 --use-bases-mask Y100n,Y100n --tiles s_5,s_6,s_7 --base-calls-format bcl-gz -s /scratch/LP6005857-DNA_A01_WF40184/Alignment/tmp.0.0.dat
    @PG ID:GATK IndelRealigner VN:1.6-22-g3ec78bd CL:knownAlleles=[] targetIntervals=/scratch/LP6005857-DNA_A01_WF40184/Alignment/LP6005857-DNA-A01_S1_chrM.intervals LODThresholdForCleaning=5.0 consensusDeterminationModel=USE_READS entropyThreshold=0.15 maxReadsInMemory=150000 maxIsizeForMovement=3000 maxPositionalMoveAllowed=200 maxConsensuses=30 maxReadsForConsensuses=120 maxReadsForRealignment=20000 noOriginalAlignmentTags=false nWayOut=null generate_nWayOut_md5s=false check_early=false noPGTag=false keepPGTags=false indelsFileForDebugging=null statisticsFileForDebugging=null SNPsFileForDebugging=null
    @PG ID:Illumina_WGS_Services VN:v2.0.2
    @PG ID:samtools PN:samtools PP:iSAAC VN:1.9 CL:samtools reheader - LP6005857-DNA_A01.bam
    @PG ID:samtools.1 PN:samtools PP:GATK IndelRealigner VN:1.9 CL:samtools reheader - LP6005857-DNA_A01.bam
    @PG ID:samtools.2 PN:samtools PP:Illumina_WGS_Services VN:1.9 CL:samtools reheader - LP6005857-DNA_A01.bam
    ```
    Your help is greatly appreciated.
  • bhandsakerbhandsaker Member, Broadie, Moderator admin

    Ah, yes. You have to be careful with hg19/b37 because there were several different reference genomes used. (The community has gotten its act together around having a single reference for hg38.)
    I suspect the root of your problem is that your bam header lists the chromosomes as "chr1", "chr2", etc. whereas the reference you are using lists them as "1", "2", etc.
    You should first make sure which reference your reads were aligned against.
    Then try to use a Genome STRiP reference bundle that matches exactly.
    The available bundles are here: ftp://ftp.broadinstitute.org/pub/svtoolkit/reference_metadata_bundles/

  • bhandsakerbhandsaker Member, Broadie, Moderator admin

    You will need to rerun SVPreprocess, not just the CNV pipeline.

Sign In or Register to comment.