4
$\begingroup$

I am trying to perform variant calling on a BAM file generated through STAR version STAR_2.6.0b for wheat genome using GATK haplotypecaller as follows:

gatk HaplotypeCaller -I sorted.bam -R wheat.fa -O st.vc
Using GATK jar /home/geslab/bin/gatk-package-4.0.4.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/geslab/bin/gatk-package-4.0.4.0-local.jar HaplotypeCaller -I sorted.bam -R wheat.fa -O st.vc                                                                                                                                                                                       
02:34:13.450 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/home/geslab/bin/gatk-package-4.0.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so                                                 
02:34:13.600 INFO  HaplotypeCaller - ------------------------------------------------------------                                                                                                                                
02:34:13.600 INFO  HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.0.4.0                                                                                                                                                 
02:34:13.600 INFO  HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/                                                                                                               
02:34:13.601 INFO  HaplotypeCaller - Initializing engine                                                                                                                                                                         
02:34:40.721 INFO  HaplotypeCaller - Done initializing engine                                                                                                                                                                    
02:34:48.175 INFO  HaplotypeCallerEngine - Disabling physical phasing, which is supported only for reference-model confidence output                                                                                             
02:34:48.598 INFO  HaplotypeCaller - Shutting down engine                                                                                                                                                                        
[June 7, 2018 2:34:48 AM PKT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.59 minutes.                                                                                      
Runtime.totalMemory()=10305929216                                                                                                                                                                                                
java.lang.IllegalArgumentException: samples cannot be empty

When I searched for this error the most proposed solution was to validate BAM file as suggested here. The ValidateSamFile gave following output for the bam file:

gatk ValidateSamFile -I sorted.bam -M SUMMARY
Using GATK jar /home/geslab/bin/gatk-package-4.0.4.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/geslab/bin/gatk-package-4.0.4.0-local.jar ValidateSamFile -I sorted.bam -M SUMMARY
02:53:27.498 INFO  NativeLibraryLoader - Loading libgkl_compression.so   from jar:file:/home/geslab/bin/gatk-package-4.0.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
[Thu Jun 07 02:53:27 PKT 2018] ValidateSamFile  --INPUT sorted.bam --MODE   SUMMARY  --MAX_OUTPUT 100 --IGNORE_WARNINGS false --VALIDATE_INDEX true --INDEX_VALIDATION_STRINGENCY EXHAUSTIVE --IS_BISULFITE_SEQUENCED false --MAX_OPEN_TEMP_FILES 8000 --SKIP_MATE_VALIDATION false --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
INFO    2018-06-07 02:54:08     SamFileValidator        Validated Read    10,000,000 records.  Elapsed time: 00:00:36s.  Time for last 10,000,000:   36s.  Last read position: TGACv1_scaffold_080127_1DS:28,522
INFO    2018-06-07 02:54:44     SamFileValidator        Validated Read    20,000,000 records.  Elapsed time: 00:01:12s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_602746_7DL:127,507
INFO    2018-06-07 02:55:20     SamFileValidator        Validated Read    30,000,000 records.  Elapsed time: 00:01:48s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_641946_U:17,161
INFO    2018-06-07 02:55:55     SamFileValidator        Validated Read    40,000,000 records.  Elapsed time: 00:02:23s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_471755_6AL:43,097
INFO    2018-06-07 02:56:31     SamFileValidator        Validated Read    50,000,000 records.  Elapsed time: 00:02:59s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_434075_5DL:29,039
INFO    2018-06-07 02:57:07     SamFileValidator        Validated Read    60,000,000 records.  Elapsed time: 00:03:34s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_062491_1DL:40,466
INFO    2018-06-07 02:57:42     SamFileValidator        Validated Read    70,000,000 records.  Elapsed time: 00:04:10s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_558158_7AL:13,072
INFO    2018-06-07 02:58:18     SamFileValidator        Validated Read    80,000,000 records.  Elapsed time: 00:04:46s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_435811_5DL:29,008
INFO    2018-06-07 02:58:53     SamFileValidator        Validated Read    90,000,000 records.  Elapsed time: 00:05:21s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_572030_7AS:9,114


## HISTOGRAM    java.lang.String
Error Type      Count
ERROR:MATE_NOT_FOUND    107771
ERROR:MISSING_READ_GROUP        1
WARNING:MISSING_TAG_NM  96856336
WARNING:RECORD_MISSING_READ_GROUP       96856336

In order to remove the errors shown by ValidateSamFile I used FixMateInformation as suggested here. as shown below(sharing part of file):

gatk FixMateInformation -I sorted.bam -O mate.bam -MC
Using GATK jar /home/geslab/bin/gatk-package-4.0.4.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/geslab/bin/gatk-package-4.0.4.0-local.jar FixMateInformation -I sorted.bam -O mate.bam -MC
06:13:13.934 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/home/geslab/bin/gatk-package-4.0.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
[Thu Jun 07 06:13:19 PKT 2018] FixMateInformation  --INPUT sorted.bam --OUTPUT mate.bam --ADD_MATE_CIGAR true  --ASSUME_SORTED false --IGNORE_MISSING_MATES true --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
INFO    2018-06-07 06:13:24     FixMateInformation      Sorting input into queryname order.
INFO    2018-06-07 06:23:26     FixMateInformation      Sorting by queryname complete.
INFO    2018-06-07 06:23:26     FixMateInformation      Output will be sorted by coordinate
INFO    2018-06-07 06:23:26     FixMateInformation      Traversing query name sorted records and fixing up mate pair information.
INFO    2018-06-07 06:23:33     FixMateInformation      Processed      1,000,000 records.  Elapsed time: 00:00:06s.  Time for last 1,000,000:    6s.    Last read position: TGACv1_scaffold_361122_4DS:70,952
INFO    2018-06-07 06:23:42     FixMateInformation      Processed     2,000,000 records.  Elapsed time: 00:00:15s.  Time for last 1,000,000:    8s.    Last read position: TGACv1_scaffold_044571_1BL:711
INFO    2018-06-07 06:23:50     FixMateInformation      Processed     3,000,000 records.  Elapsed time: 00:00:23s.  Time for last 1,000,000:    8s.   Last read position: TGACv1_scaffold_433135_5DL:36,647
 .
 .
 .
 .
 [Thu Jun 07 06:40:50 PKT 2018] picard.sam.FixMateInformation done.    Elapsed time: 27.62 minutes.
 Runtime.totalMemory()=8004304896
 Tool returned:
 0

After performing FixMateInformation I again validated the BAM file using ValidateSamFile but the errors are still the same as shown above which shows Mate not found, which means error is not fixed.

How these errors can be fixed? So that I can perform variant calling using BAM file?

Improve this question
$\endgroup$

3 Answers 3

Reset to default
1
$\begingroup$

The primary error isn't the missing mate, but likely the missing read groups. Use Picard AddOrReplaceReadGroups to ensure your BAM files have read groups and then run those through the haplotype caller.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Should I run AddOrReplaceReadGroups with default parameters as discussed here or there is something that needs to be extracted and provided as parameters for AddOrReplaceReadGroups? $\endgroup$
    – Ammar Sabir Cheema
    Jun 8, 2018 at 10:11
  • $\begingroup$ There are no defaults for the read group information. You can just make up values to put there (usually I just use the sample name). Note that you only need the required fields. $\endgroup$
    – Devon Ryan
    Jun 8, 2018 at 10:23
0
$\begingroup$

You could add read groups: 

['STAR', '--twopassMode', 'Basic', '--genomeDir', '/data2/Fshare/FastaAndIndex/IWGSC_v1.0_STAR/',
         '--runThreadN', '10', '--limitSjdbInsertNsj', '5000000', '--outSAMtype', 'BAM', 'SortedByCoordinate', '--twopass1readsN', '-1',
         '--sjdbOverhang', '100', '--readFilesIn', sra1 + '_1.fastq', sra1 + '_2.fastq',
         '--outSAMattrRGline', 'ID:' + rg1, 'SM:' + rg1, 'PL:ILLUMINA'], shell=False)
Improve this answer
$\endgroup$
3
  • $\begingroup$ In other words, "add read groups". $\endgroup$
    – Devon Ryan
    Jun 8, 2018 at 6:47
  • $\begingroup$ @user3018491 on which version of STAR you have used the parameter SortedByCoordinate? As discussed here I tried to use this on version 2.6.0b but it failed. $\endgroup$
    – Ammar Sabir Cheema
    Jun 8, 2018 at 10:16
  • $\begingroup$ STAR Usage: STAR [options]... --genomeDir REFERENCE --readFilesIn R1.fq R2.fq Spliced Transcripts Alignment to a Reference (c) Alexander Dobin, 2009-2015 ### versions versionSTAR 020201 $\endgroup$
    – masw
    Jun 8, 2018 at 13:51
0
$\begingroup$

The problem is that you are missing read groups. You can add them directly to bam files you already have using picard tools. Something like

java -jar ./software/picard.jar AddOrReplaceReadGroups \
  I=input.bam \
  O=output.bam \
  RGLB=lib1 \
  RGPL=illumina \
  RGPU=unit1 \
  RGSM=hello
Improve this answer
$\endgroup$
1
  • $\begingroup$ Hello 罗素幸福之路, it would be nice to explain what the code you have posted means. I tried to add at least one sentence, feel free to edit it and improve it further. $\endgroup$
    – Kamil S Jaron
    May 29, 2019 at 13:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.