SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Minimum short read required for transcriptome assembly edge Bioinformatics 4 08-25-2013 10:16 PM
minimum depth of coverage for transfrag extension apadr007 Bioinformatics 6 05-15-2012 07:23 AM
minimum depth variant calling samtools/gatk m_elena_bioinfo Bioinformatics 1 12-06-2011 09:31 AM
What's maximum depth and minimum depth? slny Bioinformatics 0 03-17-2011 11:13 AM
About the read depth of coverage El Mariachi Illumina/Solexa 2 12-30-2010 01:22 AM

Reply
 
Thread Tools
Old 10-24-2011, 08:19 AM   #1
Lspoor
Member
 
Location: Fife, Scotland

Join Date: Mar 2010
Posts: 18
Question vcfutils.pl -specifying minimum read depth?



Hi everyone, I was wondering if you might be able to help me with specifying a minimum read depth when calling consensus sequences?

I'm running samtools-0.1.18 on Ubuntu 10.04. I have produced VCF files for all sites for my bacterial genome sequences (from BWA BAM files), which seem ok (based on using the following commandline: /samtools-0.1.18/samtools mpileup -E -q30 -Q30 -uf REF.fa isolate1.merged.bam | /samtools-0.1.18/bcftools/bcftools view -cg - > isolate1.all_sites.vcf), and am using the following command line in vcftools with the aim of specifying a minimum read depth of 3 when calling the genome consensus sequences:

/samtools-0.1.18/bcfutils/vcfutils.pl vcf2fq -d3 AB10.all_sites_E_q30_Q30.vcf > AB10.all_sites_E_q30_Q30_RD3.fq (which I convert to FASTA using a python script).

Whether or not I specify the option -d3 (I believe 3 is the default anyway) in the commandline, in all cases where the VCF file states that DP=1 or DP=2, rather than not calling the base, it calls a lowercase letter at these sites.

Over the course of my genomes, this amounts to approximately 48 000 bases in each genome which have been called lowercase, which includes approximately 650 variant sites. If these variant sites are of poor quality and include potentially false positive SNPs, this may be adversely affecting my downstream phylogenetic analysis.

My questions are:
a) Does anyone know if this is normal for it to insert a lowercase base rather than calling 'n' at these sites (which is what I was incorrectly assuming it would do!)
b) What are your opinions on whether I should just replace all of these sites with an 'n' to be on the conservative side?

Any help is much appreciated, SeqAnswers has been so great at helping in the past!

Best Wishes,

Laura
Lspoor is offline   Reply With Quote
Old 03-08-2012, 04:06 AM   #2
wzhjlau2009
Junior Member
 
Location: china

Join Date: Jun 2011
Posts: 5
Default

how about your results finally ? i am doing the same process as your had done . i want to know may i have a chance to have a communication with you about this area.my email:wzhjlau@163.com QQ:461739553
best wishes,
wisdom
wzhjlau2009 is offline   Reply With Quote
Old 03-08-2012, 09:52 AM   #3
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Here's the relevent part of the script:

Quote:
Options: -d INT minimum depth [$opts{d}]
-D INT maximum depth [$opts{D}]
-Q INT min RMS mapQ [$opts{Q}]
-l INT INDEL filtering window [$opts{l}]
my $_Q = $opts{Q};
my $_d = $opts{d};
my $_D = $opts{D};
my %het = (AC=>'M', AG=>'R', AT=>'W', CA=>'M', CG=>'S', CT=>'Y',
GA=>'R', GC=>'S', GT=>'K', TA=>'W', TC=>'Y', TG=>'K')

...

my ($ref, $alt) = ($t[3], $1);
my ($b, $q);
$q = $1 if ($t[7] =~ /FQ=(-?[\d\.]+)/);
if ($q < 0) {
$_ = ($t[7] =~ /AF1=([\d\.]+)/)? $1 : 0;
$b = ($_ < .5 || $alt eq '.')? $ref : $alt;
$q = -$q;
} else {
$b = $het{"$ref$alt"};
$b ||= 'N';
}
$b = lc($b);
$b = uc($b) if (($t[7] =~ /MQ=(\d+)/ && $1 >= $_Q) && ($t[7] =~ /DP=(\d+)/ && $1 >= $_d && $1 <= $_D));
So basically, the only way an N gets put in is if the FQ is < 0 (which happens when your SNP is mixed), and there's no single letter code for that mix.

Otherwise, it will be lowercase unless the MQ and DP are within the ranges you set in the options, in which case, it's uppercase.
swbarnes2 is offline   Reply With Quote
Old 05-27-2013, 02:20 AM   #4
ymwur
Member
 
Location: Taiwan

Join Date: Nov 2012
Posts: 11
Default

Hi,
So, I am wodnering what does a N in the concesus generated by mpileup measns?
I am still confused about when FQ can be samller than 0? I am not sure what does mixed SNP mean?
Does a N mean there is no read mapped to that particular site (if the bam is resulted by BWA mapping)?

Thans,

Chih-Ming
ymwur is offline   Reply With Quote
Reply

Tags
consensus, minimum read depth, vcfutils.pl

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 07:15 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO