Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

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

  • #2
    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?

    Comment


    • #3
      Hi kmcarr

      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.

      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.

      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

      Comment


      • #4
        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

        Comment


        • #5
          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).

          Comment


          • #6
            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!

            Comment


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

              Comment


              • #8
                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, 03:58 PM. Reason: forgot something

                Comment


                • #9
                  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.

                  Comment


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

                    Comment


                    • #11
                      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.

                      Comment


                      • #12
                        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?

                        Comment


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

                          Comment


                          • #14
                            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:

                            Next-generation RNA sequencing (RNA-seq) maps and analyzes transcriptomes and generates data on sequence variation in expressed genes. There are few reported studies on analysis strategies to maximize the yield of quality RNA-seq SNP data. We evaluated the performance of different SNP-calling methods following alignment to both genome and transcriptome by applying them to RNA-seq data from a HapMap lymphoblastoid cell line sample and comparing results with sequence variation data from 1000 Genomes. We determined that the best method to achieve high specificity and sensitivity, and greatest number of SNP calls, is to remove duplicate sequence reads after alignment to the genome and to call SNPs using SAMtools. The accuracy of SNP calls is dependent on sequence coverage available. In terms of specificity, 89% of RNA-seq SNPs calls were true variants where coverage is >10X. In terms of sensitivity, at >10X coverage 92% of all expected SNPs in expressed exons could be detected. Overall, the results indicate that RNA-seq SNP data are a very useful by-product of sequence-based transcriptome analysis. If RNA-seq is applied to disease tissue samples and assuming that genes carrying mutations relevant to disease biology are being expressed, a very high proportion of these mutations can be detected.

                            Comment


                            • #15
                              Good thread, this has got me thinking. Just wondering where you get this from:
                              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.

                              Comment

                              Latest Articles

                              Collapse

                              • seqadmin
                                Strategies for Sequencing Challenging Samples
                                by seqadmin


                                Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
                                03-22-2024, 06:39 AM
                              • seqadmin
                                Techniques and Challenges in Conservation Genomics
                                by seqadmin



                                The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                                Avian Conservation
                                Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                                03-08-2024, 10:41 AM

                              ad_right_rmr

                              Collapse

                              News

                              Collapse

                              Topics Statistics Last Post
                              Started by seqadmin, Yesterday, 06:37 PM
                              0 responses
                              11 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, Yesterday, 06:07 PM
                              0 responses
                              10 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 03-22-2024, 10:03 AM
                              0 responses
                              51 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 03-21-2024, 07:32 AM
                              0 responses
                              67 views
                              0 likes
                              Last Post seqadmin  
                              Working...
                              X