PCR duplicates, Rna-Seq and GATK

[choishingwan]

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?

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

[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

[whataBamBam]

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

[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?

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

[dpryan]

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

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

[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 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]

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!

[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]

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]

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

[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?