<?xml version='1.0'?><rss version="2.0" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:georss="http://www.georss.org/georss" xmlns:atom="http://www.w3.org/2005/Atom" >
<channel>
	<title><![CDATA[BOL: QV calculation in Bash !]]></title>
	<link>https://bioinformaticsonline.com/snippets/view/40386/qv-calculation-in-bash?</link>
	<atom:link href="https://bioinformaticsonline.com/snippets/view/40386/qv-calculation-in-bash?" rel="self" type="application/rss+xml" />
	<description><![CDATA[]]></description>
	
	<item>
	<guid isPermaLink="true">https://bioinformaticsonline.com/snippets/view/40386/qv-calculation-in-bash</guid>
	<pubDate>Fri, 13 Dec 2019 21:14:16 -0600</pubDate>
	<link>https://bioinformaticsonline.com/snippets/view/40386/qv-calculation-in-bash</link>
	<title><![CDATA[QV calculation in Bash !]]></title>
	<description><![CDATA[<code># $1 = vcf file
# $2 = input bam file
# $3 = output QV file

module load samtools

NUM_BP=`samtools depth $2 | perl -e &#039;$c = 0; while(&lt;&gt;){chomp; @s = split(/\t/); if(scalar(@s) &gt;= 3){$c++;}} print &quot;$c\n&quot;;&#039;`
echo &quot;num bp: &quot;$NUM_BP

NUM_SNP=`cat $1 |grep -v &quot;#&quot; | awk -F &quot;\t&quot; &#039;{if (!match($NF, &quot;0/1&quot;)) print $1&quot;\t&quot;$2&quot;\t&quot;$3&quot;\t&quot;$4&quot;\t&quot;$5&quot;\t&quot;$8}&#039; | tr &#039;;&#039; &#039; &#039; | sed s/AB=//g | awk -v WEIGHT=0 &#039;{if ($6 &gt;= WEIGHT) print $0}&#039; | awk -v SUM=0 &#039;{if (length($4) == length($5)) { SUM+=length($4); } else if (length($4) &lt; length($5)) { SUM+=length($5)-length($4); } else { SUM+=length($4)-length($5)}} END { print SUM}&#039;`
echo &quot;num snp: &quot;$NUM_SNP

perl -e &#039;chomp(@ARGV); $ns = $ARGV[0]; $nb = $ARGV[1]; print (-10 * log($ns/$nb)/log(10)); print &quot;\n&quot;;&#039; $NUM_SNP $NUM_BP &gt; $3
cat $3</code>]]></description>
	<dc:creator>Shruti Paniwala</dc:creator>
</item>

</channel>
</rss>