Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Samtools on Chromosomes Versus Genome

    Hi everyone,

    Does anybody have experience running samtools on CASAVA alignments?

    I am trying to call SNPs using samtools on a CASAVA/illumina alignment.

    CASAVA outputs an alignment of the type name_export.txt.gz

    and it automatically converts to BAM format when using their SNP caller.

    Now I want to call SNPs using Samtools.

    The alignment was done using each chromosome individually:

    genomes/Homo_sapiens/UCSC/hg18/Sequence/Chromosomes/Chr...

    However, the BAM alignment provided by CASAVA is for the whole genome, so then should I use the whole reference genome?:

    samtools mpileup -f ../genomes/Homo_sapiens/UCSC/hg18/Sequence/WholeGenomeFasta/genome.fa sorted.bam > sorted.mpileupOutput

    the problem is that, in that case the mpileup output doesn't seem to recognize it as I get:

    bash-3.2$ head -n 100 sorted.mpileupOutput
    chr1.fa 3309 N 1 ^NT C
    chr1.fa 3310 N 1 G C
    chr1.fa 3311 N 1 C F
    chr1.fa 3312 N 1 C F
    chr1.fa 3313 N 1 C F


    Any ideas??

    Thanks in advance,

    Ramiro

  • #2
    try this. mpileup doesn't do too much. bcftools is the program that actually does the SNP calling.

    Code:
    samtools mpileup -uf ../genomes/Homo_sapiens/UCSC/hg18/Sequence/WholeGenomeFasta/genome.fa sorted.bam | bcftools view -bvcg - > sorted.bcf
    bcftools view sorted.bcf > sorted.vcf
    head sorted.vcf
    /* Shawn Driscoll, Gene Expression Laboratory, Pfaff
    Salk Institute for Biological Studies, La Jolla, CA, USA */

    Comment


    • #3
      or you can try splitting the bam file by chromosome and then calling SNP on them separately.. you could also try starting from the fastq or qseq files, run alignment using something like bwa and then call variants
      --
      bioinfosm

      Comment


      • #4
        indeed...i align data using BWA when i'm looking for SNPs.
        /* Shawn Driscoll, Gene Expression Laboratory, Pfaff
        Salk Institute for Biological Studies, La Jolla, CA, USA */

        Comment


        • #5
          How to combine .vcf files into whole genome

          This discussion is very useful, but when I run mpileup plus bcftools to produce .vcf files for individual Chrm, I have found no way to combine them into whole genome .vcf file for further processing. vcf-concat seems like it might do it, but it never runs due to diverse errors. Most recent is:
          ****************************
          The use of -1 for unknown number of values is deprecated, please use '.' instead.
          FORMAT=<ID=PL,Number=-1,Type=Integer,Description="List of Phred-scaled genotype likelihoods, number of values is (#ALT+1)*(#ALT+2)/2">
          ****************************
          I am running most recent version of samtools (0.1.18) -- has anyone else noticed a mismatch in the FORMAT specs between samtools and vcftools??

          Any suggestions at all would be very welcom, this is driving me nuts.

          Thanks
          Allan Lindh

          Comment


          • #6
            sdriscoll,

            It seems that mpileup looking like:

            chr1.fa 3309 N 1 ^NT C
            chr1.fa 3310 N 1 G C
            chr1.fa 3311 N 1 C F
            chr1.fa 3312 N 1 C F

            is a bad sign, it looks to me that it is not mapping to the ref. genome properly. I did as you suggested and everything was a SNP.

            Comment


            • #7
              Double check that the chromosome names in your .bam file match the chromosome names in your fasta. And I'd remake the fa.fai with samtools faidx, see if that completes without an error.

              Comment


              • #8
                I am currently trying to deal with exactly this problem: CASAVA alignments that choke mpileup. We use mpileup output directly for analysis with our own SNP finders, and it works very well with output from a BWA-based aligner (stampy). But with CASAVA output, as stated in another response, mpileup thinks everything is a SNP, all reads with phred scores of 0. This is not the case, fortunately, as the same read data run through stampy produces an alignment file that behaves well with mpileup.

                Comment


                • #9
                  I must admit I just gave up on Casava alignments, just use BWA and I am exploring other aligners.

                  Comment


                  • #10
                    I like stampy-BWA, but it does take a lot of cpu time

                    Comment

                    Latest Articles

                    Collapse

                    • seqadmin
                      Exploring the Dynamics of the Tumor Microenvironment
                      by seqadmin




                      The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
                      07-08-2024, 03:19 PM
                    • seqadmin
                      Exploring Human Diversity Through Large-Scale Omics
                      by seqadmin


                      In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
                      06-25-2024, 06:43 AM

                    ad_right_rmr

                    Collapse

                    News

                    Collapse

                    Topics Statistics Last Post
                    Started by seqadmin, 07-19-2024, 07:20 AM
                    0 responses
                    25 views
                    0 likes
                    Last Post seqadmin  
                    Started by seqadmin, 07-16-2024, 05:49 AM
                    0 responses
                    41 views
                    0 likes
                    Last Post seqadmin  
                    Started by seqadmin, 07-15-2024, 06:53 AM
                    0 responses
                    45 views
                    0 likes
                    Last Post seqadmin  
                    Started by seqadmin, 07-10-2024, 07:30 AM
                    0 responses
                    42 views
                    0 likes
                    Last Post seqadmin  
                    Working...
                    X