arq5x

Analysis of rare disease VCF with bcftools

Big Picture - You are the genome detective.

Here's the deal.

• We have sequenced the exomes of an affected child, and their unaffected parents. The child has a rare skin condition.
• We aligned their exome data to the human reference genome.
• We called variants using GATK.
• The resulting VCF file is called trio.trim.vep.vcf.gz.

arq5x

# compile
$ g++ zlib-example.cpp -lz -o zlib-example

# run
$ ./zlib-example
Uncompressed size is: 36
Uncompressed string is: Hello Hello Hello Hello Hello Hello!

----------

arq5x

#!/usr/bin/env python
import sys
import pprint

def make_grantham_dict(grantham_mat_file):

        """
	Citation:   http://www.ncbi.nlm.nih.gov/pubmed/4843792
	Provenance: http://www.genome.jp/dbget-bin/www_bget?aaindex:GRAR740104

arq5x

#!/usr/bin/env python

from collections import namedtuple, defaultdict
from pybedtools import BedTool
import argparse

Point    = namedtuple('Point', ['id', 'pos', 'type'])
Interval = namedtuple('Interval', ['chrom', 'start', 'end'])

arq5x

#!/usr/bin/env python
"""
given a callable file from goleft depth on stdin, merge adjacent LOW_COVERAGE and CALLABLE regions.
also split regions that are larger than max_region.
this is to even out parallelism for sending regions to freebayes.
"""
from __future__ import print_function, division
import sys
from itertools import groupby
from operator import itemgetter

arq5x

chr1  10

arq5x

GTTGTACTTCGTTCAATCGGTAGGTGTTTAACCGGATGGTCACGCCTACC
CAAGCATACTTCATTCAGTCAGGCGAAATTATTGCCAGGTCGCCGCCTAC
GTTGTACTTCGTTCAGTCGGTGGTGTTTAACTGGGTCATCGCCTACCGTG
AGTAATACTTCGTTCAGTTTGTGGAAGGTAGTGTTTAACCGGTTGCTCGC
GGTATGCGCTGGTCAAATCGGAGAGTGGGTGTTTATCGGATGGATCGCTG
CGGTGCGTGCTGTGATTCATCTTTGACTGGTGTTTATGGTCGGTCTTTAC
CATTGTACTTCCCGTTCAGTTTATCAAATTTGGTGTTTATATGAACCATA
GATGGTGGTGGTGGCGGTGGCGCTGGTGTTGCATGGTGTTGCGCATTATT
GTATCGCGTTCAGTTTCGGAAGGTGGTGTTTAACCAGTTCGCCGCCTACC
GTTTTTTCGTTCATTTGGTACGGTGTTGCGCCGGTCGCCGCCTACCGTGA

arq5x

dad	mom	kid	Inheritance description
HOM_REF	HOM_REF	HOM_REF	Expected
HOM_REF	HOM_REF	HET	Mendelian violation (plausible de novo)
HOM_REF	HOM_REF	HOM_ALT	Mendelian violation (implausible de novo)
HOM_REF	HOM_ALT	HOM_REF	Mendelian violation (uniparental disomy)
HOM_REF	HOM_ALT	HET	Expected
HOM_REF	HOM_ALT	HOM_ALT	Mendelian violation (uniparental disomy)
HOM_REF	HET	HOM_REF	Expected
HOM_REF	HET	HET	Expected

arq5x

cat a.bed
chr1	10	50	10

cat b.bed
chr1	20	40	20

cat c.bed
chr1	30	33	30

# Find the sub-intervals shared and unique to each file.

arq5x

echo "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/wgEncodeBroadHmmGm12878HMM.bed.gz
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/wgEncodeBroadHmmH1hescHMM.bed.gz
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/wgEncodeBroadHmmHepg2HMM.bed.gz
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/wgEncodeBroadHmmHmecHMM.bed.gz
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/wgEncodeBroadHmmHsmmHMM.bed.gz
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/wgEncodeBroadHmmHuvecHMM.bed.gz
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/wgEncodeBroadHmmK562HMM.bed.gz
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/wgEncodeBroadHmmNhekHMM.bed.gz
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/wgEncodeBroadHmmNhlfHMM.bed.gz" \
> chromhmm-files.txt