Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Advice on tools sought - exporting consensus sequences with coverage & SNPs?

    Hi all!

    I need some tools advice.
    So; I have mapped reads in BAM format.
    I have called SNPs using samtools/bcftools.

    Are there any good tools for exporting consensus sequences from; say a vcf file; in; say fasta format; that maintain infos like coverage?

    Hope I am making sense here. Feel free to point me to appropriate threads on seqanswers; my crappy searching has not turned up anything!

  • #2
    Okay: after a lot of reading around here on seqanswers, I think I have managed to some of this with samtools - always helps to read the manual carefully.

    Only thing where I am stuck is the consensus. My current read library has okay-ish, but not great, coverage which means there are gaps in my alignment to the reference. These gaps might be there because of extensive deletion or something: but lets not even worry about that now...

    So, if I call a consensus from the vcf file; how are the regions having no coverage treated? Is the sequence just taken from the ref sequence? Can I track/visualise that somehow? Just so I know which bit in the consensus came from the actual bam (i.e. has experimental evidence from my sequencing runs) and which bit came from the reference and has no experimental evidence?

    Comment


    • #3
      Some progress.

      Apologies if it seems like I am talking to myself; but maybe at some point someone will search for the same issue and this may help...

      So; after using the vcfutils and f2q conversion in samtools; I have ended up with a consensus fastq file.

      The file contains: Upper case letters (good quality?); lower case letters (low quality?) and - most interestingly - a bunch of Ns (lower case) where my reads do not cover the reference.

      After converting this to fasta (seqtk) - which obviously looses the upper/lower case info - I have something I can work with reasonably well for now!

      2 Questions remain for me:

      1. Is there a description of what vcfutils actually does in detail somewhere? So, when does it call upper/lower case?
      2. Is there any way to force it to use the reference sequence when it does not find coverage? So; it prints the ref sequence instead of an N?

      Maybe I need another tool for that.

      Comment


      • #4
        Originally posted by TabeaK View Post
        Some progress.

        Apologies if it seems like I am talking to myself; but maybe at some point someone will search for the same issue and this may help...

        So; after using the vcfutils and f2q conversion in samtools; I have ended up with a consensus fastq file.

        The file contains: Upper case letters (good quality?); lower case letters (low quality?) and - most interestingly - a bunch of Ns (lower case) where my reads do not cover the reference.

        After converting this to fasta (seqtk) - which obviously looses the upper/lower case info - I have something I can work with reasonably well for now!

        2 Questions remain for me:

        1. Is there a description of what vcfutils actually does in detail somewhere? So, when does it call upper/lower case?
        2. Is there any way to force it to use the reference sequence when it does not find coverage? So; it prints the ref sequence instead of an N?

        Maybe I need another tool for that.
        I have an older version of samtools but vcfutils is in perl, so it's human readable. If you are using vcf2fastq, the key lines are these:

        Code:
        my ($b, $q);
        	  $q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);
        	  if ($q < 0) {
        		$_ = ($t[7] =~ /AF1=([\d\.]+)/)? $1 : 0;
        		[B]$b = ($_ < .5 || $alt eq '.')? $ref : $alt;[/B]
        		$q = -$q;
        	  } else {
        		[B]$b = $het{"$ref$alt"};
        		$b ||= 'N';[/B]
        	  }
        	  $b = lc($b);
        	  $b = uc($b) if (($t[7] =~ /MQ=(\d+)/ && $1 >= $_Q) && ($t[7] =~ /DP=(\d+)/ && $1 >= $_d && $1 <= $_D));
        $b is depending on the vcf entry, either the reference letter, or the alt letter, or a letter representing both letters, if the vcf thinks its a mixed SNP, or an N, if there's nothing at all for that letter. Those lines are the first 3 bolded lines.

        The last line determines if $b is uppercase or lowercase, which basically happens if the mapping quality is above the minimum you specified, and the coverage is between the min and max you specified.

        The other place in the program where it will make your sequence lowercase is here:

        Code:
        for my $g (@$gaps) {
        	my $beg = $g->[0] > $l? $g->[0] - $l : 0;
        	my $end = $g->[0] + $g->[1] + $l;
        	$end = length($$seq) if ($end > length($$seq));
        	substr($$seq, $beg, $end - $beg) = lc(substr($$seq, $beg, $end - $beg));
          }
        And @gaps is populated here:

        Code:
        	} elsif ($t[4] ne '.') { # an INDEL
        	  push(@gaps, [$t[1], length($t[3])]);
        	}
        So around indels, the software will make a window of lowercase letters.

        And if you have no coverage for a region, I'm not sure you necessarily want the software to just tell you that it's reference, because what if there's a deletion?

        Comment


        • #5
          If you have a haploid sample, and only one sample, and the VCF reflects that fact (e.g. you called with freebayes --ploidy 1) then you can use vcfgeno2haplo (https://github.com/ekg/vcflib/blob/m...geno2haplo.cpp) with an arbitrary window size to get a consensus sequence. It will be recorded in your "ALT" in the only line of the output VCF. This works with indels.

          If you have diploid samples, you'll have to phase them somehow. This method won't work if your samples aren't phased.

          For coverage, you can get statistics in a number of ways. I'm not sure the best for your usage.

          Comment

          Latest Articles

          Collapse

          • seqadmin
            Essential Discoveries and Tools in Epitranscriptomics
            by seqadmin




            The field of epigenetics has traditionally concentrated more on DNA and how changes like methylation and phosphorylation of histones impact gene expression and regulation. However, our increased understanding of RNA modifications and their importance in cellular processes has led to a rise in epitranscriptomics research. “Epitranscriptomics brings together the concepts of epigenetics and gene expression,” explained Adrien Leger, PhD, Principal Research Scientist...
            Yesterday, 07:01 AM
          • seqadmin
            Current Approaches to Protein Sequencing
            by seqadmin


            Proteins are often described as the workhorses of the cell, and identifying their sequences is key to understanding their role in biological processes and disease. Currently, the most common technique used to determine protein sequences is mass spectrometry. While still a valuable tool, mass spectrometry faces several limitations and requires a highly experienced scientist familiar with the equipment to operate it. Additionally, other proteomic methods, like affinity assays, are constrained...
            04-04-2024, 04:25 PM

          ad_right_rmr

          Collapse

          News

          Collapse

          Topics Statistics Last Post
          Started by seqadmin, 04-11-2024, 12:08 PM
          0 responses
          57 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 04-10-2024, 10:19 PM
          0 responses
          53 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 04-10-2024, 09:21 AM
          0 responses
          45 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 04-04-2024, 09:00 AM
          0 responses
          55 views
          0 likes
          Last Post seqadmin  
          Working...
          X