Page Comparison

...

Expand

title	One solution...

This linux one-liner should give you a snapshot of data sufficient to figure it out:

Code Block
cat all.samtools.vcf \| head -10000 \| awk '{if ($6>500) {print $2"\t"$10"\t"$11"\t"$12}}' \| grep "0/0" \| sed s/':'/' '/g \| awk '{print $2"\t"$5"\t"$8}' \| tail -100 \| sort \| uniq -c

Expand

title	Here's the output and discussion

Code Block

tacc:/scratch/01057/sphsmith/human_variation$ cat all.samtools.vcf | head -10000 | awk '{if ($6>500) {print $2"\t"$10"\t"$11"\t"$12}}' | grep "0/0" | sed s/':'/' '/g | awk '{print $2"\t"$5"\t"$8}' | tail -100 | sort | uniq -c 
     12 0/0	0/1	0/0
      5 0/0	0/1	0/1
      3 0/1	0/0	0/0
      4 0/1	0/0	0/1
      8 0/1	0/0	1/1
     43 0/1	0/1	0/0
     24 0/1	1/1	0/0
      1 1/1	0/1	0/0

Here are the steps going into this command:

1) Dump the contents of all.samtools.vcf

2) Take the first 10,000 lines

3) If the variant quality score is greater than 500, then print fields 2 (SNP position), 10, 11, and 12 (the 3 genotypes).

4) Filter for only lines that have at least one homozygous SNP (exercise to the reader to understand why...)

5) Break the genotype call apart from other information about depth: "sed" turns the colons into spaces so that awk can just print the genotype fields.

6) Take the last 100 lines, sort them, then count the unique lines

Here is my interpretation of the data:

1) This method effectively looks at a very narrow genomic region, probably within a homologous recombination block.

2) The most telling data: the child will have heterozygous SNPs from two homozygous parents.

3) So all this data is consistent with column 1 (NA12878) being the child:

	 12 0/0	0/1	0/0
      5 0/0	0/1	0/1
      4 0/1	0/0	0/1
      8 0/1	0/0	1/1
     43 0/1	0/1	0/0
     24 0/1	1/1	0/0

"Outlier" data are:

      3 0/1	0/0	0/0

      1 1/1	0/1	0/0

This is, in fact, the correct assessment - NA12878 is the child.

Look at the differences between single- and multiple-sample SNP calls, and between Samtools/Bcftools SNP calls and GATK SNP calls

How many SNPs were called in each case?

Expand

title	Try it your own way, or try this one-line command...

 for file in `ls *.vcf`; do echo "File: $file `cat $file | grep -v '^#' | wc -l`"; done

What's the overlap between all the single-sample SNP calls aggregated together and the multi-sample SNP calls?

Expand

Code Block

cat NA128*samtools.vcf | grep -v '^#' | awk '{print $2}' | sort | uniq > all.single.samtools.snps
cat all.samtools.vcf | grep -v '^#' | awk '{print $2}' | sort | uniq > all.mutiple.samtools.snps
comm all.single.samtools.snps all.mutiple.samtools.snps | awk 'BEGIN {print "Only in single\tOnly in multiple\tIn both"; FS="\t"} {i[NF]++} END {print i[1]"\t"i[2]"\t"i[3]}'

In theory, GATK does a much better job ruling out false positives. But are there some SNPs GATK calls with high confidence that Samtools doesn't call at all?

Expand

title	Here are four commands to the answer...

Code Block

grep -v '^#' all.GATK.vcf | awk '{$2=$2"AAAAAAA"; print}' | sort -k 2 > g.v
grep -v '^#' all.samtools.vcf | awk '{$2=$2"AAAAAAA"; print}' | sort -k 2 > s.v
join -a 1 -1 2 -2 2 g.v s.v | awk '{if (NF==12) {print}}' | grep PASS > g.v.pass.only
cat g.v.pass.only | sort -k 6 -n -r | head

What's going on here:

1) Yank out all the variant calls (comments start with '#') and add a string of "AAAAAAA" to them to make "sort" do it's job in a way that "join" will later like

2) Join the two files using their chromosome position as the join field, but also include any lines from the GATK file that DON'T match the samtools file. Use "awk" to figure out which ones came only from GATK (they are missing a bunch of fields from the samtools variant calls), look only for those that GATK has labeled "PASS" and write them to a file.

3) Sort the resultant file on the variant quality value - take the top 10 lines.

You will note that many of these are complex variants, particularly insertions, so it's not too surprising that GATK does better. But here's a SNP that GATK does much better on:

chr21:34278313

It has an interesting quantitative signature though... you might want to look at it in IGV.

Other notes on bcftools

bcftools has many other sub-commands such as performing association tests and estimating allele frequency spectrums.

...

Versions Compared

Old Version 8

New Version 9

Key

Look at the differences between single- and multiple-sample SNP calls, and between Samtools/Bcftools SNP calls and GATK SNP calls

Other notes on bcftools