[Heisman]

I am not affiliated with VAAST in any way, but I have used it extensively and absolutely love it. I don't want to derail this thread in any way but I can certainly answer questions about it.

[ulz_peter]

thanks a lot ulz_peter!
Could you please clarify which steps of your pipeline are specifically for diploid genomes in order I can customize for my purposes?

I didn't find it yet, but there was a statement on the GATK homepage that the options descirbed there (which are basically pretty mcuh the same as I use) only work for diploid genomes and expected shifts of allele frequency must be adressed. So the question is: what are you planning to do: find rare alleles within some strains, sequence a genetically homogeneous strain...

[liu_xt005]

Problem with SNP-calling

Following ulz_peter's original doc, I have some problem when doing the SNP-calling.

java -Xmx4g -jar /path/GenomeAnalysisTK-1.1-35-ge253f6f/GenomeAnalysisTK.jar \
-glm BOTH \
-R hg18.fa \
-T UnifiedGenotyper \
-I myinput.marked.realigned.fixed.recal.bam \
-D dbsnp132_hg18.txt \
-o myoutput.snps.vcf \
-metrics snps.metrics \
-stand_call_conf 50.0 \
-stand_emit_conf 10.0 \
-dcov 1000 \
-A DepthOfCoverage \
-A AlleleBalance \
-L hg18_exonIntervals.bed

This "-L" option does not work.
I got the hg18_exonIntervals.bed from UCSC as ulz_peter's original doc shows.
I run the SNP-calling without the "-L" line.
Then the variant quality score recalibration step does not work, generating an empty output.tranches file.

Can somebody help me out? Thanks a lot.

[ulz_peter]

What is the error message when you specify the -L argument?

I actually stopped using the Variant quality Score recalibration as it often did not work out for me (I never work on more than 2 exomes at a time).

I out the version withouth the recalibration on the SEQanswers Wiki/How-To section. You may have a look there, as I will update that in the future and stop uploading newer versions of the PDF file...

[raonyguimaraes]

Following ulz_peter's original doc, I have some problem when doing the SNP-calling.

java -Xmx4g -jar /path/GenomeAnalysisTK-1.1-35-ge253f6f/GenomeAnalysisTK.jar \
-glm BOTH \
-R hg18.fa \
-T UnifiedGenotyper \
-I myinput.marked.realigned.fixed.recal.bam \
-D dbsnp132_hg18.txt \
-o myoutput.snps.vcf \
-metrics snps.metrics \
-stand_call_conf 50.0 \
-stand_emit_conf 10.0 \
-dcov 1000 \
-A DepthOfCoverage \
-A AlleleBalance \
-L hg18_exonIntervals.bed

This "-L" option does not work.
I got the hg18_exonIntervals.bed from UCSC as ulz_peter's original doc shows.
I run the SNP-calling without the "-L" line.
Then the variant quality score recalibration step does not work, generating an empty output.tranches file.

Can somebody help me out? Thanks a lot.

Did you try with -B:targetIntervals,BED hg18_exonIntervals.bed?

By the way I couldn't figure out how to use this on version 1.2

[mirabilia]

I didn't find it yet, but there was a statement on the GATK homepage that the options descirbed there (which are basically pretty mcuh the same as I use) only work for diploid genomes and expected shifts of allele frequency must be adressed. So the question is: what are you planning to do: find rare alleles within some strains, sequence a genetically homogeneous strain...

The first part of the analysis, until the local realignment around indels, is of course suitable for all genomes I have. Quality score recalibration and SNP calling need customization and I'm thinking on... Actually, I'd like to catch rare alleles in sequenced bacterial populations, like repeats, tandem repeats and structural variants, but I suspect it's a really difficult task to discriminate between real rare variations and noise introduced by such data.
Obviously any kind of suggestion, it's really appreciate!

[liu_xt005]

-L option problem solved

What is the error message when you specify the -L argument?

I actually stopped using the Variant quality Score recalibration as it often did not work out for me (I never work on more than 2 exomes at a time).

I out the version withouth the recalibration on the SEQanswers Wiki/How-To section. You may have a look there, as I will update that in the future and stop uploading newer versions of the PDF file...

Dear ulz_peter and raonyguimaraes,

Thanks VERY MUCH to both of you.
The problem seems to be solved by removing the random and hap intervals from the .bed file.

ulz_peter,
raonyguimaraes posted a similar pipeline by Gayle Philip.
SNPs and Indels are recalibrated/filtered separately and combined after.
I am trying that, and think that it is a good idea to exclude Indels from Gaussian models.

[MQ-BCBB]

Great resource!!!
Just a few comments:
In the picard/MarkDuplicates.jar, the option 'CREATE_INDEX=true' should be added.
With respect to adding read group information, instead of using the bwa sampe -r option, picard AddOrReplaceReadGroups.jar is an easier way to go as it tells you which options are required. Thanks for sharing!

[ulz_peter]

You're absolutely right on the first one.
For the sampe -r option: I'd like to keep it in the bwa part as another step handling the BAM file takes a lot of time than just doing it once... but I could add that as an alternative

[emilyjia2000]

reference dictionary

When I tried to use GATK to do the local realignment according to ulz_peter's instruction, one error message occurred: Invalid command line: Failed to load reference dictionary. Could anybody let me know where to get this reference dictionary? how to use it in the command line?

Thanks in advance.

[raonyguimaraes]

Post your command line here... and the full error message

[Heisman]

When I tried to use GATK to do the local realignment according to ulz_peter's instruction, one error message occurred: Invalid command line: Failed to load reference dictionary. Could anybody let me know where to get this reference dictionary? how to use it in the command line?

Thanks in advance.

Look here: http://www.broadinstitute.org/gsa/wi...ference_genome

[liu_xt005]

A quick comment: I would suggest people to do the SAM -> BAM conversion using Picard.
In fact, samtools sort generates 79 temporary files for one sample.

[emilyjia2000]

Thanks for all of your quick response. I used the command line:

java -Xmx4g -jar /path/to/GenomeAnalysisTK.jar -T RealignerTargetCreator -R hg19.fa -o output.interval -I /path/to/reorder_dedup.bam

I already copied the ucsc.hg19.dict in the same directory.

When I run this command, the error message:

##### ERROR ------------------------------------------------------------------------------------------
##### ERROR A USER ERROR has occurred (version 1.2-26-g43b0c98):
##### 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 wiki for extensive documentation http://www.broadinstitute.org/gsa/wiki
##### ERROR Visit our forum to view answers to commonly asked questions http://getsatisfaction.com/gsa
##### ERROR
##### ERROR MESSAGE: Invalid command line: Failed to load reference dictionary
##### ERROR ------------------------------------------------------------------------------------------

[raonyguimaraes]

If I'm not wrong ucsc.hg19.dict should be on the same directory of the reference file.