Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Does Samtools mpileup command require a reference fasta?

    I performed an experiment to track mutation accumulation over time in a fungal population. I sequenced several individual genomes from different time points and would like to extract the SNP sites only.

    I currently have transformed my raw reads into *.sorted.bam files. I would like to use mpileup without a reference .fasta file because I'm only interested in the SNPs among *my* genomes, not between my genomes and the published reference genome. Is this possible?

    When I use the following commands without a supplied "-f ref.fasta" switch:
    samtools mpileup -ug sorted.bam1..3 | bcftools view -bcvg - > mpileup_noref.bcf
    bcftools view mpileup_noref.bcf | vcfutils varFilter > mpileup_noref.vcf

    I get zero SNPs in the .vcf file, though the .bcf file is enormous.

    When I use the commands:
    samtools mpileup -ugf ref.fasta sorted.bam1..3 | bcftools view -bcvg - > mpileup_ref.bcf
    bcftools view mpileup_ref.bcf | vcfutils varFilter > mpileup_ref.vcf

    I get 1700 SNPs in the .vcf, and the .bcf is much smaller than the no_ref.bcf file.

    Does anyone have any idea what is causing the discrepancy here? I would expect to find more SNPs between my samples and the reference.fasta than I would just among my samples, which are all related.

    Thank you!

  • #2
    Just use the reference fasta. I'm not sure that bcftools is able to work with no reference fasta (though the regular pileup output can) Your vcf will not include all the sites where every sample matches the reference, it will only show the sites where it believes there are differences. That's what the -v switch is for.

    Comment


    • #3
      Thank you If my samples are different from each other they will also be different from the reference, so I should be able to see them in the with_reference.vcf.

      However, what if all of my samples are the same at one particular site, but different from the reference? It seems like I would recover irrelevant information this way, when I am only interested in sites that are different among my individuals and not between them and the ref. Is there a simple way to filter these variants out?

      Comment


      • #4
        You just filter out SNPs found in every sample.

        Or, if your sample numbers are modest, and you are going to be doing more samples from this family, use the SNP data to tweak your reference genome to match the consensus, then realign everything to that consensus, and recall the SNPs.

        Comment


        • #5
          Originally posted by swbarnes2 View Post
          You just filter out SNPs found in every sample.
          How do I do this? Any suggestion? I can't seem to find a solution on how to do this.

          Comment


          • #6
            Are your coding skills modest? I wrote my own script to analyze the .vcf output to find guys who were all the same within the population, but different from the reference (and therefore not particularly interesting to me).

            Comment


            • #7
              Yes, my coding skills are poor. Is it possible to grep for any keyword related to homozygotic variants in the vcf-file?

              Comment


              • #8
                I'm not sure how you could use grep... I'm much more comfortable in perl This is a nice little practice problem for scripting, though, so you should try it on your own at some point.

                That said, working through other people's code was really helpful to me when I was first learning, so I wrote a version out for you. I think it works.
                Code:
                #!/usr/bin/env perl
                # Roxana Capper 28 May 2013, University of Texas at Austin
                
                # -- program description and required arguments
                
                unless ($#ARGV == 1)
                	{
                	print "\ntakes a .vcf file and prints only SNPs that are polymorphic\n";
                	print "WITHIN individuals, not just between all samples and the reference.\n";
                        print "\nUsage:\t script input.vcf polymorphs_only.vcf\n";
                        print "\n"; exit;
                        }
                
                open(VCF, "$ARGV[0]");
                open(OUT, ">$ARGV[1]");
                $one_allele=0; $mult_alleles=0; $total_snps=0;
                
                while (<VCF>)
                	{
                	chomp;
                	$count=0;
                	undef(%genotype);
                	if ($_=~/^#/) {print OUT "$_\n";next;} #print all the header stuff to the new file
                	$line = $_;
                	@columns = split(/\t/,$line);
                	foreach $sample (@columns[9..$#columns])
                		{
                		@ind=split(/:/,$sample); ($allele1, $allele2) = split(/\/|\|/,$ind[0]);
                		if (($allele1 | $allele2) eq ".") {next;}
                		$genotype{$allele1}++;$genotype{$allele2}++;
                		}
                	foreach $variant (keys %genotype) {$count++};  # how many alleles are seen for this SNP?
                	if ($count==1) {$one_allele++};
                	if ($count > 1) {$mult_alleles++; print OUT "$line\n";}  #if polymorphic, print out
                	$total_snps++;
                	}
                
                #output some stats to std out:
                print "There were $total_snps in $ARGV[0].\n";
                print "Of these, $one_allele were monomorphic among the samples\n";
                print "$mult_alleles were printed to $ARGV[1].";

                Comment


                • #9
                  You don't even have to code; if your files are not titanic, you can drop them into Excel, and filter them there. I use conditional formatting to change all the "1/1" into red font, "0/1" into purple, and "0/0" blue. You can then use filters to filter by color. You'd want to get rid of all the entries that are all red.

                  I suppose you could grep for "0/0", that would mean that at least one sample of your bunch matches the reference, making it a meaningful SNP.

                  Comment


                  • #10
                    Roxana, I just tried out your perl script, and it worked perfectly! (didn't have time earlier)
                    Thank you so much! Very useful! I will definitely have more use of your code. Thanks again!

                    Comment

                    Latest Articles

                    Collapse

                    • 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
                    • seqadmin
                      The Impact of AI in Genomic Medicine
                      by seqadmin



                      Artificial intelligence (AI) has evolved from a futuristic vision to a mainstream technology, highlighted by the introduction of tools like OpenAI's ChatGPT and Google's Gemini. In recent years, AI has become increasingly integrated into the field of genomics. This integration has enabled new scientific discoveries while simultaneously raising important ethical questions1. Interviews with two researchers at the center of this intersection provide insightful perspectives into...
                      02-26-2024, 02:07 PM

                    ad_right_rmr

                    Collapse

                    News

                    Collapse

                    Topics Statistics Last Post
                    Started by seqadmin, 03-14-2024, 06:13 AM
                    0 responses
                    33 views
                    0 likes
                    Last Post seqadmin  
                    Started by seqadmin, 03-08-2024, 08:03 AM
                    0 responses
                    72 views
                    0 likes
                    Last Post seqadmin  
                    Started by seqadmin, 03-07-2024, 08:13 AM
                    0 responses
                    80 views
                    0 likes
                    Last Post seqadmin  
                    Started by seqadmin, 03-06-2024, 09:51 AM
                    0 responses
                    68 views
                    0 likes
                    Last Post seqadmin  
                    Working...
                    X