update norm/uniform

variants/norm_annot.sh

+#!/bin/bash
+#union.sh
+
+usage() {
+  echo "-h Help documentation for gatkrunner.sh"
+  echo "-r  --Reference Genome: GRCh38 or GRCm38"
+  echo "-p  --Prefix for output file name"
+  echo "-v  --VCF File" 
+  echo "Example: bash union.sh -p prefix -r /path/GRCh38"
+  exit 1
+}
+OPTIND=1 # Reset OPTIND
+while getopts :r:p:v:h opt
+do
+    case $opt in
+        r) index_path=$OPTARG;;
+        p) pair_id=$OPTARG;;
+	v) vcf=$OPTARG;;
+        h) usage;;
+    esac
+done
+function join_by { local IFS="$1"; shift; echo "$*"; }
+shift $(($OPTIND -1))
+baseDir="`dirname \"$0\"`"
+
+source /etc/profile.d/modules.sh
+module load bedtools/2.26.0 samtools/1.6 bcftools/1.6 snpeff/4.3q 
+
+prefix="${vcf%.vcf.gz}"
+perl $baseDir\/uniform_vcf_gt.pl $vcf
+bgzip -f ${prefix}.uniform.vcf
+j=${prefix}.uniform.vcf.gz
+tabix -f $j
+bcftools norm -m - -O z -o ${prefix}.norm.vcf.gz $j

variants/uniform_vcf_gt.pl

@@ -9,6 +9,12 @@ open VCF, "gunzip -c $vcf|" or die $!;
 while (my $line = <VCF>) {
     chomp($line);
     if ($line =~ m/#/) {
+	if ($line =~ m/#CHROM/) {
+	    print OUT "##FORMAT=<ID=AO,Number=A,Type=Integer,Description=\"Alternate allele observation count\">\n";
+	    print OUT "##FORMAT=<ID=RO,Number=1,Type=Integer,Description=\"Reference allele observation count\">\n";
+	    print OUT "##FORMAT=<ID=AD,Number=R,Type=Integer,Description=\"Allelic depths for the ref and alt alleles in the order listed\">\n";
+	    print OUT "##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"Approximate read depth (reads with MQ=255 or with bad mates are filtered)\">\n";
+	}
 	print OUT $line,"\n";
 	next;
     }
@@ -34,11 +40,6 @@ while (my $line = <VCF>) {
       foreach my $i (0..$#deschead) {
 	  $gtdata{$deschead[$i]} = $gtinfo[$i];
       }
-      if ($gtdata{DP} == 0 || $gtdata{GT} eq './.') {
-	  push @newgts, '.:.:.:.:.';
-	  $missingGT ++;
-	  next FG;
-      }
       if ($gtdata{AD}){
 	  ($gtdata{RO},@alts) = split(/,/,$gtdata{AD});
 	  $gtdata{AO} = join(",",@alts);
@@ -60,6 +61,11 @@ while (my $line = <VCF>) {
       if ($gtdata{DP} && $gtdata{DP} < 5) {
 	  $missingGT ++;
       }
+      if ($gtdata{DP} == 0 || $gtdata{GT} eq './.') {
+	  push @newgts, '.:.:.:.:.';
+	  $missingGT ++;
+	  next FG;
+      }
       push @newgts, join(":",$gtdata{GT},$gtdata{DP},$gtdata{AD},$gtdata{AO},$gtdata{RO});
   }
     next if ($missingGT == scalar(@gts));