ASEReadCounter ouputs only header without error

Hi,

I have a problem similar to this:
HaplotypeCaller output header and one position recode without error

but with GATK4 ASEReadCounter. I have VCF file created using Samtools mpileup, compressed with bgzip, duplicate variants removed by bcftools norm and VCF file indexed by tabix. RNA-Seq BAM file was aligned to hg38 with STAR and sorted and indexed with Samtools. I also used AddOrReplaceReadGroups and ValidateSamFile for the BAM file.

Command:

module load gatk-env/4.0.2.1
gatk ASEReadCounter \--input sample_readgroups.bam \--variant sample.vcf.gz \--disable-tool-default-read-filters true \--reference hg38.fa \--output sample.ASEReadCounter.rtable \--java-options '-DGATK_STACKTRACE_ON_USER_EXCEPTION=true'

What happens:

----------------------------------------------------------------------------
  Setting up new environment, removing all currently loaded modules
----------------------------------------------------------------------------
  Loading new modules:
    r-env/3.2.5    rstudio    java/oracle/1.8    picard/2.13.2
----------------------------------------------------------------------------
  Setting up new environment, removing all currently loaded modules
----------------------------------------------------------------------------
  Loading new modules:
    gcc/4.9.3    intelmpi/5.1.1    mkl/11.3.0    r-app/3.2.5
r-app R-3.2.5 environment loaded

rstudio v0.99.1196 environment loaded

Start RStudio with rstudio

For interactive work consider using taito-shell.csc.fi

Loading application Oracle Java 1.8 
picard 2.13.2 environment loaded

Due to MODULEPATH changes the following have been reloaded:
  1) intelmpi/5.1.1  2) mkl/11.3.0

10:15:09.824 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/homeappl/appl_taito/bio/GATK4/gatk-4.0.2.1/gatk-package-4.0.2.1-local.jar!/com/intel/gkl/native/libgkl_compression.so
10:15:10.259 INFO  ASEReadCounter - ------------------------------------------------------------
10:15:10.260 INFO  ASEReadCounter - The Genome Analysis Toolkit (GATK) v4.0.2.1
10:15:10.260 INFO  ASEReadCounter - For support and documentation go to https://software.broadinstitute.org/gatk/
10:15:10.261 INFO  ASEReadCounter - Executing as [email protected] on Linux v2.6.32-642.15.1.el6.x86_64 amd64
10:15:10.261 INFO  ASEReadCounter - Java runtime: Java HotSpot(TM) 64-Bit Server VM v1.8.0_151-b12
10:15:10.261 INFO  ASEReadCounter - Start Date/Time: July 5, 2018 10:15:09 AM GMT
10:15:10.261 INFO  ASEReadCounter - ------------------------------------------------------------
10:15:10.262 INFO  ASEReadCounter - ------------------------------------------------------------
10:15:10.264 INFO  ASEReadCounter - HTSJDK Version: 2.14.3
10:15:10.264 INFO  ASEReadCounter - Picard Version: 2.17.2
10:15:10.264 INFO  ASEReadCounter - HTSJDK Defaults.COMPRESSION_LEVEL : 1
10:15:10.264 INFO  ASEReadCounter - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
10:15:10.265 INFO  ASEReadCounter - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
10:15:10.265 INFO  ASEReadCounter - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
10:15:10.265 INFO  ASEReadCounter - Deflater: IntelDeflater
10:15:10.265 INFO  ASEReadCounter - Inflater: IntelInflater
10:15:10.265 INFO  ASEReadCounter - GCS max retries/reopens: 20
10:15:10.265 INFO  ASEReadCounter - Using google-cloud-java patch 6d11bef1c81f885c26b2b56c8616b7a705171e4f from https://github.com/droazen/google-cloud-java/tree/dr_all_nio_fixes
10:15:10.265 INFO  ASEReadCounter - Initializing engine
10:15:11.286 INFO  FeatureManager - Using codec VCFCodec to read file file:///sample.vcf.gz
10:15:11.442 INFO  ASEReadCounter - Done initializing engine
10:15:11.448 INFO  ProgressMeter - Starting traversal
10:15:11.448 INFO  ProgressMeter -        Current Locus  Elapsed Minutes        Loci Processed      Loci/Minute
10:15:21.460 INFO  ProgressMeter -          chr1:368742              0.2                117000         701228.6
10:15:31.997 INFO  ProgressMeter -          chr1:629976              0.3                213000         621958.3
10:15:42.018 INFO  ProgressMeter -          chr1:698491              0.5                236000         463199.2
10:15:44.969 WARN  ASEReadCounter - Ignoring site: variant is not het at postion: chr1:1302087
10:15:45.124 WARN  ASEReadCounter - Ignoring site: cannot run ASE on non-biallelic sites: [VC Unknown @ chr1:1319980 Q0.00 of type=MIXED alleles=[G*, <*>, A] attr={BQB=1, DP=16, I16=[1, 4, 2, 6, 417, 38245, 629, 53651, 100, 2000, 160, 3200, 110, 2600, 191, 4603], MQ0F=0, MQB=1, MQSB=1, QS=[0.384615, 0.615385, 0], RPB=0.908075, SGB=-0.651104, VDB=0.00836812} GT=PL   85,0,48,100,72,141 filters=
10:15:45.127 WARN  ASEReadCounter - Ignoring site: cannot run ASE on non-biallelic sites: [VC Unknown @ chr1:1320130 Q0.00 of type=MIXED alleles=[G*, <*>, T] attr={BQB=0.477935, DP=17, I16=[3, 4, 5, 3, 470, 31682, 494, 31808, 140, 2800, 160, 3200, 110, 2306, 125, 2683], MQ0F=0, MQB=1, MQSB=1, QS=[0.466667, 0.533333, 0], RPB=0.921243, SGB=-0.651104, VDB=0.890472} GT=PL  83,0,74,104,98,164 filters=
10:15:52.045 INFO  ProgressMeter -         chr1:1566872              0.7                758000        1120279.8
10:15:56.002 WARN  ASEReadCounter - Ignoring site: cannot run ASE on non-biallelic sites: [VC Unknown @ chr1:1796907 Q0.00 of type=MIXED alleles=[C*, <*>, T] attr={BQB=0.5, DP=6, I16=[3, 1, 0, 2, 277, 19213, 124, 7738, 80, 1600, 40, 800, 51, 695, 28, 392], MQ0F=0, MQB=1, MQSB=0.861511, QS=[0.666667, 0.333333, 0], RPB=0.75, SGB=-0.453602, VDB=0.02} GT=PL 19,0,53,31,59,78 filters=
.
.
.
(continues the same way)
.
.
.
10:56:25.174 WARN  ASEReadCounter - Ignoring site: cannot run ASE on non-biallelic sites: [VC Unknown @ chr22:50445452 Q0.00 of type=MIXED alleles=[C*, <*>, T] attr={BQB=0.991201, DP=53, I16=[7, 19, 6, 14, 1786, 128860, 1271, 83753, 520, 10400, 400, 8000, 501, 11051, 386, 8524], MQ0F=0, MQB=1, MQSB=1, QS=[0.565217, 0.434783, 0], RPB=0.95936, SGB=-0.692067, VDB=0.171326} GT=PL  115,0,137,193,197,255 filters=
10:56:25.317 INFO  ProgressMeter -       chr22:50466875             41.2             375434000        9105591.3
10:56:37.922 INFO  ProgressMeter -    KI270733.1:177566             41.4             375720000        9066332.5
10:56:47.925 INFO  ProgressMeter -        chrX:17782656             41.6             377480000        9072304.7
10:56:58.022 INFO  ProgressMeter -        chrX:48266117             41.8             379323000        9079875.6
10:57:08.024 INFO  ProgressMeter -        chrX:56859067             41.9             380706000        9076761.4
10:57:18.025 INFO  ProgressMeter -        chrX:75064330             42.1             382619000        9086261.8
10:57:28.050 INFO  ProgressMeter -       chrX:104532016             42.3             384291000        9089900.6
10:57:38.056 INFO  ProgressMeter -       chrX:133362446             42.4             386383000        9103474.1
10:57:48.187 INFO  ProgressMeter -       chrX:154401268             42.6             388123000        9108235.1
10:57:58.188 INFO  ProgressMeter -         chrY:1392302             42.8             388731000        9086958.6
10:58:08.784 INFO  ProgressMeter -    GL000220.1:115358             43.0             390465000        9089967.3
10:58:18.815 INFO  ProgressMeter -    KI270744.1:125108             43.1             390951000        9065996.4
10:58:28.827 INFO  ProgressMeter -    GL339449.2:846918             43.3             393140000        9081616.5
10:58:38.832 INFO  ProgressMeter -     KI270832.1:10337             43.5             395213000        9094489.1
10:58:49.419 INFO  ProgressMeter -    KI270853.1:403868             43.6             397072000        9100299.4
10:58:59.456 INFO  ProgressMeter -    KI270853.1:665159             43.8             397231000        9069173.3
10:59:09.465 INFO  ProgressMeter -    KI270853.1:942620             44.0             397407000        9038766.6
10:59:19.467 INFO  ProgressMeter -   KI270857.1:1171038             44.1             398301000        9024882.4
10:59:29.480 INFO  ProgressMeter -    KI270875.1:189541             44.3             400159000        9032866.6
10:59:39.725 INFO  ProgressMeter -    KI270897.1:311024             44.5             400922000        9015300.9
10:59:54.064 INFO  ProgressMeter -    KI270897.1:449611             44.7             400979000        8968387.6
11:00:04.076 INFO  ProgressMeter -   GL000251.2:3109525             44.9             401940000        8956457.8
11:00:14.089 INFO  ProgressMeter -   KI270908.1:1250971             45.0             403993000        8968849.4
11:00:24.597 INFO  ProgressMeter -   GL000252.2:4519995             45.2             405265000        8962242.8
11:00:34.598 INFO  ProgressMeter -   GL949749.2:1057485             45.4             406585000        8958412.1
11:00:44.600 INFO  ProgressMeter -     GL000255.2:23082             45.6             407776000        8951774.4
11:00:54.620 INFO  ProgressMeter -    GL949751.2:181561             45.7             408847000        8942501.6
11:01:04.727 INFO  ProgressMeter -     KI270938.1:81508             45.9             410141000        8937873.7
11:01:26.276 INFO  ProgressMeter -            chrM:2078             46.2             410300000        8871901.2
11:01:53.093 INFO  ProgressMeter -            chrM:3078             46.7             410301000        8787001.9
11:02:07.769 INFO  ProgressMeter -            chrM:4078             46.9             410302000        8741233.7
11:02:19.157 INFO  ProgressMeter -            chrM:5078             47.1             410303000        8706051.4
11:02:32.930 INFO  ProgressMeter -            chrM:6078             47.4             410304000        8663873.3
11:03:00.646 INFO  ProgressMeter -            chrM:7078             47.8             410305000        8580202.6
11:03:43.980 INFO  ProgressMeter -            chrM:8078             48.5             410306000        8452562.9
11:04:15.210 INFO  ProgressMeter -            chrM:9078             49.1             410307000        8362911.1
11:04:41.972 INFO  ProgressMeter -           chrM:10078             49.5             410308000        8287588.3
11:04:57.510 INFO  ProgressMeter -           chrM:11078             49.8             410309000        8244483.9
11:05:26.815 INFO  ProgressMeter -           chrM:12078             50.3             410310000        8164379.3
11:05:38.066 INFO  ProgressMeter -           chrM:14078             50.4             410312000        8134069.1
11:05:48.378 INFO  ProgressMeter -           chrM:15078             50.6             410313000        8106469.4
11:05:58.376 INFO  ASEReadCounter - No reads filtered by: AllowAllReadsReadFilter
11:05:58.376 INFO  ProgressMeter -           chrM:16078             50.8             410314491        8079898.7
11:05:58.376 INFO  ProgressMeter - Traversal complete. Processed 410314491 total loci in 50.8 minutes.
11:05:58.376 INFO  ASEReadCounter - Shutting down engine
[July 5, 2018 11:05:58 AM GMT] org.broadinstitute.hellbender.tools.walkers.rnaseq.ASEReadCounter done. Elapsed time: 50.81 minutes.
Runtime.totalMemory()=2045640704
Using GATK jar /homeappl/appl_taito/bio/GATK4/gatk-4.0.2.1/gatk-package-4.0.2.1-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=1 -DGATK_STACKTRACE_ON_USER_EXCEPTION=true -jar /homeappl/appl_taito/bio/GATK4/gatk-4.0.2.1/gatk-package-4.0.2.1-local.jar ASEReadCounter --input sample_readgroups.bam --variant sample.vcf.gz --disable-tool-default-read-filters true --reference hg38.fa --output sample.ASEReadCounter.rtable

Output rtable contains only header:

contig  position    variantID   refAllele   altAllele   refCount    altCount    totalCount  lowMAPQDepth    lowBaseQDepth   rawDepth    otherBases  improperPairs


What should I do?

The same thing happens also when I enable default read filters.

Here is a my vcf:

##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##samtoolsVersion=1.3.1+htslib-1.3.1
##samtoolsCommand=samtools mpileup -v -u -f hg38.fa -l sample.bed -o sample.vcf sample_sorted.bam
##reference=file:///hg38.fa
##contig=<ID=chr1,length=248956422>
##contig=<ID=KI270706.1,length=175055>
##contig=<ID=KI270707.1,length=32032>
##contig=<ID=KI270708.1,length=127682>
##contig=<ID=KI270709.1,length=66860>
##contig=<ID=KI270710.1,length=40176>
##contig=<ID=KI270711.1,length=42210>
##contig=<ID=KI270712.1,length=176043>
##contig=<ID=KI270713.1,length=40745>
##contig=<ID=KI270714.1,length=41717>
##contig=<ID=chr2,length=242193529>
##contig=<ID=KI270715.1,length=161471>
##contig=<ID=KI270716.1,length=153799>
##contig=<ID=chr3,length=198295559>
##contig=<ID=GL000221.1,length=155397>
##contig=<ID=chr4,length=190214555>
##contig=<ID=GL000008.2,length=209709>
##contig=<ID=chr5,length=181538259>
##contig=<ID=GL000208.1,length=92689>

.
.
.
(continues the same way)
.
.
.
##contig=<ID=KI270930.1,length=200773>
##contig=<ID=KI270931.1,length=170148>
##contig=<ID=KI270932.1,length=215732>
##contig=<ID=KI270933.1,length=170537>
##contig=<ID=GL000209.2,length=177381>
##contig=<ID=chrM,length=16569>
##ALT=<ID=*,Description="Represents allele(s) other than observed.">
##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
##INFO=<ID=IDV,Number=1,Type=Integer,Description="Maximum number of reads supporting an indel">
##INFO=<ID=IMF,Number=1,Type=Float,Description="Maximum fraction of reads supporting an indel">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
##INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias for filtering splice-site artefacts in RNA-seq data (bigger is better)",Version="3">
##INFO=<ID=RPB,Number=1,Type=Float,Description="Mann-Whitney U test of Read Position Bias (bigger is better)">
##INFO=<ID=MQB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality Bias (bigger is better)">
##INFO=<ID=BQB,Number=1,Type=Float,Description="Mann-Whitney U test of Base Quality Bias (bigger is better)">
##INFO=<ID=MQSB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality vs Strand Bias (bigger is better)">
##INFO=<ID=SGB,Number=1,Type=Float,Description="Segregation based metric.">
##INFO=<ID=MQ0F,Number=1,Type=Float,Description="Fraction of MQ0 reads (smaller is better)">
##INFO=<ID=I16,Number=16,Type=Float,Description="Auxiliary tag used for calling, see description of bcf_callret1_t in bam2bcf.h">
##INFO=<ID=QS,Number=R,Type=Float,Description="Auxiliary tag used for calling">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
##bcftools_normVersion=1.4+htslib-1.4
##bcftools_normCommand=norm -d any -O z -o sample.vcf.gz sample.vcf.gz; Date=Thu Jul  5 10:33:29 2018
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample_sorted.bam
chr1    1302087 .       G       <*>     0       .       DP=3;I16=0,2,0,0,139,9661,0,0,40,800,0,0,49,1201,0,0;QS=1,0;MQ0F=0      PL      0,6,36
chr1    1319980 .       G       A,<*>   0       .       DP=16;I16=1,4,2,6,417,38245,629,53651,100,2000,160,3200,110,2600,191,4603;QS=0.384615,0.615385,0;VDB=0.00836812;SGB=-0.651104;RPB=0.908075;MQB=1;MQSB=1;BQB=1;MQ0F=0    PL      85,0,48,100,72,141
chr1    1320130 .       G       T,<*>   0       .       DP=17;I16=3,4,5,3,470,31682,494,31808,140,2800,160,3200,110,2306,125,2683;QS=0.466667,0.533333,0;VDB=0.890472;SGB=-0.651104;RPB=0.921243;MQB=1;MQSB=1;BQB=0.477935;MQ0F=0       PL      83,0,74,104,98,164
.
.
.
(etc.)

I am out of ideas. Thank you so much for your time.

Best Answer

Answers

  • Hi @shlee ,

    Thank you for the answer! I see now why it did not work. I will give it a new try using your suggestions and the human data from GATK Resource Bundle. Those multiple alt alleles were probably generated by Samtools mpileup so I will also need to check that pipeline again. The symbolic alternate alleles are not actually necessary for my analysis.

    Thank you! Let's see if it will work now. :)

  • shleeshlee CambridgeMember, Broadie, Moderator admin

    Let us know how it goes @SofiaRandelin.

  • Hi @shlee ,

    I tried with the VCF file from GATK Resource Bundle and ASEReadCounter is now ignoring many more sites and the output table has only header. There are multiple lines like this:

    14:03:50.586 WARN  ASEReadCounter - Ignoring site: variant is not het at postion: chrX:153783936
    

    I also tried the latest version GATK 4.0.5.2. Both of my runs had some lines/sites that were not ignored (below) but they are not written in the output file. Should these positions be in the output table:

    14:04:08.747 INFO  ProgressMeter -        chrY:18327568             55.2             389599000        7061426.0
    

    Could the problem be in my BAM file or in the way I processed it?

Sign In or Register to comment.