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:
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?
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…..
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?
Comment