I am working on a project which requires me to identify Variants in a marine species using RNAseq data.
I have assembled the transcriptome (using all available RNASeq data sets for the experiment) for my organism (no Reference is available).
I then referred to the link 'http://gatkforums.broadinstitute.org/gatk/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference' and created the '.dict' file and indexed the reference.
java -Xmx3g -jar /share/apps/picard/1.115/CreateSequenceDictionary.jar R=seaurchin_topHit_Ids.fasta O=seaurchin_topHit_Ids.dict samtools faidx seaurchin_topHit_Ids.fasta
I have the following RNASeq data-sets available: Two growth Conditions : Control (C) vs Bubble (B) At two Growth Times: 8h and 48h
RNASeq data is generated in replicates (R1 and R2)
So in all I have the following read-sets C1-8h , C2-8h C1-48h, C2-48h
B1-8h, B2-8h B1-48h, B2-48h
The questions we might want to ask are : Identify variants : Q1)Across conditions : Condition C (Control) vs Condition B Q2)Across times : Growth times 8h vs 48h
I am following the best practice guidelines for GATK for RNASeq data..
1) Mapping individual data-sets using bwa (I will repeat the pipeline using STAR later, if I can get through all steps)..I have been using bwa for a while, so felt more comfortable using the tool.
2) sort the individual 'sam files' and convert into 'bam'
3) Mark Duplicates for individual bams to get the '.dedup.bam files'
4) Add Read groups I am keeping the RGSM tag common across datasets belonging to the same sample and separating the replicates by the RGID tag.. e.g. java -Xmx3g -jar /share/apps/picard/1.115/AddOrReplaceReadGroups.jar \ I= B1-BB-8h.dedup.bam \ O= B1-BB-8h.dedup_RG.bam \ SORT_ORDER=coordinate \ RGID=B1-BB-8h.dedup \ RGLB=BB8h \ RGPL=illumina \ RGSM=BB8h \ CREATE_INDEX=True \ RGPU=run_barcode
java -Xmx3g -jar /share/apps/picard/1.115/AddOrReplaceReadGroups.jar \ I= B2-BB-8h.dedup.bam \ O= B2-BB-8h.dedup_RG.bam \ SORT_ORDER=coordinate \ RGID=B2-BB-8h.dedup \ RGLB=BB8h \ RGPL=illumina \ RGSM=BB8h \ CREATE_INDEX=True \ RGPU=run_barcode
java -Xmx3g -jar /share/apps/picard/1.115/AddOrReplaceReadGroups.jar \ I= B1-BC-8h.dedup.bam \ O= B1-BC-8h.dedup_RG.bam \ SORT_ORDER=coordinate \ RGID=B1-BC-8h.dedup \ RGLB=BB8h \ RGPL=illumina \ RGSM=BC8h \ CREATE_INDEX=True \ RGPU=run_barcode
Can someone confirm if the above steps look alright?
I am a bit unsure as to what I do further!!
Should I merge all dedup.bam files and do the downstream processing ? Since I have the read groups marked (associating samples but separating the conditions and growth times), can I expect to identify SNPs/indels etc for Q1 and Q2 ?
I will appreciate if someone can guide me through this,
I have DNASeq data of an offspring "C" which is cross of parents "A" and "B". Reference genomes of both "A" and "B" are available. I am interested in getting the variants, which should be in the "A" homozygous region. I am confused which reference to use and how to proceed. Your advice will be very valuable. Thank you.
Hi there, According to your experience, what do you suggest for calling variants in tumor samples (exome) - HaplotypeCaller or UnifiedGenotyper? Have you seen any (major or minor) difference with both the callers on same exome tumor samples? What's your advice. Do you have any reference publication(s) where HaplotypeCaller is used on tumor samples? Thanks Raj
I want to know what's the best way to use VariantEval to get statistics for each sample in a multisample VCF file. If I call it like this:
java -jar GenomeAnalysisTK.jar \
-R ucsc.hg19.fasta \
-T VariantEval \
-o multisample.eval.gatkreport \
--eval annotated.combined.vcf.gz \
where annotated.combined.vcf.gz is a VCF file that contains ~1Mio variants for ~800 samples I get statistics for all samples combined, e.g.
#:GATKTable:CompOverlap:The overlap between eval and comp sites
CompOverlap CompRod EvalRod JexlExpression Novelty nEvalVariants ...
CompOverlap dbsnp eval none all 471704 191147
CompOverlap dbsnp eval none known 280557 0 CompOverlap dbsnp eval none novel 191147 191147
But I would like to get one such entry per sample. Is there an easy way to do this?