This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Variant calling using GATK::ERROR

Dear All,

I am working to setup pipelines for variants calling. I used SOAPnuke, HISAT, for pre-process. Finally, I am using GATK to call snps/variants. I have tried a lot to figure out my error, please advice........

Here is my process:

  1. Removing adaptors and cleaning other parts using SOAPnuke

    ../SOAPnuke filter -l 10 -q 0.5 -n 0.05 -1 trt-unassembled.fastq -C clean_trt.fastq -o result -T 8

  2. Index and alignment with hisat2-2.2.0

    ./hisat2-build ..sequence.fastq ../PEDV

    ./hisat2 -x ../PEDV -U ../try.fastq -S ../trt.sam

  3. SAM to BAM

    samtools view -Sb ../trt.sam >../trt.bam

    samtools sort ../trt.bam >../sorted-trt.bam

    samtools index ../sorted-trt.bam

  4. making .dict and .fai

    java -jar picard.jar CreateSequenceDictionary R= ../sequence.fasta O= ../sequence.dict

    samtools faidx ../sequence.fasta

  5. SNPs Calling using GATK

    ./gatk HaplotypeCaller -R ../sequence.fasta -I ../sorted-trt.bam.bai -O ../variants-trt.vcf

-------ERROR: A USER ERROR has occurred: Input files reference and reads have incompatible contigs: No overlapping contigs found. reference contigs = [LGE3336] reads contigs = []

snp rna-seq sequencing
• 2,422 views • link
written by Kumar Kumar

Here is the .dict file content:

@HD VN:1.5
@SQ SN:KF267450.1   LN:28038    M5:96f71e94736a55d625aba2a04a7d0e11 UR:file:/home/kumarm/allan-work/bowtie2-trt/sequence.fasta

.fai file content

KF267450.1  28038   12  70  71
• 0 views • link
updated by GenoMax GenoMax • written by Kumar Kumar

Please use ADD COMMENT/ADD REPLY when responding to existing posts to keep threads logically organized. SUBMIT ANSWER is for new answers to original question.

• 0 views • link
written by GenoMax GenoMax

2 answers

./gatk HaplotypeCaller -R ../sequence.fasta -I ../sorted-trt.bam.bai -O ../variants-trt.vcf

is not sorted-trt.bam? the *bai file is the index

• 0 views • link
written by JC JC

try -I ../sorted-trt.bam instead of -I ../sorted-trt.bam.bai

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Still I am not getting any vcf output file.

command:

./gatk HaplotypeCaller -R ../sequence.fasta -I ../sorted-try.bam -O ../variants-trt.vcf

MESSAGE:

Using GATK jar /home/kumarm/allan-work/gatk-4.1.7.0/gatk-package-4.1.7.0-local.jar
Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/kumarm/allan-work/gatk-4.1.7.0/gatk-package-4.1.7.0-local.jar HaplotypeCaller -R /home/kumarm/allan-work/hisat-trt/sequence.fasta -I /home/kumarm/allan-work/hisat-trt/sorted-try.bam -O /home/kumarm/allan-work/hisat-trt/variants-trt.vcf
10:53:35.225 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/home/kumarm/allan-work/gatk-4.1.7.0/gatk-package-4.1.7.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
May 12, 2020 10:53:35 AM shaded.cloud_nio.com.google.auth.oauth2.ComputeEngineCredentials runningOnComputeEngine
INFO: Failed to detect whether we are running on Google Compute Engine.
10:53:35.425 INFO  HaplotypeCaller - ------------------------------------------------------------
10:53:35.427 INFO  HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.1.7.0
10:53:35.427 INFO  HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
10:53:35.427 INFO  HaplotypeCaller - Executing as kumarm@microbiome on Linux v4.15.0-99-generic amd64
10:53:35.427 INFO  HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v11.0.1+13-LTS
10:53:35.428 INFO  HaplotypeCaller - Start Date/Time: May 12, 2020 at 10:53:35 AM CDT
10:53:35.428 INFO  HaplotypeCaller - ------------------------------------------------------------
10:53:35.428 INFO  HaplotypeCaller - ------------------------------------------------------------
10:53:35.431 INFO  HaplotypeCaller - HTSJDK Version: 2.21.2
10:53:35.431 INFO  HaplotypeCaller - Picard Version: 2.21.9
10:53:35.431 INFO  HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
10:53:35.431 INFO  HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
10:53:35.431 INFO  HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
10:53:35.432 INFO  HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
10:53:35.432 INFO  HaplotypeCaller - Deflater: IntelDeflater
10:53:35.432 INFO  HaplotypeCaller - Inflater: IntelInflater
10:53:35.432 INFO  HaplotypeCaller - GCS max retries/reopens: 20
10:53:35.432 INFO  HaplotypeCaller - Requester pays: disabled
10:53:35.433 INFO  HaplotypeCaller - Initializing engine
10:53:35.621 INFO  HaplotypeCaller - Done initializing engine
10:53:35.629 INFO  HaplotypeCallerEngine - Disabling physical phasing, which is supported only for reference-model confidence output
10:53:35.638 INFO  HaplotypeCaller - Shutting down engine
[May 12, 2020 at 10:53:35 AM CDT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.01 minutes.
Runtime.totalMemory()=2113929216
java.lang.IllegalArgumentException: samples cannot be empty
        at org.broadinstitute.hellbender.utils.Utils.validateArg(Utils.java:727)
        at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.ReferenceConfidenceModel.<init>(ReferenceConfidenceModel.java:119)
        at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCallerEngine.initialize(HaplotypeCallerEngine.java:223)
        at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCallerEngine.<init>(HaplotypeCallerEngine.java:166)
        at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller.onTraversalStart(HaplotypeCaller.java:196)
        at org.broadinstitute.hellbender.engine.GATKTool.doWork(GATKTool.java:1046)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:139)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:191)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:210)
        at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:163)
        at org.broadinstitute.hellbender.Main.mainEntry(Main.java:206)
        at org.broadinstitute.hellbender.Main.main(Main.java:292)
• 0 views • link
updated by GenoMax GenoMax • written by Kumar Kumar

I am using HISAT2 for alignment (.sam) and converting .sam to .bam with following command...

Index and alignment with hisat2-2.2.0

./hisat2-build ..sequence.fastq ../PEDV

./hisat2 -x ../PEDV -U ../try.fastq -S ../trt.sam
• 0 views • link
updated by GenoMax GenoMax • written by Kumar Kumar

Please let me know how to prepare the command for adding read group.

Here is the .dict file content...

@HD VN:1.5
@SQ SN:KF267450.1   LN:28038    M5:96f71e94736a55d625aba2a04a7d0e11 UR:file:/home/kumarm/allan-work/bowtie2-trt/sequence.fasta
• 0 views • link
updated by GenoMax GenoMax • written by Kumar Kumar

I took information from .dict to prepare command to add read group for analysis using HISAT2. command:

./hisat2 --rg-id "@RG\tID:KF267450.1\tLN:28038\tM5:96f71e94736a55d625aba2a04a7d0e11" -f -x ../hisat-trt/PEDV -U ../hisat-trt/clean_trt.fasta -S ../hisat-trt/trt.sam

Is this correct command?

-- Here is the log process of alignment at HISAT2. 

15065158 reads; of these

15065158 (100.00%) were unpaired; of these

15065081 (100.00%) aligned 0 times

77 (0.00%) aligned exactly 1 time

0 (0.00%) aligned >1 times

0.00% overall alignment rate

• 0 views • link
written by Kumar Kumar

I have figured out the issue. I added read group like this.

./hisat2 --rg ID:KF267450.1 --rg SM:501 --rg PL:ILLUMINA --rg LB:lib-501 --rg PU:C3MF6ACXX.1.NoIndex

• 0 views • link
written by Kumar Kumar

I tried after adding read group. However I am still not getting vcf file.

I tried following commands:

RG="@RG\tID:XX.L001\tSM:1\tPL:ILLUMINA\tLB:lib501\tPU:XX.1.NoIndex"

./hisat2 --rg-id --rg-id $RG -f -x ../PEDV -U ../try.fasta -S ../trt.sam

• 0 views • link
written by Kumar Kumar

Log in to answer this question.

More posts like this