Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • mpileup consensus behaviour

    I am wondering if anyone has any experience calling a consensus genotype from mpileup data. I am curious how the suggested command from samtools documentation treats indels and what it does with 'heterozygous' indels.

    Code:
    samtools mpileup -uf ref.fa aln.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq
    Thanks

  • #2
    If anyone is curious - the answer to my question is that this consensus will include NO indels. Any lines in the VCF/BCF that contain indels -- heterozygous or homozygous -- will be ignored.

    I wrote a parser for the vcf created by this command:

    Code:
    $samtools mpileup -uf ref.fa aln.bam | bcftools view -cg -
    that will generate a consensus sequence where heterozygous indels are ignored and heterozygous sites are represented with ambiguity codes. It will also output 2 additional consensus sequences in which heterozygous SNPs and indels are included with gaps to keep the sequences aligned.
    Last edited by rwness; 06-02-2011, 07:32 AM.

    Comment


    • #3
      I've been sticking with the old pileup -c, so I can get the heterozygous sites with ambiguity codes and genotype quality scores. How were you able to get this out of the mpileup/vcf output?

      Comment


      • #4
        The new way posted above, using the all-points vcf and vcfutils.pl will get make ambiguous letters in the consensus. Here's the relevant part of the code:

        my %het = (AC=>'M', AG=>'R', AT=>'W', CA=>'M', CG=>'S', CT=>'Y',
        GA=>'R', GC=>'S', GT=>'K', TA=>'W', TC=>'Y', TG=>'K');


        $q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);
        if ($q < 0) {
        $_ = ($t[7] =~ /AF1=([\d\.]+)/)? $1 : 0;
        $b = ($_ < .5 || $alt eq '.')? $ref : $alt;
        $q = -$q;
        } else {
        $b = $het{"$ref$alt"};
        $b ||= 'N';
        }
        In English, the first line means "check the FQ"; the FQ is negative when the SNP is homozygous, and positive when it's mixed, and the bigger the absolute value, the more confident the SNP.

        So if it's < 0, it does the first part of code (it checks against the AF1, and if the AF1 is > 0.5, which it should be for a homozygous SNP, it sets $b as the alternate letter, if for some reason the AF1 is < .5, it sets $b as the old reference letter.)

        If the FQ is positive, then the SNP should be mixed, and it concatenates the two letters, and checks the hash above to know what single letter to set $b to. Then it adds $b to the growing sequence.

        $q, which is derived from the FQ, ends up being the quality score, though it gets tweaked a little; it adds 33.449 to the figure, then converts it to a letter, capping it at a quality of 126.

        $q = int($q + 33 + .499);
        $q = chr($q <= 126? $q : 126);

        Gaps are handled as they were in the old program, where they are NOT added in, there is just a window of lowercase letters around them. Personally, I made a little perl script, and I feed it the genome, and a conservatively filtered list of SNPs, and I put the changes in that way.

        Comment


        • #5
          Thanks for the nice explanation, swbarnes2. Just want to add that the key purpose of generating consensus is to distinguish regions that cannot be called and regions that are called confidently as homozygous reference. Heterozygotes found in the consensus sequence actually have lower accuracy because fewer filters were applied.

          Comment


          • #6
            I am trying to understand how to set $q. In the line

            $q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);

            Does it mean get $q from command line (1st paramter) in case there is an FQ value?

            Thanks very much
            Last edited by nupurgupta; 07-24-2012, 09:59 AM.

            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
            8 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, Yesterday, 06:07 PM
            0 responses
            8 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, 03-22-2024, 10:03 AM
            0 responses
            49 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, 03-21-2024, 07:32 AM
            0 responses
            66 views
            0 likes
            Last Post seqadmin  
            Working...
            X