Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • multisample variant calling and consensus sequences

    Hi,

    I have used BWA to map reads from ~45 individuals of my favorite species to a reference (haploid) mitochondrial DNA sequence. My goal is to get a consensus sequence for each individual, which I then want to use in a skyline plot of effective population size in BEAST.

    In the past I have used samtools/bcftools/vcfutils.pl to create a consensus sequence. Using these commands it seems like you make a pileup, call variant sites, and then convert to fastq.

    A couple of my individuals have low depth, and I have been considering how to handle them (without throwing them out). I thought perhaps I could follow the general steps above and do a multisample variant call, filter the variants, and then make the fastq for each individual.

    I have tried several different things and none have worked. Here is one example:

    samtools-1.2 mpileup -ugf Ref.fasta Individual1.bam Individual2.bam Individual3.bam Individual4.bam Individual5.bam | bcftools-1.2 call -vmO z –S samples.txt -o Variants.vcf.gz

    bcftools-1.2 filter -e 'QUAL<30 | MIN(DP)<10' Variants.vcf.gz >Variants_Filter.vcf

    vcfutils-1.2.pl subsam Variants_Filter.vcf Individual1.bam > Individual1.vcf

    vcfutils-1.2.pl vcf2fq -d 2 Individual1.vcf > Individual1.fastq

    This gives:
    Code:
    @Individual1
    nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnAnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnMnnnn…..
    where the A is call of the alternate allele in this individual and the M is actually a well-supported call of the reference sequence at this site (not sure why it is given as an ambiguity). I guess I have lost the information about the non-variant sites along the way.

    Any suggestions for a better way to go about handling my low depth individuals?

    Beyond trying the multisample variant calling, I also tried to get positions in the consensus sequence with depth less than some minimum to be output as N, but the best I can do is get them output in lower case from vcfutils (presumably I would have to use regex to change them to N)

    I have thought about simply using the SNPs in BEAST (as opposed to the consensus sequence) but for effective population sizes, it seems like the relative numbers of variable and invariable sites might have some information?
    Last edited by Brian Bushnell; 08-07-2015, 08:43 PM. Reason: formatting

  • #2
    What kind of depth do you have? Mapping the reads to mito with very high sensitivity, keeping those that map, and assembling them with a normal assembler should work.

    Comment


    • #3
      Most are > 10x, but I have two at 5x and two at 2x. These lower depth individuals are from ancient DNA so I merged the reads and then mapped. I can do tests about the impact of missing data, but I want to make sure the positions that are called are as accurate as possible.

      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, Today, 11:09 AM
      0 responses
      14 views
      0 likes
      Last Post seqadmin  
      Started by seqadmin, 07-19-2024, 07:20 AM
      0 responses
      146 views
      0 likes
      Last Post seqadmin  
      Started by seqadmin, 07-16-2024, 05:49 AM
      0 responses
      120 views
      0 likes
      Last Post seqadmin  
      Started by seqadmin, 07-15-2024, 06:53 AM
      0 responses
      111 views
      0 likes
      Last Post seqadmin  
      Working...
      X