RADseq data and GATK

mlpimslermlpimsler University of AlabamaMember

Hello!
I am pretty new to bioinformatics- mostly just have taken one class.

I have been handed some RADseq data and we want to do some SNP variant calling. I am trying to make sure that I use GATK correctly- at present we only have about 250 libraries, but we should end up with at least twice that many, possibly three times.

I have aligned the reads to a reference genome, sorted and indexed them, and am now ready for HaplotypeCaller. We do not have SNP location information, so I was going to follow the Troubleshooting guide by using GATK to de novo call SNPs, feed that information to recalibrate based on quality scores, and then re-run Haplotype caller. I did not do the dedup step due to the nature of the data.

I am trying to make sure that I am using the correct options. What I have is:

java -jar GenomeAnalysisTK.jar \
-T HaplotypeCaller \
-R reference.fa \
-I preprocessed_reads.bam \
--genotyping_mode DISCOVERY \
-stand_emit_conf 10 \
-stand_call_conf 30 \
-o raw_variants.vcf

Any help would be greatly appreciated

Tagged:

Best Answers

Answers

  • mlpimslermlpimsler University of AlabamaMember
    edited September 2015

    I have, at present, 250 separate libraries. I will have about 1000 libraries by the end of the project.

    I have a script that goes through every file in a for loop and does all of the processing up to creating the .g.vcf files for each RAD-tag library.

    I have not been able to successfully batch-call all of the .g.vcf files that have been created in combineGVCF, however, without creating a line for each file individuals. I have tried using the wildcard symbol (*) I have seen in other help threads (eg. http://gatkforums.broadinstitute.org/discussion/5469/how-to-speed-up-combinegvcfs-seems-unfeasibly-slow ) but it didn't work.

    All of my files are named: "{libraryname}-raw_variants.g.vcf". I'm testing out five of the libraries at first to get used to the pipeline before pushing forward with all of the data.

    I have five files (library1-raw_variants.g.vcf, library2-raw_variants.g.vcf, etc) all in the directory I am calling the command from, and the reference .fasta and associated indexes and dictionaries in the parent directory.

    java -jar GenomeAnalysisTK.jar \
    -T CombineGVCFs \
    -R ../reference.fasta \
    -drf DuplicateRead \
    --variant *raw_variants.g.vcf \
    -o $DIR.g.vcf

    The output I get from the program includes the following comment:

    ERROR A USER ERROR has occurred (version 3.4-46-gbc02625):
    ERROR MESSAGE: Invalid argument value 'library2-raw_variants.g.vcf' at position 8.
    ERROR Invalid argument value 'library3-raw_variants.g.vcf' at position 9.
    ERROR Invalid argument value 'library4-raw_variants.g.vcf' at position 10.
    ERROR Invalid argument value 'library5-raw_variants.g.vcf' at position 11.

    None of the documentation here, or provided with the --help option is clear about how whether each individual variant must be on it's own line.

    Any help would be appreciated.

  • mlpimslermlpimsler University of AlabamaMember

    Oh, great, okay. I will create that list file and see how it goes. Thanks!

  • LuisaBLuisaB SwitzerlandMember

    Hi,

    I am analyzing single-end RAD-seq data and I have a question regarding the steps after the mapping. I used Bowtie2 to align my reads to the reference genome and I think in the resulting .bam files (one for each individual) all reads are kept, I mean those that mapped in a unique position, those that mapped in multiple positions and those that didn't map at all.
    My question is the following: should I keep in my .bam files only the reads that map uniquely and remove the others? Or the reads that map multiple times and those that didn't map at all have such a low mapping quality that they will not be taken into account by the tools I will use downstream (RealignerTargetCreator, IndelRealigner, BaseRecalibrator, HaplotypeCaller)?

    Thank you very much!

    Luisa

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @LuisaB
    Hi Luisa,

    You can run ValidateSamFile on your bam file. If there are no errors, the GATK tools will take care of the unmapped reads or reads that do not map to only one position. http://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile

    -Sheila

  • LuisaBLuisaB SwitzerlandMember

    Hi @Sheila,

    thank you for your reply.
    What do you mean when you say that "the GATK tools will take care" of these reads? Do you mean they will not be taken into account?

    Thanks again!

    Luisa

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @LuisaB
    Hi Luisa,

    The tools have their own filters to simply filter out the reads that might cause an issue (they will not taken into account :smile: ). You can see which filters are applied in the individual tool documentation. For more information on the read filters, you can look them up under Read Filters here: https://www.broadinstitute.org/gatk/guide/tooldocs/

    -Sheila

  • LuisaBLuisaB SwitzerlandMember

    Hi @Sheila,

    I used ValidateSamFile on some of my .bam files and the output says "No errors found", so they should be ok.

    But there is still something unclear to me. Here is the end of the standard error produced by HaplotypeCaller (I ran it on one chromosome on 520 individuals all together - not in -ERC mode because I was not aware of this possibility yet):

    INFO 13:11:15,232 MicroScheduler - 11991669 reads were filtered out during the traversal out of approximately 39370204 total reads (30.46%)
    INFO 13:11:15,232 MicroScheduler - -> 0 reads (0.00% of total) failing DuplicateReadFilter
    INFO 13:11:15,232 MicroScheduler - -> 0 reads (0.00% of total) failing FailsVendorQualityCheckFilter
    INFO 13:11:15,232 MicroScheduler - -> 11991669 reads (30.46% of total) failing HCMappingQualityFilter
    INFO 13:11:15,232 MicroScheduler - -> 0 reads (0.00% of total) failing MalformedReadFilter
    INFO 13:11:15,232 MicroScheduler - -> 0 reads (0.00% of total) failing MappingQualityUnavailableFilter
    INFO 13:11:15,232 MicroScheduler - -> 0 reads (0.00% of total) failing NotPrimaryAlignmentFilter
    INFO 13:11:15,233 MicroScheduler - -> 0 reads (0.00% of total) failing UnmappedReadFilter
    INFO 13:11:16,399 GATKRunReport - Uploaded run statistics report to AWS S3

    Since in my .bam files there should be uniquely mapping reads, reads mapping in multiple locations and unmapped reads, I would expect the filters NotPrimaryAlignmentFilter and UnmappedReadFilter to filter out some percentage of the reads, but this doesn't seem to happen. Only HCMappingQualityFilter filters out some reads.
    Do you have any idea of the reason?

    Thanks again,

    Luisa

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @LuisaB
    Hi Luisa,

    How do you know the reads that should be filtered out exist in your bam file? Can you tell me the counts of reads you think should be filtered out? It is possible those reads might have been removed from the data in a previous step.

    -Sheila

  • mlpimslermlpimsler University of AlabamaMember

    I was observing about 11% of reads per library being filtered out in the indel identification steps (-T RealignerTargetCreator). This is approximately the percentage of reads in each library which multiply mapped or didn't align at all.

    INFO 15:49:33,278 MicroScheduler - 178046 reads were filtered out during the traversal out of approximately 1591954 total reads (11.18%)
    INFO 15:49:33,279 MicroScheduler - -> 0 reads (0.00% of total) failing BadCigarFilter
    INFO 15:49:33,279 MicroScheduler - -> 0 reads (0.00% of total) failing BadMateFilter
    INFO 15:49:33,279 MicroScheduler - -> 0 reads (0.00% of total) failing DuplicateReadFilter
    INFO 15:49:33,280 MicroScheduler - -> 0 reads (0.00% of total) failing FailsVendorQualityCheckFilter
    INFO 15:49:33,280 MicroScheduler - -> 0 reads (0.00% of total) failing MalformedReadFilter
    INFO 15:49:33,280 MicroScheduler - -> 0 reads (0.00% of total) failing MappingQualityUnavailableFilter
    INFO 15:49:33,281 MicroScheduler - -> 176102 reads (11.06% of total) failing MappingQualityZeroFilter
    INFO 15:49:33,281 MicroScheduler - -> 1944 reads (0.12% of total) failing NotPrimaryAlignmentFilter
    INFO 15:49:33,281 MicroScheduler - -> 0 reads (0.00% of total) failing Platform454Filter
    INFO 15:49:33,281 MicroScheduler - -> 0 reads (0.00% of total) failing UnmappedReadFilter

    However, similar to Luisa, each file ALSO had about 20- 25% of reads filtered out in the .g.vcf creation (-T HaplotypeCaller).

    INFO 17:11:55,619 MicroScheduler - 367557 reads were filtered out during the traversal out of approximately 1491136 total reads (24.65%)
    INFO 17:11:55,620 MicroScheduler - -> 0 reads (0.00% of total) failing BadCigarFilter
    INFO 17:11:55,620 MicroScheduler - -> 0 reads (0.00% of total) failing FailsVendorQualityCheckFilter
    INFO 17:11:55,620 MicroScheduler - -> 366074 reads (24.55% of total) failing HCMappingQualityFilter
    INFO 17:11:55,620 MicroScheduler - -> 0 reads (0.00% of total) failing MalformedReadFilter
    INFO 17:11:55,621 MicroScheduler - -> 0 reads (0.00% of total) failing MappingQualityUnavailableFilter
    INFO 17:11:55,621 MicroScheduler - -> 1483 reads (0.10% of total) failing NotPrimaryAlignmentFilter
    INFO 17:11:55,621 MicroScheduler - -> 0 reads (0.00% of total) failing UnmappedReadFilter

    In my case, research is on a non-model organism for which genomic data is not available, so we are using a closely related sequences sibling species. However, it should be noted that this genus demonstrates significant conservation in synteny and sequence with the other sequenced species in this genus.

    However, I am still concerned that I am losing ~30% of the reads through the GATK process.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @mlpimsler You can lower the threshold for the HC mapping quality filter using the -mmq argument if you want to increase the proportion of reads that get used, since in your case that's the main culprit for filtering reads.

  • LuisaBLuisaB SwitzerlandMember

    Hi @Sheila,

    I know the reads are still in the bam files because I compared the number of reads I have for one individual in the fastq file before mapping, and the number of reads I have in the bam file (with samtools flagstat), and they are the same. For example, I have 1604521 reads totally in the fastq file and the output of samtools flagstat on the bam file is the following:

    1604521 + 0 in total (QC-passed reads + QC-failed reads)
    0 + 0 duplicates
    1360892 + 0 mapped (84.82%:-nan%)
    0 + 0 paired in sequencing
    0 + 0 read1
    0 + 0 read2
    0 + 0 properly paired (-nan%:-nan%)
    0 + 0 with itself and mate mapped
    0 + 0 singletons (-nan%:-nan%)
    0 + 0 with mate mapped to a different chr
    0 + 0 with mate mapped to a different chr (mapQ>=5)

    The total number is the same (1604521) and 1360892 corresponds to the sum of the reads that align uniquely and the reads that align in multiple locations. The numbers are still the same also in the bam file I get after the IndelRealignment.

    I don't know precisely the count of the reads that should be filtered out in the example I posted before. But from the results of the mapping, I know that, for each of my individuals, on average 55% of the reads map uniquely, 30% map in multiple locations and 15% don't map.

    Thanks for your help,

    Luisa

    Issue · Github
    by Sheila

    Issue Number
    186
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    vdauwera
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @LuisaB,

    I think I know what's happening here. There can be some overlap between what some filters use as filtering criteria, and depending on the order in which the filters are applied, if reads that could be filtered by two different filters get filtered out by one, they won't be seen (or counted) by the other. In this case, the HCMappingQuality filter is being applied before the other filters, and since it operates on any read that has MQ lower than a given threshold (20 by default), including 0, it is filtering out reads that have low or zero MQ because they map poorly or non-uniquely.

    For the subset of reads that don't map at all, I think there might be some special-casing that makes the tool not run on them at all so they don't even get seen by the filters. I think I remember some debate about whether that's the right way to handle them or not. I'll look into it. But I don't think there's anything here you need to worry about.

  • LuisaBLuisaB SwitzerlandMember

    Hi @Geraldine_VdAuwera,

    thank you for your answer, it makes totally sense for HaplotypeCaller (also if I am not sure the percentage of reads filtered out by the HCMappingQualityFilter matches exactly my expectations). But I'm not sure this is the reason for other tools: when I run for example the RealignerTargetCreator on all my individuals on all the chromosomes, only 6.92% of the reads are filtered out because failing MappingQualityZeroFilter, but the other filters (NotPrimaryAlignmentFilter and UnmappedReadFilter) didn't filter out any read. It can be that the reads that should be filtered out by these two filters have already been filtered by the first one, but I would expect a much higher percentage than 6.92% (according to the percentages given by Bowtie2).
    And if I run UnifiedGenotyper, no reads are filtered out by all the filters.
    I would really like to be sure that only good reads are taken into account in the analysis.

    Thanks a lot,

    Luisa

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    If you want to know exactly what should get filtered by each filter independently, you can use PrintReads with the -rf argument to specify one filter at a time and generate filtered bams. But to be honest I think it would be a waste of time -- the read filtering is well tested and we're pretty confident that the engine is doing the right thing. I agree it does seem confusing on the face of it, but there's usually a very simple explanation for this sort of thing, and unfortunately we really can't take the time to look at your stats in detail. I understand this is probably not a very satisfying answer but it's all I can do at the moment.

  • BegaliBegali GermanyMember
    edited June 2017

    @Geraldine_VdAuwera
    @Sheila
    hi

    I have one question as I am working on RADSeq and I do not have knownsite_indel so it is nesseray to do indel-realignment step or not because I call my variants directly after preprocesses steps ( map my reads.sample by BWA with mem for single read to obtain sam .then sort index my bam.file ) , also I can not apply BSQR so skip those will not affect my result at the end ..aim of my project to discover Varaints at the end and link them to phenotypes of plants for multiple samples ..

    Thanks so much

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin
    Indel realignment is no longer necessary if you use a caller that includes a graph assembly step like HaplotypeCaller.

    BQSR should still be done, and it is possible to bootstrap a set of known variants if you don't have one. By skipping it, you pass on the opportunity to improve your data and the quality of your results. However we have no way to quantify the effect since it varies from one sequencing run to the next. So it's up to you to decide what you want to do.
  • BegaliBegali GermanyMember

    hi
    @Geraldine_VdAuwera

    Thanks so much for your information and reply .. :smile: )

  • BegaliBegali GermanyMember

    Hi @Geraldine_VdAuwera

    Thanks so much and I am sorry to write alot as I am new to deal with anaylsis seq and those tools can you please kindly give me your hints for my confusion
    I was calling variants and got my vcf.file by four methods (samtools with mpileup , UnifiedGenotyper ,HaplotypeCaller , in GVCF mode ) so with haplotyecaller as u mentioned no need for indel realignment however for othres methosd it is still necessary ????? ..
    if yes I need to apply for others methods (samtools with mpileup , UnifiedGenotyper , in GVCF mode ) so in case as you said also I do not have known variants I mean miss this parameter known indels.vcf , indelrealignments should work well even if no known variants available (plant Arabidopsis thaliana ) .. another question also for BQSR in two its steps (RealignerTargetCreator and IndelRealigner ) still I do not have known indels.vcf folder for my data can I apply BQSR without it it still work perfect ...also BQSR step should apply before each methods for calling variants (samtools with mpileup , UnifiedGenotyper ,HaplotypeCaller , in GVCF mode ) because my project goal is to call variants for RADSeq for plant (197 samples) then compare which method indicates correct and high quality and good for RADSeq and our organism at the end I will plot my original vcf.file before filtering to determine thresholds which I would like to remove actually one with lowerest score that what I understood then filter them by only this option for RADSeq hard filtering ...

    Thanks so much for your reply in advance ,,,

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Begali
    Hi,

    Because HaplotypeCaller does a local reassembly, we do not recommend doing Indel Realignment. For other variant callers such as Unified Genotyper which does not do a local reassembly, we do recommend doing Indel Realignment. Have a look at this blog post for more information.

    You can still perform Indel Realignment without a known indels file, but it will increase the runtime.

    For BQSR, have a look at this article under "I'm working on a genome that doesn't really have a good SNP database yet. I'm wondering if it still makes sense to run base quality score recalibration without known SNPs." Note that you can use the known indels you get from BQSR in Indel Realignment.

    -Sheila

  • JenniJenni Member

    Hi,
    I'm just using GATK for the first time and am still fairly new to bioinformatics. I have a set of ddRAD-Seq data, which I'm wanting to use HaplotypeCaller for.
    So far I have aligned it to my reference genome (in >81000 contigs) using BWA MEM and am now waiting on a run of the Haplotype caller. I don't have any 'known' SNPs, other than those which have been called by Julian Catchen's Stacks program. And, these are likely to be all of the variants considered worthy of being called.
    So I have a question about the BQSR: As I have no input file I could call variants using HaplotypeCaller and then use the variants called as input for BQSR, run this and then re-run HaplotypeCaller, is this correct?
    But, my question is, as I don't fully understand it all, why would I want to do this? Is BQSR adjusting quality scores of all bases other than known variants, in order that HaplotypeCaller can then look at these bases (unknown, and assumed errors) and call real variants from them? Or is it more to do with accepting or rejecting reads to use in the analysis based on overall read quality?
    Thanks!
    Jenni

  • JenniJenni Member

    Also, could I just ask a very basic question... I have input all of my 62 bam files (each for a different individual), one after the other in the same script, but have just one output file. If I was instead to ask to generate gVCF files (using -ERC) and then combine the gVCFs and re-analyse jointly to output one vcf file, should I run each sample separately. Ie. submit 62 jobs?
    Thanks!

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Jenni
    Hi Jenni,

    As I have no input file I could call variants using HaplotypeCaller and then use the variants called as input for BQSR, run this and then re-run HaplotypeCaller, is this correct?

    Yes, that is correct. Have a look at this thread.

    I hope this article will answer your other questions in your first post.

    For your second post, have a look at this article and this article.

    -Sheila

Sign In or Register to comment.