PCR duplicates, Rna-Seq and GATK

[whataBamBam]

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).

[kmcarr]

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:

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?

[whataBamBam]

Hi kmcarr

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.

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.

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]

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

[Kennels]

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).

[whataBamBam]

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!

[vishnuamaram]

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.

[whataBamBam]

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

[vishnuamaram]

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.

[dpryan]

Usually you would use RealignerTargetCreator to create an interval file for use by IndelRealigner.

[sindrle]

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.

[whataBamBam]

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]

sindrle no I didn't.. this thread satisfied me that marking dupes is not the correct thing to do for Rna-Seq

[sindrle]

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%3Adoi%2F10.1371%2Fjournal.pone.0058815

[bruce01]

Good thread, this has got me thinking. Just wondering where you get this from:

...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.