We've moved!
This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!

Questions about input files

Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin
This discussion was created from comments split from: What input files does the GATK accept / require?.

Comments

  • zugzugzugzug Member ✭✭

    Can you provide a link for "Preparing the essential GATK input files: the reference genome?"

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Added in the article-- thanks for reporting this missing link.

  • zugzugzugzug Member ✭✭

    It's also missing under Guide -> Introductory Materials -> What input files does the GATK accept?

    Thanks for the quick response!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    That's actually the same article (the Guide "pulls" articles from the forum), but it's a cached version that hasn't been refreshed yet -- will be soon though.

    I'm here to help :)

  • mikemike Member

    Hi, Geraldine:

    the document as above mentioned: The contig ordering in the reference you used must exactly match that of one of the official references canonical orderings. These are defined by historical karotyping of largest to smallest chromosomes, followed by the X, Y, and MT. The order is thus 1, 2, 3, ..., 10, 11, 12, ... 20, 21, 22, X, Y, MT. However, I noticed that one of training data at VQSR step: dbsnp_135.hg19.vcf has chrM at beginning, and so I run into issue at VQSR step. When I changed my reference chromosome order as chrM, chr1, chr2,…chr10,chr11,…….chr22,chrX, chrY, then it works fine. Could you comment on this behavior?

    Thanks,

    Mike

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Mike,

    The order given in the document refers to the ordering we use preferentially, which corresponds to the b37 version of the human reference genome. The version of dbsnp you used is ordered to suit the UCSC-style reference genome, hg19 (which you can recognize in the filename).

    You should always choose your reference version very carefully at the start of your analysis. If your data is unaligned, you can usually choose whatever you feel like. But once your data is already aligned to a reference, you should pick the corresponding reference, and stick with it throughout your analysis. That means that when it comes to using dbsnp, for example, you have to use the right version that matches your reference (again, the filename will tell you which reference it corresponds to). We provide both versions in our resource bundle.

    If you don't do this, or try to mix & match halfway through, you're gonna have a bad time.

  • mikemike Member

    Thanks for the input! Mike

  • When building picard-style Interval Lists, must my targeted intervals be sorted according to the chromosome which they target? Or can I list the targeted intervals in any order?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Off the top of my head I would say that they don't have to be ordered, but I'll have to check. You can do a very simple test by feeding in a few unordered intervals; you'll see quickly if the program chokes or not...

  • ebanksebanks Broad InstituteMember, Broadie, Dev ✭✭✭✭

    That's right - the GATK will sort intervals on the fly.

  • sej1985sej1985 Member

    Do you have a link pointing to this:::
    See [here] for converting .ped files to VCF

  • sboylesboyle Bay Area, CAMember

    I have a question about the chromosome ordering in the reference file. When I go into your bundle (/bundle/1.5/hg19) and look at the ucsc.hg19.fasta file located there I find that it begins with >chrM and not >chr1. On this page it mentions that our reference file should go from 1-22, X, Y, M. I have assembled my own genome.fa from the ensemble chromosome files that follows your suggested order and am creating an index file, but I wonder whether I should just be using the available ensemble genome.fa and precomputed BWA index file that follows the same order of your ucsc.hg19.fasta? Do you have any suggestions to offer on these two options? Thanks in advance for the assistance!

  • ebanksebanks Broad InstituteMember, Broadie, Dev ✭✭✭✭

    I fixed the text above to differentiate between b3x and hg1x references.

  • mikemike Member
    edited November 2012

    Hi,

    what if my reference has some contigs like

    >chr4_gl000194_random
    >chr4_gl000193_random
    >chr9_gl000200_random
    

    which is not in your dbSNP file, does it matter how they are ordered as long as put them all at the end of the reference file after all the regular chromosomes such as chrM, chr1, chr2, chr21, chr22, chrX, chrY?

    Any suggestions?

    Thanks!

    Mike

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    That should be fine, Mike.

  • mikemike Member

    Dear Geraldine_VdAuwera:

    Thanks for the comments!

    Mike

  • Hi Geraldine,

    FYI, the link to VCF poster specifications is missing (and I was quite curious about it, I must say ;) )

    best
    Andrea

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Andrea, thanks for pointing this out. It seems we have a couple of missing links there, sorry about that. I'll try to restore them asap; in the meantime the best resources for learning about VCFs are at the 1000 Genomes project website and the VCFTools website, which I've added links for in that paragraph.

  • vyellapavyellapa Member

    I was trying to use the Select Variants to filter variants over a specified interval list. I understand that the ".interval_list" is picard interval style. Usually I use the header from the bam as the header for this file. If I do not have the bam, is there a way to create .interval_list?

    Thank you,
    Teja

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin
    edited February 2013

    Hi Teja,

    The intervals list doesn't need a header. Just write your intervals, one per line, in the following format:
    1:10000-20000

    See the documentation / FAQs on input formats for more information.

  • GarethLGarethL Member

    Hi Geraldine

    I'm just about to start using GATK, but you say the reference must only contain the standard bases, ATCG. What about N's?

    Our reference is full of iupac codes, but if N's are ok I will just swap them all for N's and see how we go.

    thanks

    Gareth

  • CarneiroCarneiro Charlestown, MAMember admin

    N is totally fine.

  • GarethLGarethL Member
  • igcocoleigcocole Member

    Hi,
    You mention that "no non-standard bases (W, for example) are tolerated": why does the GATK bundle b37 fasta file contains R and M bases in chr 3 (human_g1k_v37.fasta.gz file) if this is mentioned as an issue?

    Can that influence the analysis somehow?

    Thank you!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @igcocole,

    This is a case of outdated doc -- the GATK now handles all IUPAC bases, though by "handles" I mean "ignores" (whereas the previous behavior was to blow up). I will correct the doc, thanks for pointing this out.

  • drunkcoderdrunkcoder Member

    Hi, Geraldine
    I use newest hg19 from UCSC as reference, there are contigs like chr6_apd_hap1.fa, chr8_gl000197_random.fa and unplaced contigs like chrUn_gl000221.fa ... what order should these contigs be placed? or just put them all at the end of the reference file after all the regular ?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @drunkcoder,

    At the most basic level, the important thing is to have your bam sorted in the same order as the reference you use. The convention we use is to sort the reference lexicographically.

    In practice however -- keep in mind that if you're using the latest version of the reference, you may not be able to use the resource files we provide (snps, known indels etc). If the new version just has some additional contigs then it may be compatible as long as the extra contigs are placed at the end. But if the reference has been patched in the regular contigs then our resources might not be compatible.

  • jlinjlin Member

    I had the error about the non-iupac "10" character, but i have to start with bam files that were aligned to another reference file and thus can't use the ones in the gatk bundle. You said that we could clean up the reference file ourselves... how could I do that? Thanks in advance!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    The simplest way is to write a script that will recognize and remove the bad characters. We do not provide guidance on how to do that, since it's a general programming question that is not related to the GATK, but perhaps someone in the community might be able to share some advice. Otherwise I would recommend asking your IT department or a colleague with relevant experience.

  • jlinjlin Member

    I talked to some people with more programming experience than me, and they can't figure out how to clean up the reference file. I had already tried common methods of converting windows to unix files (e.g. dos2unix) [although my reference files had never been opened by a windows program].

    Since I can't realign my files to the GATK bundle reference file, I guess I'm going to have to go to an older version of GATK to bypass this error, although I know that's not ideal.

    Just as a suggestion, could you look at changing the error message so that it indicates exactly which file has the non-iupac base (i.e. the fasta file, the index file, or the dictionary file) and which line the non-iupac base is found on? That would give us a place to start in cleaning up the file.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Jilin,

    The reference file is just a text file, so you don't really need to convert anything. The general idea is to run a simple character substitution function to remove the bad characters.

    This problem is not going to go away if you use an older version of GATK -- the program is always going to choke on the bad characters.

    If you're working with standard human data you could just revert your bams, remap them to a clean reference and proceed normally.

    Finally, the error message does state what file the problem is in. If you post the full error message you got I can help you find where this is written.

  • jlinjlin Member

    Geraldine,

    Your help would be appreciated! I can easily substitute characters with sed, but I'm not sure what characters I'm looking for. I searched for '10,' but the only 10s that come up are for chromosome number (and bp position in the index file).

    The older version of GATK did actually run on my reference file... so maybe I'm interpreting this error wrong?

    How long would it take to revert bams and remap them to the ref file from the GATK bundle?

    The full error message was as follows:

    ERROR ------------------------------------------------------------------------------------------
    ERROR A USER ERROR has occurred (version 2.6-2-ge03a5e9):
    ERROR The invalid arguments or inputs must be corrected before the GATK can proceed
    ERROR Please do not post this error to the GATK forum
    ERROR
    ERROR See the documentation (rerun with -h) for this tool to view allowable command-line arguments.
    ERROR Visit our website and forum for extensive documentation and answers to
    ERROR commonly asked questions http://www.broadinstitute.org/gatk
    ERROR
    ERROR MESSAGE: Bad input: We encountered a non-standard non-IUPAC base in the provided reference: '10'
    ERROR ------------------------------------------------------------------------------------------

    The command was: java -jar GenomeAnalysisTK.jar -R hg19.genome.new.fa -T UnifiedGenotyper -I bamlist.list -L chr1.bed -o chr1.vcf -U

    Thanks for your help!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hmm, which older version did you use?

    I have to correct what I said earlier -- simple character substitution isn't going to work since the problem is in the file encoding, not content -- my mistake, I was thinking of a different problem.

    I haven't ever needed to do this myself so I can't give you exact instructions, but this should be a good starting point: http://stackoverflow.com/questions/64860/best-way-to-convert-text-files-between-character-sets

    The error message does state that it's your reference file which has a problem:

    We encountered a non-standard non-IUPAC base in the provided reference: '10'

    Reverting and remapping can take quite a while if you have a lot of data. If you figure out how to solve this encoding problem, you will be able to deal with it quickly if it happens again. I recommend trying to fix your reference. But it's your choice.

  • jlinjlin Member

    Thanks for the link, I will check it out. The older version is GATK 1.4.

    Yes, I new remapping would take time that we really don't have for this project unfortunately, so I will press on with the older version if I can't fix the reference file!

    Thanks.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    It's up to you, but I have to caution you that using version 1.4 is a really bad idea. The software has improved enormously since then and I can guarantee you would have much better results with the latest version.

  • frankoffrankof Member

    Dear Geraldine,
    I downloaded the human_g1k_v37 from Broad bundle and I created the database for gsnap aligner. After obtaining my .bam (sorted and indexd and after duplicate removal) GATK returns me a chromosome order error. it is referreder only to non canonical (chromosome like GL000207.1, GL000226.1). How can I order my bam file? Is the human_g1k_v37 sorted and indexed with samtools? Usually I use samtools to sort and index my bam.
    Thanks in advance,
    Francesco

  • oliviajmoliviajm Member

    Hi,

    I've noticed that the intervals given in the .interval_list file are grouped when they overlap or when the start of one interval is contiguous to the end of another one. So I have 200 lines in my interval_list file, but I get only 100 in the sample_interval_summary.
    Is there a way to avoid this behaviour to get the 200 intervals separately in the sample_interval_summary file?

    Olivia

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @frankof,

    The reference we provide in our bundle is sorted and indexed according to the rules laid out in our documentation, you shouldn't need to do any additional sorting if you followed the instructions we provide. Can you post the text of the error you got?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @oliviajm,

    No, the interval merging behavior cannot be overridden, sorry.

  • frankoffrankof Member

    Hi Geraldine.
    Thank for your prompt response
    the error is: ##### ERROR MESSAGE: Input files reads and reference have incompatible contigs: Relative ordering of overlapping contigs differs, which is unsafe.

    ERROR reads contigs = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, X, Y, MT, GL000191.1, GL000192.1, GL000193.1, GL000194.1, GL000195.1, GL000196.1, GL000197.1, GL000198.1, GL000199.1, GL000200.1, GL000201.1, GL000202.1, GL000203.1, GL000204.1, GL000205.1, GL000206.1, GL000207.1, GL000208.1, GL000209.1, GL000210.1, GL000211.1, GL000212.1, GL000213.1, GL000214.1, GL000215.1, GL000216.1, GL000217.1, GL000218.1, GL000219.1, GL000220.1, GL000221.1, GL000222.1, GL000223.1, GL000224.1, GL000225.1, GL000226.1, GL000227.1, GL000228.1, GL000229.1, GL000230.1, GL000231.1, GL000232.1, GL000233.1, GL000234.1, GL000235.1, GL000236.1, GL000237.1, GL000238.1, GL000239.1, GL000240.1, GL000241.1, GL000242.1, GL000243.1, GL000244.1, GL000245.1, GL000246.1, GL000247.1, GL000248.1, GL000249.1]
    ERROR reference contigs = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, X, Y, MT, GL000207.1, GL000226.1, GL000229.1, GL000231.1, GL000210.1, GL000239.1, GL000235.1, GL000201.1, GL000247.1, GL000245.1, GL000197.1, GL000203.1, GL000246.1, GL000249.1, GL000196.1, GL000248.1, GL000244.1, GL000238.1, GL000202.1, GL000234.1, GL000232.1, GL000206.1, GL000240.1, GL000236.1, GL000241.1, GL000243.1, GL000242.1, GL000230.1, GL000237.1, GL000233.1, GL000204.1, GL000198.1, GL000208.1, GL000191.1, GL000227.1, GL000228.1, GL000214.1, GL000221.1, GL000209.1, GL000218.1, GL000220.1, GL000213.1, GL000211.1, GL000199.1, GL000217.1, GL000216.1, GL000215.1, GL000205.1, GL000219.1, GL000224.1, GL000223.1, GL000195.1, GL000212.1, GL000222.1, GL000200.1, GL000193.1, GL000194.1, GL000225.1, GL000192.1]
    ERROR ------------------------------------------------------------------------------------------

    Now I'm trying to reorder my bam with picard RerderSam.jar. I hope it will be the right way. What do you suggest?
    Many thanks,
    Francesco

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hmm, I'm not sure why the sorting order would be different if you mapped your reads to that reference, but using ReorderSAM should fix it, yes.

  • KathKath Member

    Can anyone tell me whether the human_g1k_b37 reference is masked? I can't seem to find out from 1000 genomes...

    Kath

  • JackJack Member

    Hi, there. I'm using SelectVariants to get the concordance of the variants got from GATK UnifiedGenotyper and SAMtools mpileup, but GATK kept complaining "Invalid command line: No tribble type was provided on the command line and the type of the file could not be determined dynamically. Please add an explicit type tag :NAME listing the correct type from among the supported types:

    ERROR Name FeatureType Documentation
    ERROR BCF2 VariantContext http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_variant_bcf2_BCF2Codec.html
    ERROR VCF VariantContext http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_variant_vcf_VCFCodec.html
    ERROR VCF3 VariantContext http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_variant_vcf_VCF3Codec.html"

    I have rerun UnifiedGenotyper and use the newly generated vcf file, but got the same error, I don't know why, can anyone give some ideas?

  • JackJack Member

    I have validate my vcf file, nothing seems wrong, I wonder how this happens ?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Jack, can you post your version and full command-line please?

  • JackJack Member

    I used GATK version 2.4-9. Here is my command line:
    for generating variant file,
    java -jar /home/Software/GenomeAnalysisTK-2.4-9/GenomeAnalysisTK.jar -R ref.fa -T UnifiedGenotyper -I sample.recal.bam -o INDEL.vcf --dbsnp dbSNP.vcf -glm INDEL -rf BadCigar
    and for SelectVariants,
    java -jar /home/Software/GenomeAnalysisTK-2.4-9/GenomeAnalysisTK.jar -R ref.fa -T SelectVariants --variant sample.gatk.raw.vcf --concordance sample.samtools.raw.vcf -o sample.concordance.raw.vcf
    I just don't know what's wrong...

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Jack,

    I'm not sure, the GATK should recognize automatically that your files are VCF. You can try to force recognition by writing --variant:VCF sample.gatk.raw.vcf for example.

    Are you using an older version of GATK? Can you tell from the error message exactly which file it is complaining about?

  • pdexheimerpdexheimer Member ✭✭✭✭

    I would concentrate on the samtools VCF. This is just a hypothesis, but it wouldn't surprise me greatly if it were a valid VCF version 4 file, and GATK requires 4.1 - or some similar situation

  • JackJack Member
    edited August 2013

    I would concentrate on the samtools VCF.

    I think I find the problem, it's the problem of samtools, now I have fixed it after rerunning, with a few arguments modified.

    Thank you, Geraldine and pdexheimer!

    Post edited by Geraldine_VdAuwera on
  • DavidRiesDavidRies Member ✭✭

    If anyone is still interested in it, I had the same problem and tried different aproaches to fix this:

    @jlin said:
    Geraldine,

    Your help would be appreciated! I can easily substitute characters with sed, but I'm not sure what characters I'm looking for. I searched for '10,' but the only 10s that come up are for chromosome number (and bp position in the index file).

    The older version of GATK did actually run on my reference file... so maybe I'm interpreting this error wrong?

    How long would it take to revert bams and remap them to the ref file from the GATK bundle?

    The full error message was as follows:

    ERROR ------------------------------------------------------------------------------------------
    ERROR A USER ERROR has occurred (version 2.6-2-ge03a5e9):
    ERROR The invalid arguments or inputs must be corrected before the GATK can proceed
    ERROR Please do not post this error to the GATK forum
    ERROR
    ERROR See the documentation (rerun with -h) for this tool to view allowable command-line arguments.
    ERROR Visit our website and forum for extensive documentation and answers to
    ERROR commonly asked questions http://www.broadinstitute.org/gatk
    ERROR
    ERROR MESSAGE: Bad input: We encountered a non-standard non-IUPAC base in the provided reference: '10'
    ERROR ------------------------------------------------------------------------------------------

    The command was: java -jar GenomeAnalysisTK.jar -R hg19.genome.new.fa -T UnifiedGenotyper -I bamlist.list -L chr1.bed -o chr1.vcf -U

    Finally, I just imported the reference in byopython as fasta, and then exported it again, which fixed the problem.

    `#!/usr/bin/env python

    from Bio import SeqIO

    handle = open("reference.fa", "rU")

    input_seq_iterator = SeqIO.parse(handle,"fasta")

    output_handle = open("reference_clean.fa", "w")

    SeqIO.write(input_seq_iterator, output_handle, "fasta")

    output_handle.close()
    `

    Best,

    David

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Thanks for reporting your solution, David -- and three cheers for Biopython!

  • longyanglongyang qingdaoMember

    Hi Geraldine!
    I'm new here to study gatk, and I'm very appreciated for your work. I have a question that I want to do snp calling for Marsupenaeus japonicus, and I have some NGS sequences, but I can't find appropriate reference sequences in the hige, then what should I do for this step? May I use the sequences in my hand itself as reference sequence? If not, is gatk the right software for me? Thanks a lot for your patient and labor, have a good day!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @longyang,

    Assuming you have several genomes from different individuals, you can do de novo assembly of one of them and use that as a reference to call variants on the others with GATK. Good luck!

  • longyanglongyang qingdaoMember
    edited December 2013

    Hi @Geraldine_VdAuwera

    So happy to receive your reply! Thanks for your advise, and I'll follow it. Have a good vacation!

  • mfletchermfletcher DEMember

    Hi @Geraldine_VdAuwera‌ (or whoever else is around!),

    Just a quick question about ROD formats for known SNP sites for VQSR. I'm trying a set of SNPchip data and the genotype calls are in the following csv format:

    JAX00000001,A,T,chr1,3013441,T,T,T,T,T,T,A,A,A,A,T,T,T,A,A,T,T,T,T,T,T,A,A,A,A,A,A,A,A,A,A,T,T,T,A,A,A,T,A,T,T,T,T,T,A,A,A,A,T,T,T,A,A,T,T,N,A,T,A,A,A,T,A,A,T,T,T,T,T,A,T,A,A,A,A,T,T,A,A,A,A,A,T,A,A,A,T,T,T,T,T,T,T,A,T,A,A,A,T,T,T,T,T,T,T,N,T,T,T,T,T,T,N,T,A,A,N,A,N,T,T,A,T,N,T,A,T,T,T,T,N,N,T,T,T,A,T,A,T,T,N,N,T,T,T,A,T,T,A,T,A,N,T,T,T,T,T,T,T,T,A,T,T,T,T,T,T,T,T,T,T,T,T,T,T,T,T,T,T,T,T,T,N,T,A,N,N,T,T,A,T,T,T,T,T,T,T,T
    Where col1=JAX SNP ID, col2= reference allele, col3 = alt allele, col4 = chrom, col5 = base position, and the following ~100 columns are per-mouse strain genotype calls.

    My question is: which ROD format should I use?

    My initial instinct would be to try and mangle the data into vcf format, but vcfs contain a lot of additional quality data that I just don't have.

    So I guess I should be using the "UCSC formatted dbSNP" instead, but I just wanted to check that this refers to the format described on the UCSC wiki here: http://genomewiki.ucsc.edu/index.php/DbSNP_Track_Notes#Subset_of_NCBI_fields_used_to_build_snpNNN_track

    Finally, do all of those UCSC dbSNP columns need to be populated, or will VariantRecalibrator be happy with what I have (the columns up to and including "observed")?

    Thanks very much!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @mfletcher, I'm not sure as I've never used UCSC dbsnp myself, but this loooks like it should work. If it doesn't, VQSR will let you know quickly (and in more detail than you really want) that it's not happy. If so let me know and we'll try to figure it out.

  • pdexheimerpdexheimer Member ✭✭✭✭

    @mfletcher - UCSC dbSNP format is a "sites-only" format, you won't be able to get the per-strain genotypes into it. Of course, VQSR doesn't care about the calls, so that's alright.

    It's pretty easy to make a minimal VCF, though - the rule of thumb is that any data you don't have can be replaced with a ".". So the VCF line for your example might be

    chr1  3013441  JAX00000001  A  T  .  PASS  .
    

    This is a sites-only record, and has no value for QUAL or for the INFO annotations. I would recommend converting to VCF rather than the UCSC record because support for VCF format will probably be much more stable over time

  • mfletchermfletcher DEMember
    edited June 2014

    Hi @Geraldine_VdAuwera‌, yes you're totally correct - as always the error message had way more detail than I really wanted...!

    `

    ERROR ------------------------------------------------------------------------------------------
    ERROR A USER ERROR has occurred (version 3.1-1-g07a4bf8):
    ERROR
    ERROR This means that one or more arguments or inputs in your command are incorrect.
    ERROR The error message below tells you what is the problem.
    ERROR
    ERROR If the problem is an invalid argument, please check the online documentation guide
    ERROR (or rerun your command with --help) to view allowable command-line arguments for this tool.
    ERROR
    ERROR Visit our website and forum for extensive documentation and answers to
    ERROR commonly asked questions http://www.broadinstitute.org/gatk
    ERROR
    ERROR Please do NOT post this error to the GATK forum unless you have really tried to fix it yourself.
    ERROR
    ERROR MESSAGE: Invalid command line: No tribble type was provided on the command line and the type of the file could not be determined dynamically. Please add an explicit type tag :NAME listing the correct type from among the supported types:
    ERROR Name FeatureType Documentation
    ERROR BCF2 VariantContext (this is an external codec and is not documented within GATK)
    ERROR VCF VariantContext (this is an external codec and is not documented within GATK)
    ERROR VCF3 VariantContext (this is an external codec and is not documented within GATK)
    ERROR ------------------------------------------------------------------------------------------

    `

    Interestingly there's no UCSC dbSNP option there... Suggestions? (Two comments: I commented out the header line of this file with a #, and saved it as a tab-delimited *txt file)

    Thanks again!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Ah, you're going to need to convert your file upfront, as it seems the VQSR isn't capable of reading in non-VCF/BCF formats natively. Have a look at the VariantsToVCF tool: http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_sting_gatk_walkers_variantutils_VariantsToVCF.html

  • asgroasgro RomaMember

    Hi @Geraldine_VdAuwera‌ ,
    Sorry but I don't understand very well this "Warning: this file format is 0-based for the start coordinates, so coordinates taken from 1-based formats should be offset by 1.". It means that I have to add (+1) to the start column possitions of my bed file or (-1) to the end column ?

    Thank you!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @asgro,

    It depends; if you are generating the bed file yourself based on a list of intervals that is in a different format, you need to make sure the intervals are adjusted appropriately. But if the bed file was pre-existing (e.g. if it's the intervals file supplied with the sequencing experiment) then you don't need to worry about it.

  • wangdongliwangdongli Member

    hi!@Geraldine_VdAuwera ,when I use gatk to call SNP.I met a problem puzzled me for some times, I hope you can help me ,the error message is:Lexicographically sorted human genome sequence detected in reads after I run the command:java -Xmx4g -jar /home/wangdongli/GATK/GenomeAnalysisTK.jar -R hg.19.fa -T RealignerTargetCreator -o my.realn.intervals -I my.bam.then I tried the Resortsam tool to resort my.bam file,and I saw the read as follows:
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchr8146364022
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchr8_gl000196_random38914
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchr8_gl000197_random37175
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchr9141213431
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchr9_gl000198_random90085
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchr9_gl000199_random169874
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchr9_gl000200_random187035
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchr9_gl000201_random36148
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchrM16571
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchrUn_gl000211166566
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchrUn_gl000212186858
    '''
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchrX155270560
    INFO 2017-03-09 09:41:01 ReorderSam SN=%s LN=%d%nchrY59373566
    so can you see where is the problem? and I use the hg19.fa as refererence,should I also resort the hg19.fa file? and how can I do it?

  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @wangdongli
    Hi,

    Can you please post your BAM file header and FASTA dict file?

    Thanks,
    Sheila

Sign In or Register to comment.