SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Ion torrent PCR duplicates DavyK Bioinformatics 0 08-29-2013 03:14 AM
Ignore PCR Duplicates in CollectTargetedPcrMetrics fongchun Bioinformatics 5 02-16-2013 03:51 PM
PCR duplicates aquleaf RNA Sequencing 5 10-30-2012 11:40 AM
PCR duplicates questions slny Bioinformatics 8 06-07-2011 04:06 AM
how critical is the filtering of potential PCR duplicates? julien General 3 03-26-2010 10:24 AM

Reply
 
Thread Tools
Old 09-11-2013, 12:06 AM   #1
whataBamBam
Member
 
Location: Italy

Join Date: May 2013
Posts: 27
Default PCR duplicates, Rna-Seq and GATK

Hi all..

I gather from this thread on biostars.

http://www.biostars.org/p/55648/

Linking these threads on biostars and here..

http://www.biostars.org/p/14283/

http://www.biostars.org/p/47229/

http://seqanswers.com/forums/showthread.php?t=6854

..and from some other reading around the forums.. that there's a general consensus NOT to remove PCR duplicates when doing Rna-Seq (at least at the stage of differential expression analysis) - for the perfectly sensible reasons discussed in those threads.

Previously with Dna resequencing and SNP calling, PCR duplicate removal was part of my standard preprocessing pipeline. Fine this is different - this is Rna we are dealing with now.

Now I am interested in doing some SNP calling. I'm planning to use Samtools and also GATK in order to achieve a consensus set of SNPs by taking the intersection of the results from the two tools (since this has worked well for me in the past).

Now my question is this..

I know from past experience that GATK's UnifiedGenotyper won't actually allow you to run it on Bam files which have not had PCR duplicates marked. So I'm wondering what my best option is here?

What I'm thinking is to use the alignments without duplicate removal for differential expression analysis.. and then use a second set which have had the dupes taken out for calling SNPs.

Is that a good idea? does that make better sense for the SNP calling regardless of GATK's fussiness over input files? It does seem to make sense to me that we want identically mapped reads for the differential expression analysis but not for the variant caling. However I don't want to follow an illogical procedure / do bad science just because of GATKs requirements on input files i.e. I want my decision to be based on science not getting tools to run - I could always just use a different SNP caller for consensus (i.e. Freebayes).
whataBamBam is offline   Reply With Quote
Old 09-11-2013, 09:38 AM   #2
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,178
Default

Duplicate removal is valid when removing PCR duplicates; you do not want to remove duplicate reads which arose from independent fragments. Duplicate marking/removal programs can not distinguish between these two types of duplicates. They have the built in assumption that any duplicate found is the result of PCR duplication. This is a reasonable assumption if your reads are from a genomic DNA library. It is NOT a valid assumption for RNA-Seq data. For RNA-Seq it is more likely that observed duplicates are from independent cDNAs from highly abundant transcripts. You still do not want to remove duplicates from RNA-Seq data even if you are doing SNP analysis with it.

I do have a question about this comment:

Quote:
I know from past experience that GATK's UnifiedGenotyper won't actually allow you to run it on Bam files which have not had PCR duplicates marked.
How does the UnifiedGenotyper "know" whether or not the input BAM file has had MarkDuplicates run on it? (I'm asking, I honestly don't know.) Unless it rechecks the input for duplicates the only way it would know is by finding reads with the duplicate flag set. Let's imagine you run MarkDuplicates on a BAM file which had no duplicates at all in it (this is thought experiment, just go with it). There would be no evidence recorded in the file that MarkDuplicates was run on that particular BAM. Would UnifiedGenotyper refuse to accept this file?
kmcarr is offline   Reply With Quote
Old 09-11-2013, 09:37 PM   #3
whataBamBam
Member
 
Location: Italy

Join Date: May 2013
Posts: 27
Default

Hi kmcarr

Quote:
Originally Posted by kmcarr View Post
Duplicate removal is valid when removing PCR duplicates; you do not want to remove duplicate reads which arose from independent fragments. Duplicate marking/removal programs can not distinguish between these two types of duplicates. They have the built in assumption that any duplicate found is the result of PCR duplication. This is a reasonable assumption if your reads are from a genomic DNA library. It is NOT a valid assumption for RNA-Seq data. For RNA-Seq it is more likely that observed duplicates are from independent cDNAs from highly abundant transcripts.
Great - that's what I thought from read around the forums. Yeah for sure the program has no idea - it's just blindly tagging reads that map in the same place.

Quote:
Originally Posted by kmcarr View Post
You still do not want to remove duplicates from RNA-Seq data even if you are doing SNP analysis with it.
Perfect - that's what I needed to know, thankyou.

Quote:
Originally Posted by kmcarr View Post
How does the UnifiedGenotyper "know" whether or not the input BAM file has had MarkDuplicates run on it? (I'm asking, I honestly don't know.) Unless it rechecks the input for duplicates the only way it would know is by finding reads with the duplicate flag set. Let's imagine you run MarkDuplicates on a BAM file which had no duplicates at all in it (this is thought experiment, just go with it). There would be no evidence recorded in the file that MarkDuplicates was run on that particular BAM. Would UnifiedGenotyper refuse to accept this file?
It just looks at the @PG tags in the Bam header - if you just add an @PG flag for a duplicate removal program it stops complaining. I think that's right.. I did play around with 'hacking' Bam files a bit for various things GATK moans about to get it to accept them (like I say it's fussy). Mostly just to work out exactly what it's doing. Pretty sure I did the test I mention here.. In any case I'll try that again now so I can let you know. If I want these reads even for SNP calling then that sounds the way to go. That was what I needed to know.. whether they are still needed for the SNPs so thanks.

Yeah it's just a bit annoying they have that limitation because messing around with the Bams like that isn't very publishable. For example for some previous work I changed some stuff in the Bam to make it accept alignments from BFAST (it wants you to use BWA). Fine for my own purposes but not very reportable. Still it's information I wouldn't have had otherwise. I only want this consensus SNP set for doing some parameter setting anyway. I search out values for some parameters by looking when the SNP sets are at maximum convergence. Any comments on whether that is a crazy fool method are welcome I'm experimenting here
whataBamBam is offline   Reply With Quote
Old 09-11-2013, 09:40 PM   #4
whataBamBam
Member
 
Location: Italy

Join Date: May 2013
Posts: 27
Default

Oh sorry - the last bit of your question

Yeah there would be evidence, because every program in the pipeline that works on a Bam adds an @PG flag... or in any case mostly they do (I believe they 'should' do / are supposed to)

So downstram programs know what upstream programs have touched the Bam even if they do nothing but add that flag
whataBamBam is offline   Reply With Quote
Old 10-01-2013, 11:16 PM   #5
Kennels
Senior Member
 
Location: Sydney

Join Date: Feb 2011
Posts: 149
Default

Hi
I am doing a similar analysis to you, and have run BWA-MEM aligned files, and GATK v2.2-3.

Are you sure that GATK Unified Genotyper tool doesn't run when duplicates are not marked? I could run it even when I do and don't use Picard MarkDuplicates on my input .bam file (i.e. duplicates marked or not).
Kennels is offline   Reply With Quote
Old 10-02-2013, 01:20 AM   #6
whataBamBam
Member
 
Location: Italy

Join Date: May 2013
Posts: 27
Default

Quote:
Originally Posted by Kennels View Post
Hi
I am doing a similar analysis to you, and have run BWA-MEM aligned files, and GATK v2.2-3.

Are you sure that GATK Unified Genotyper tool doesn't run when duplicates are not marked? I could run it even when I do and don't use Picard MarkDuplicates on my input .bam file (i.e. duplicates marked or not).
Hi Kennels

I'm quoting this from memory of my previous project..

- It could of been Picard tools that was refusing to process without dupe marking and I just remember wrong (as that was part of my pipeline)

or

- They could of changed the program since then.. I know they are updating that thing all the time.. I asked about the ploidy setting over on GATK forum and that had been added since the publication of the GATK paper

if it's working then it's working.. great

I haven't tried it myself yet this time around. Haven't got as far as SNP calling yet.. far too many other issues to worry about on this project for now!
whataBamBam is offline   Reply With Quote
Old 10-02-2013, 06:45 AM   #7
vishnuamaram
Member
 
Location: india

Join Date: Jun 2013
Posts: 42
Default

Hey guys, (kennels, whatabambam)

Kindly any of you respond and let me know,

what are the proper commands to be used for GATK- indel realigner, BQSR and variant calling.

When ever i run a GATK commands, it throws an error.

Kindly help me.

Thank you,
Vishnu.
vishnuamaram is offline   Reply With Quote
Old 10-02-2013, 04:56 PM   #8
whataBamBam
Member
 
Location: Italy

Join Date: May 2013
Posts: 27
Default

What are the errors your getting? Possibly this is the sort of stuff I was talking about.. it's fussy about input files. I've only used UnifiedGenotyper but I think a lot of the requirements are for all GATK walkers.

To be fair though it gives quite informative error messages.. if your just using it wrong it's probably pretty easy to figure out what's wrong

Last edited by whataBamBam; 10-02-2013 at 04:58 PM. Reason: forgot something
whataBamBam is offline   Reply With Quote
Old 10-02-2013, 11:20 PM   #9
vishnuamaram
Member
 
Location: india

Join Date: Jun 2013
Posts: 42
Default

Hey whatabambam,

this is the cd i used:-
java -Xmx30g -jar GenomeAnalysisTK-2.6-5-gba531bd/GenomeAnalysisTK.jar -T IndelRealigner \ -R REF_GENOME_hg19/hg19_karyotypic.fa -I WG_PBMC_SAM9_CORDSORT_RMDUp_KARYTP.BAM -targetIntervals
intervalListFromRTC.intervals \ -o WG_PBMC_SAM9_indelrealigned.bam

here is the error:-
##### ERROR MESSAGE: Couldn't read file /data/odity/Project_Blood-GNPC-464/Sample_WG-9/RAWFASTQ_CAT_FILES/intervalListFromRTC.intervals because The interval file does not exist

How exactly should we generate this interval file. Truly no idea. Let me know.

Thanking you in tons,
Vishnu.
vishnuamaram is offline   Reply With Quote
Old 10-03-2013, 04:33 AM   #10
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Usually you would use RealignerTargetCreator to create an interval file for use by IndelRealigner.
dpryan is offline   Reply With Quote
Old 10-03-2013, 01:16 PM   #11
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

Hi!
Did you finally mark duplicates (with Picard)?
Im doing the same thing as you, but I marked my duplicates.

I will maybe run it all again do compare the result.
sindrle is offline   Reply With Quote
Old 10-07-2013, 07:10 PM   #12
whataBamBam
Member
 
Location: Italy

Join Date: May 2013
Posts: 27
Default

Hi sorry vishnuamaram

I posted a reply to your question (or thought I did) but my browser crashed in posting it..

Did you fix this? It just sounds like your path was incorrect yes?
whataBamBam is offline   Reply With Quote
Old 10-07-2013, 07:11 PM   #13
whataBamBam
Member
 
Location: Italy

Join Date: May 2013
Posts: 27
Default

sindrle no I didn't.. this thread satisfied me that marking dupes is not the correct thing to do for Rna-Seq
whataBamBam is offline   Reply With Quote
Old 10-08-2013, 04:43 AM   #14
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

Ok, I did mark duplicates and my results looks good, but maybe Ill run everything one more time not marking duplicates to check if even better .

Have you read this? Its what I based on:

http://www.plosone.org/article/info%...l.pone.0058815
sindrle is offline   Reply With Quote
Old 10-08-2013, 06:55 AM   #15
bruce01
Senior Member
 
Location: .

Join Date: Mar 2011
Posts: 157
Default

Good thread, this has got me thinking. Just wondering where you get this from:
Quote:
Originally Posted by kmcarr View Post
...They [PCR duplicate removal software] have the built in assumption that any duplicate found is the result of PCR duplication. This is a reasonable assumption if your reads are from a genomic DNA library. It is NOT a valid assumption for RNA-Seq data. For RNA-Seq it is more likely that observed duplicates are from independent cDNAs from highly abundant transcripts. You still do not want to remove duplicates from RNA-Seq data even if you are doing SNP analysis with it.
Even if there are 90% duplicates from independent cDNA fragments how can you determine that any given duplicate is not a PCR duplicate? Essentially, if you are being conservative, is the correct approach not to assume all duplicates are from PCR and remove them (ie reduce type 1 error)? Or do you screen all potential SNPs using Sanger anyway, so get as many SNPs as possible and then validate them?

In a nutshell: for DE analysis the transcripts are 'functional' and so we should retain the duplicates; for SNP analysis we want a consensus of the structure of our sequence.
bruce01 is offline   Reply With Quote
Old 10-09-2013, 03:08 AM   #16
choishingwan
Member
 
Location: Hong Kong

Join Date: Feb 2012
Posts: 21
Default

I do have a set of data containing both Exome Seq and RNA Seq. It will be interesting to see how they correlate in terms of calling. Specifically, it will be interested to see how mark duplication will affect the concordance between the RNA samples from the DNA samples. However, I do have a question: is it valid to use tools like GATK to perform snp calling on RNA Samples? To my knowledge, it seems like GATK might contain certain prior specifically designed for exome sequencing (or whole genome sequencing). Wouldn't that also affect the call concordance in the final data?
choishingwan 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 06:59 PM.


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