SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Ion Torrent $1000 Genome!? Benchtop Ion Proton Sequencer aeonsim Ion Torrent 88 10-28-2012 04:50 AM
Ion Torrent data quality impressions? rcorbett Ion Torrent 15 03-27-2012 11:43 AM
GATK and Ion Torrent data NextGenSeb Bioinformatics 2 02-12-2012 08:36 PM
Technical Brief: Processing Ion Torrent Data with RTG Investigator Stewart Noyce Vendor Forum 0 10-19-2011 12:08 PM
Ion Torrent PGM data with Mosaik-aligner Magnus Bioinformatics 0 05-19-2011 03:59 AM

Reply
 
Thread Tools
Old 05-03-2012, 11:04 AM   #1
aituka
Member
 
Location: brussels

Join Date: Mar 2012
Posts: 13
Default pipeline for ion torrent data

Hi,
I just received one ion torrent data to test (sample.fastq). I wanted to have the community's point of view on what I am doing.

1/ aligning the reads
-------------------------
I use:
$bwa -aln hg19ref.fa sample.fastq > sample.sai
$bwa samse -r "@RG\tID:IDa\tSM:sample\tPL:PL" hg19ref.fa sample.sai sample.fastq > sample.sam
$java -Xmx4g -jar picard/SortSam.jar SO=coordinate INPUT=sample.sam OUTPUT=sample.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true

looking to flagstat, I get:
428427 in total
0 QC failure
0 duplicates
373491 mapped (87.18%)

RMK/questions:
* should I rather use bwa -bwasw ?? or bowtie2??
* is the platform header information included in the sam used somewhere?

2/ gatk1.3 guideline: realigning reads, etc.
-----------------------------------------------

ooooo Marking PCR duplicates BAM ooooo
java -Xmx4g -jar picard/MarkDuplicates.jar INPUT=sample.bam OUTPUT=sample.marked.bam METRICS_FILE=metrics VALIDATION_STRINGENCY=LENIENT

ooooo Indexing MARKED BAM ooooooooo '
samtools index sample.marked.bam

ooooo Realignment around indels (1) create Indels Table ooooo
java -Xmx4g -jar gatk.jar -T RealignerTargetCreator -R hg19ref.fa -filterMBQ -nt 16 --known GATK_1000G_PHASE1_INDELS --known GATK_MILLS1000G_GOLD_INDELS -o indeltablelist -I sample.marked.bam

ooooo Realignment around indels (2) realigns reads around those targets ooooo
java -Xmx4g -jar gatk.jar -I sample.marked.bam -R hg19ref.fa -filterMBQ -known GATK_1000G_PHASE1_INDELS -known GATK_MILLS1000G_GOLD_INDELS -T IndelRealigner -targetIntervals indeltablelist -o sample.marked.realigned.bam

ooooo Fix mate for pair end reads oooooooo '
java -Xmx4g -jar picard/FixMateInformation.jar INPUT=sample.marked.realigned.bam
OUTPUT=sample.marked.realigned.fixed.bam SO=coordinate VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true

ooooo Quality score recalibration (1) Countcovariates oooooooo
java -Xmx4g -jar gatk.jar -l INFO -R hg19ref.fa -nt 16 -knownSites:dbsnp,VCF dbsnp_135.b37.vcf -I sample.marked.realigned.fixed.bam -T CountCovariates -cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate -cov DinucCovariate -recalFile sample.recal_data.csv


RMK: This was failing with gatk 1.5 (complaining that the platform was unknown - wanted illumina or 454 or solid). It seems it works with gatk1.6

ooooo Quality score recalibration (2) Tablerecalibration oooooooo '
java -Xmx4g -jar gatk.jar -l INFO -R hg19ref.fa -I sample.marked.realigned.fixed.bam -T TableRecalibration --out sample.marked.realigned.fixed.recal.bam -recalFile sample.recal_data.csv


Does it looks ok?
Then calling the variant using the rest of the pipeline guideline.

any comments?
thanks
tuka
aituka is offline   Reply With Quote
Old 05-09-2012, 10:36 PM   #2
kenietz
Member
 
Location: Singapore

Join Date: Nov 2011
Posts: 85
Default

hi aituka,
i was trying to figure out how to use GATK with torrent data as well. More or less i ended up with something similar to your pipeline. But my customers gave me some info about expected deletions which i could not find using that pipeline. Turned out to be a problem with markduplicate step. After the marking lots of SNP and Indels disappear from the vcf even though one can see the expected deletion in the bam file when using IGV.

So for now after i add readgroup, sort and index the bam file i directly go to UnifiedGenotyper and VariantFiltration from GATK. So i skip all recalibration and etc. Interestingly the results with and without recalibration,both excluding markduplicates are almost identical.

Btw, i use 'bwa bwasw' to align my data to human.
kenietz is offline   Reply With Quote
Old 05-10-2012, 12:20 AM   #3
colinmolter
Member
 
Location: belgium

Join Date: Jul 2011
Posts: 11
Default

Hi tuka,
I do agree with kenietz. The markduplicate trick might remove too much reads. The problem relates to this post: http://seqanswers.com/forums/showthread.php?t=6854 (Removing duplicates is it really necessary?). When sequencing the whole exome, you don't expect too much duplicates. Removing them might be a good option to avoid PCR duplicates. However, in targeted sequencing, it might be natural to have a lot of duplicates. Removing them might not be the good option.

In an example of ion PGM data, I got:
429,424 reads, among which I could align (bwa bwasw) 383,091 (89.21%) reads.
But there was 371,738 duplicates !

As for Kenietz, results with or without duplicates were identical.
colinmolter is offline   Reply With Quote
Old 05-10-2012, 12:36 AM   #4
kenietz
Member
 
Location: Singapore

Join Date: Nov 2011
Posts: 85
Default

@colinmolter:
exactly the same situation. i have this targeted reseq with coverage of 1000-2000x and around 95% of the reads were marked duplicates.
Its sad that the documentation on GATK is not enough. But unfortunately that is the situation one can not describe all possibilities and hence sometimes when one follows the general guide one is getting lost. Is good that there are many forums around tho

Cheers
kenietz is offline   Reply With Quote
Old 05-17-2012, 10:53 PM   #5
kenietz
Member
 
Location: Singapore

Join Date: Nov 2011
Posts: 85
Default

Hi guys,
what argument did you use for UnifiedGenotyper and VariantFiltration?
kenietz is offline   Reply With Quote
Old 05-17-2012, 11:46 PM   #6
colinmolter
Member
 
Location: belgium

Join Date: Jul 2011
Posts: 11
Default

Quote:
Originally Posted by kenietz View Post
Hi guys,
what argument did you use for UnifiedGenotyper and VariantFiltration?
I used something like this:

Code:
java -Xmx4g -jar gatk.jar -T UnifiedGenotyper -R bwa6.1/h_g1k_v37.fasta -L targetintervals.bed -nt 16 -A AlleleBalance -A DepthOfCoverage -stand_call_conf 30.0 -stand_emit_conf 10 -glm BOTH --dbsnp dbsnp_135.b37.vcf -o DS.SNV.all.vcf -metrics DS.SNV.all.vcf.metric -I s1.bam -I s2.bam
what about yours?
any comments?
colin
colinmolter is offline   Reply With Quote
Old 05-17-2012, 11:55 PM   #7
kenietz
Member
 
Location: Singapore

Join Date: Nov 2011
Posts: 85
Default

Hi colin,
i followed a guide for GATK for Illumina which i found on this site. But had to modify it just a tiny bit.

`java -jar /opt/GATK/GenomeAnalysisTK.jar -R $refgen -T UnifiedGenotyper -I realigned.recal.bam --dbsnp /mnt/hd/GATK_recource_bundle/hg19/dbsnp_135.hg19.vcf.reordered -o gatk_var.raw.vcf --num_threads 8 -L pbrm1_intervals.bed --genotype_likelihoods_model BOTH --metrics_file snp.metrics -stand_call_conf 30 -stand_emit_conf 10 -A DepthOfCoverage -A AlleleBalance`;

`java -Xmx4g -jar /opt/GATK/GenomeAnalysisTK.jar -R $refgen -T VariantFiltration --variant gatk_var.raw.vcf -o gatk_var_filtered.vcf --clusterWindowSize 10 --filterExpression "MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1)" --filterName "HARD_TO_VALIDATE" --filterExpression "DP < 5 " --filterName "LowCoverage" --filterExpression "QUAL < 30.0 " --filterName "VeryLowQual" --filterExpression "QUAL > 30.0 && QUAL < 50.0 " --filterName "LowQual" --filterExpression "QD < 1.5 " --filterName "LowQD" --filterExpression "SB > -10.0 " --filterName "StrandBias"`;

Seems that results from GATK and VariantCaller plugin from torrent suite differ
Now trying to make the standalone variantcaller to run. Partial success for now tho
kenietz is offline   Reply With Quote
Old 05-18-2012, 10:41 AM   #8
adaptivegenome
Super Moderator
 
Location: US

Join Date: Nov 2009
Posts: 437
Default

Has anyone tried TMAP over BWA?
adaptivegenome is offline   Reply With Quote
Old 05-18-2012, 12:15 PM   #9
aggp11
Member
 
Location: Wisconsin

Join Date: Jun 2011
Posts: 87
Default

@genericforms

I am sorry, but I am not sure if your question is more general or specifically towards this problem.

In case it is more general, then yes I have tried both TMAP and BWA on our datasets and it seems like there is not much of a difference if you select the right BWA algorithm (aln-samse / bwasw).. As TMAP selects the appropriate algorithm based on the characteristics (most probably read length).. unless off-course TMAP decides to run SSAHA, which would give different results.

Thanks,
Praful
aggp11 is offline   Reply With Quote
Old 05-18-2012, 12:17 PM   #10
adaptivegenome
Super Moderator
 
Location: US

Join Date: Nov 2009
Posts: 437
Default

Quote:
Originally Posted by aggp11 View Post
@genericforms

I am sorry, but I am not sure if your question is more general or specifically towards this problem.

In case it is more general, then yes I have tried both TMAP and BWA on our datasets and it seems like there is not much of a difference if you select the right BWA algorithm (aln-samse / bwasw).. As TMAP selects the appropriate algorithm based on the characteristics (most probably read length).. unless off-course TMAP decides to run SSAHA, which would give different results.

Thanks,
Praful
Yes, I was curious in your experience how BWA-SW matches up to TMAP...
adaptivegenome is offline   Reply With Quote
Old 05-20-2012, 05:40 PM   #11
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Quote:
Originally Posted by aggp11 View Post
@genericforms

I am sorry, but I am not sure if your question is more general or specifically towards this problem.

In case it is more general, then yes I have tried both TMAP and BWA on our datasets and it seems like there is not much of a difference if you select the right BWA algorithm (aln-samse / bwasw).. As TMAP selects the appropriate algorithm based on the characteristics (most probably read length).. unless off-course TMAP decides to run SSAHA, which would give different results.

Thanks,
Praful

Simulate some Ion 100-200bp human data and compare the two software for yourself (included are the simulator, simulation results evaluator, and plotting software
. Feel free to create a new post to give feedback
nilshomer is offline   Reply With Quote
Old 05-21-2012, 04:08 PM   #12
adaptivegenome
Super Moderator
 
Location: US

Join Date: Nov 2009
Posts: 437
Default

We are working on it Nils. In a moment of weakness I got lazy and wanted to know if someone had already looked into this...
adaptivegenome is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 09:06 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO