Header Leaderboard Ad

Collapse

How Varscan process the pileup data?

Collapse

Announcement

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

  • How Varscan process the pileup data?

    Hi,

    if any of you has used Varscan to find somatic mutations given normal and tumor pileup files, maybe you can give me some light about how it processes the data.

    For instance, given the following position in the tumor pileup file:

    Code:
    chr1	173824	T	21	GG...*..*.G.......GG.	
    				##9!!!!;!<#;<<=<<<!#<
    Varscan says that 8 tumor reads supports the 'T' and 2 tumor reads supports the 'G'. I do not see how it is fitted!

    Another random example:

    Code:
    chr1	30028	c	9	,tt,,,T,,	
    				5::778?.4
    Varscan says that in the previous position, there are 3 reads for nucleotide 'T' and 5 for 'C' (why not 6??). There are these kind of discordances in almost all the positions I've checked, so I'm pretty confused.

    thanks!
    david

  • #2
    VarScan includes a minimum base quality threshold requirement in order to count reads. By default, this is 15 or 20 depending on the command you use.

    In the first example, some bases have a base quality represented by '#', indicating a Phred score of 3. In the second example, one reference base has a quality represented by ".", which is 14. Both of these fall below the minimum base quality threshold in VarScan.

    Comment


    • #3
      Thank you for your response, dkoboldt.

      I indeed thought in the possibility that some 'default quality filtering' was being applied, but (for instance) in the first example, all G's should be not taken into account since their qualities are either '#' or '!'. However, Varscan stats say that for that position there are two G's. This does not seem to fit.

      There is plenty of things like that, and I am not able to find a 'pattern' of quality filtering by checking the pileup positions as compared to the corresponding VArscan entry...

      thanks a lot,
      cheers
      David
      Last edited by david.tamborero; 01-04-2012, 06:24 PM.

      Comment


      • #4
        Just for those that can be interested in the varscan filtering according to the base quality: apparently Varscan has a bug that arises when certain characters appear in the field containing the bases (for instance, the '*'). This provokes a shift in the 'base-base quality' correspondance and then the counting of bases supporting each allele becomes wrong.

        Therefore, the base quality filtering performed by Varscan fails for those reads containing these characters. I've written to the authors, hope it helps.

        cheers,
        david

        Comment


        • #5
          Hello David,
          I intended to try Varscan with normal and tumoral files, so thanks for mentionning the problem. Please, keep us informed if you got news from the authors. I would like to try the tool when the problem is solved!
          Do you know if the problem is also present when using Varsan with only one file?

          Comment


          • #6
            No idea, I only have used the somatic command. You can try the other Varscan commands (as the pileup2snp) with a dummy pileup file and play with the --min-avg-qual parameter to check it out.

            By the way, I need that Varscan performs the quality filtering because the samtools mpileup command also does it wrong (the -Q argument does not work). Seems that the base quality filtering is not that important!

            cheers,
            David

            Comment


            • #7
              Ok, thanks!

              Comment


              • #8
                I tried VarScan in the somatic mode and I also experienced this problem.
                For example, VarScan gives as output the following number of reads:

                - for the normal sample: 183 (reference) and 4 (variant)
                - for the tumor sample: 8 (reference) and 14 (variant)

                This results in a somatic mutation. I checked my bam files in IGV and I found:
                - for the normal sample: 185 (reference) and 165 (variant)
                - for the tumor sample: 8 (reference) and 359 (variant)

                Of course, some reads can be deleted because they have low quality, but I doubt that it could be the cases for so many reads... This variation is rather a LOH rather than a somatic mutation.
                And I've got this kind of results several times...

                I wonder how this problem affect the SNPs detection and how wrong is my analysis...

                Do you know if this problem will be taken into account?

                Comment


                • #9
                  I guess than before running varscan you have converted the bam file to pileup file, right?

                  If so, such 'lost' reads could be explained by the samtools mpileup command. Check its arguments: the mapping_quality filtering works (as far as i remember, and as opposite to the base_quality filtering), therefore some reads can be removed due to this parameter. The other source of reads removal during bam to pileup conversion can be the 'anomalous read pairs': check the '-A' argument of the mpileup command to see if reads counts are more consistent with IGV data (which opens the bam file).

                  Comment


                  • #10
                    Originally posted by david.tamborero View Post
                    I guess than before running varscan you have converted the bam file to pileup file, right?

                    If so, such 'lost' reads could be explained by the samtools mpileup command. Check its arguments: the mapping_quality filtering works (as far as i remember, and as opposite to the base_quality filtering), therefore some reads can be removed due to this parameter. The other source of reads removal during bam to pileup conversion can be the 'anomalous read pairs': check the '-A' argument of the mpileup command to see if reads counts are more consistent with IGV data (which opens the bam file).
                    Exactly, I have converted my bam files into pileup files to use VarScan.

                    samtools mpileup -f ~/fasta/hg19.fasta /data/patient1/s_garma-296_converted_sorted.bam > /data/patient1/s_garma-296_converted_sorted.pileup
                    The default arguments are:
                    -q INT skip alignments with mapQ smaller than INT [0]
                    -Q INT skip bases with baseQ/BAQ smaller than INT [13]
                    I checked on IGV the PHRED of some of the reads at this position and they were good (>30). So the missing reads are probably not due to base quality. I don't know for the mapping quality... How can I check it?

                    As you suggested me, I re-run samtools mpileup with -A:
                    samtools mpileup -f ~/fasta/hg19.fasta /data/patient1/s_garma-296_converted_sorted.bam > /data/patient1/s_garma-296_converted_sorted.pileup
                    And I have still the same kind of problems: LOH considered as somatic mutations because of missing reads.
                    I am not sure what is doing the -A option... Does it add some anomalous reads?
                    Last edited by Jane M; 02-06-2012, 08:52 AM.

                    Comment


                    • #11
                      Jane,

                      It will also help if you provide the pileup output for 1-2 examples of this base dropout - we can quickly look at what's in it to see if VarScan is not counting anything.

                      Please note that the base quality parsing issue was resolved as of VarScan v2.2.8. If any of you encounter it, please let me know!

                      Yours,

                      Dan Koboldt

                      Comment


                      • #12
                        Thank you for your answer Dan,

                        Here is the output of the pileup file at the position where I have the problem I mentioned:
                        With VarScan somatic:
                        - for the normal sample: 183 (reference G) and 4 (variant A)
                        - for the tumor sample: 8 (reference) and 14 (variant)

                        With IGV:
                        - for the normal sample: 185 (reference) and 165 (variant)
                        - for the tumor sample: 8 (reference) and 359 (variant)

                        [PCJane patient1]$ grep 186380515 s_garma-fibros_converted_sorted.pileup
                        chr4 186380515 G 350 .$,$a$,$,$a$a$.Aa,aA.aaa,..,a,.....aa,,,a,,AAaaaA.AA,,Aa,,aa...a,..a,,.AAaaa,a...,,a,,,..A,a,a,.,,.Aaaaa,,,.,Aa.,Aa,,Aa,AAa,aa..a,,a..a,aAa,aA.Aaa.aaAAAA.A,a,aA,A,,,..aa,aaa.A.aaa..aa...A,,,aA.AA..Aa,,aa...a,aaaAA.AaAa,a,AAA.,a,AAA....a..,aa,AA,A.a..aA,aaa.,,aA..AAa,aA.,aA,,AA,,,,,AaAAA,A,,,A.aaa..aa,AAa,,,..A,,,,A.a,a,aa,,....aaa,....,,,,..A..,aa......^~A^~, [email protected]!!B([email protected][email protected]:FF75HH55JIJ5GJJ9HH733555H5GII:J3JIJBH4J3J3JIJJJ53333JIJEJ43IJ43JJ43J334F33EJ4JJ3JI3J333J;3J433J333363J4J3F43J6JGIJJ3:J333?3J333JJ33JGJ4JJJ84J33JJ33GJ33JJJ3J:3346J3343J3J>33JJ*J334IJGJ4JJJ34G33F3J4JJ46J433JDJ33IJ333I33JJ33JJ43JIJJJ43334[email protected]CCCC%E
                        ^C
                        [PCJane patient1]$
                        I intend to test Varscan in the "simple mode" today to check if I have this issue too.

                        When using the CASAVA pipeline on my normal sample, I got:

                        186380515 chr4 (G)155 (A) 142
                        By comparison between, IGV, CASAVA and Varscan in the normal sample, I have:
                        -for the reference: 185, 155, 183
                        -for the variant: 165, 142, 4

                        Moreover, I tried JointSNVMix, which let the direct comparison between normal and tumoral sample. At this specific position, I have nothing. But at several sites, where I see this problem in VarScan ,my results (number of read counts) with JointSNV are closed to the ones in IGV...
                        Last edited by Jane M; 03-02-2012, 01:39 AM.

                        Comment


                        • #13
                          I ran VarScan in the simple mode on the normal sample, using pileup2snp:

                          java -Xmx20g -jar VarScan.v2.2.8.jar pileup2snp /data/patient1/s_garma-fibros_converted_sorted.pileup --min-coverage 10 --min-avg-qual 25 > /data/patient1/output_varscan_snp_garma.snp
                          At the position where I have a problem in read counts, I got:
                          Code:
                          chr4  	186380515	(G)  184	   (A) 5
                          I have the same problem in both mode. It's very strange... It cannot be related to the version of VarScan since I'm using the latest one, v2.2.8...

                          I attach here the screenshots from IGV at this specific position.
                          Attached Files

                          Comment


                          • #14
                            Hello,

                            I just wanted to let you know that there are buggs in VarScan, in both default and somatic modes. There are other people experiencing the same problem like me.
                            To conclude, I give up VarScan but I hope the problem will be taken into account, because appart this issue, this seems to be an interesting tool!

                            Comment


                            • #15
                              Hello,

                              If some users are experiencing the same issue and read this subject, the problem comes from samtools. Add the -B or -E option to disable the BAQ computation.

                              Comment

                              Latest Articles

                              Collapse

                              • seqadmin
                                A Brief Overview and Common Challenges in Single-cell Sequencing Analysis
                                by seqadmin


                                ​​​​​​The introduction of single-cell sequencing has advanced the ability to study cell-to-cell heterogeneity. Its use has improved our understanding of somatic mutations1, cell lineages2, cellular diversity and regulation3, and development in multicellular organisms4. Single-cell sequencing encompasses hundreds of techniques with different approaches to studying the genomes, transcriptomes, epigenomes, and other omics of individual cells. The analysis of single-cell sequencing data i...

                                01-24-2023, 01:19 PM
                              • seqadmin
                                Introduction to Single-Cell Sequencing
                                by seqadmin
                                Single-cell sequencing is a technique used to investigate the genome, transcriptome, epigenome, and other omics of individual cells using high-throughput sequencing. This technology has provided many scientific breakthroughs and continues to be applied across many fields, including microbiology, oncology, immunology, neurobiology, precision medicine, and stem cell research.

                                The advancement of single-cell sequencing began in 2009 when Tang et al. investigated the single-cell transcriptomes
                                ...
                                01-09-2023, 03:10 PM
                              • seqadmin
                                AVITI from Element Biosciences: Latest Sequencing Technologies—Part 6
                                by seqadmin
                                Element Biosciences made its sequencing market debut this year when it released AVITI, its first sequencer. The AVITI System uses avidity sequencing, a novel sequencing chemistry that delivers higher quality data, decreases cycle times, and requires lower reagent concentrations. This new instrument reportedly features lower operating and start-up costs while maintaining quality sequencing.

                                Read type and length
                                AVITI is a short-read benchtop sequencer that also offers an innovative...
                                12-29-2022, 10:44 AM

                              ad_right_rmr

                              Collapse
                              Working...
                              X