Attention:
The frontline support team will be unavailable to answer questions on April 15th and 17th 2019. We will be back soon after. Thank you for your patience and we apologize for any inconvenience!

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 ✭✭✭✭✭

    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.