Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • vcfutils.pl -specifying minimum read depth?



    Hi everyone, I was wondering if you might be able to help me with specifying a minimum read depth when calling consensus sequences?

    I'm running samtools-0.1.18 on Ubuntu 10.04. I have produced VCF files for all sites for my bacterial genome sequences (from BWA BAM files), which seem ok (based on using the following commandline: /samtools-0.1.18/samtools mpileup -E -q30 -Q30 -uf REF.fa isolate1.merged.bam | /samtools-0.1.18/bcftools/bcftools view -cg - > isolate1.all_sites.vcf), and am using the following command line in vcftools with the aim of specifying a minimum read depth of 3 when calling the genome consensus sequences:

    /samtools-0.1.18/bcfutils/vcfutils.pl vcf2fq -d3 AB10.all_sites_E_q30_Q30.vcf > AB10.all_sites_E_q30_Q30_RD3.fq (which I convert to FASTA using a python script).

    Whether or not I specify the option -d3 (I believe 3 is the default anyway) in the commandline, in all cases where the VCF file states that DP=1 or DP=2, rather than not calling the base, it calls a lowercase letter at these sites.

    Over the course of my genomes, this amounts to approximately 48 000 bases in each genome which have been called lowercase, which includes approximately 650 variant sites. If these variant sites are of poor quality and include potentially false positive SNPs, this may be adversely affecting my downstream phylogenetic analysis.

    My questions are:
    a) Does anyone know if this is normal for it to insert a lowercase base rather than calling 'n' at these sites (which is what I was incorrectly assuming it would do!)
    b) What are your opinions on whether I should just replace all of these sites with an 'n' to be on the conservative side?

    Any help is much appreciated, SeqAnswers has been so great at helping in the past!

    Best Wishes,

    Laura

  • #2
    how about your results finally ? i am doing the same process as your had done . i want to know may i have a chance to have a communication with you about this area.my email:[email protected] QQ:461739553
    best wishes,
    wisdom

    Comment


    • #3
      Here's the relevent part of the script:

      Options: -d INT minimum depth [$opts{d}]
      -D INT maximum depth [$opts{D}]
      -Q INT min RMS mapQ [$opts{Q}]
      -l INT INDEL filtering window [$opts{l}]
      my $_Q = $opts{Q};
      my $_d = $opts{d};
      my $_D = $opts{D};
      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')

      ...

      my ($ref, $alt) = ($t[3], $1);
      my ($b, $q);
      $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';
      }
      $b = lc($b);
      $b = uc($b) if (($t[7] =~ /MQ=(\d+)/ && $1 >= $_Q) && ($t[7] =~ /DP=(\d+)/ && $1 >= $_d && $1 <= $_D));
      So basically, the only way an N gets put in is if the FQ is < 0 (which happens when your SNP is mixed), and there's no single letter code for that mix.

      Otherwise, it will be lowercase unless the MQ and DP are within the ranges you set in the options, in which case, it's uppercase.

      Comment


      • #4
        Hi,
        So, I am wodnering what does a N in the concesus generated by mpileup measns?
        I am still confused about when FQ can be samller than 0? I am not sure what does mixed SNP mean?
        Does a N mean there is no read mapped to that particular site (if the bam is resulted by BWA mapping)?

        Thans,

        Chih-Ming

        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
        12 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
        68 views
        0 likes
        Last Post seqadmin  
        Working...
        X