Daniel:Notebook/Haplotyping/MergeVCF

From ZhangLabWiki
Jump to navigation Jump to search

Creating Reference VCFs[edit]

Back to Notebook


  • The reference data used can be found in this directory on genome miner: /home/kunzhang/genomeSeq/Data/PGP1_CGI_WGS
  • The final vcfs (most current) can be found here: /oasis/tscc/scratch/djacobse/haplotypeSeqData/vcfs/merged_wgs_vcfs
    • The vcfs have names merged.chr*.hg19.liftover.reorder.vcf

Data Format[edit]

It is originally in the output format of Complete Genomics as a result of their Whole Genome Sequencing (WGS) of PGP1. There are three cell lines,

  1. PGP1FP8; Passage 8 PGP1 fibroblasts
  2. PGP1L; PGP1
  3. PGP1CDIiPS; induced Pluripotent Stem Cell line

The alignment was done using NCBI build 36 (making it several versions older than hg19, which I use). It also uses dbsnp build 130, which is also behind the times. Both of these were addressed before using it with my current data.

File Format[edit]

The data used was taken from the files var-***.tsv, which are all around 1-2GB. These files have columns that look like this:

File:Cgiwgs.example1.png

For the most part, this information is not too helpful. However, every now and then you find rows with data like this:

File:Cgiwgs.SNPexample.png

This shows two areas with SNPs. The first SNP, at position 1089, calls a SNP in both the P1 and P2 chromosome (3rd column 1 and 2, respectively). However, the SNP at 1101 calls ref for allele 1 and snp for allele 2. Note that alleles 1 and 2 are interchangable without phasing data. Therefore, the first snp called is homozygous, and not helpful for phasing, whereas the 2nd snp is heterozygous and is useful for haplotype generating.

Processing Reference Data[edit]

The reference data was processed in several steps. First, it was run through MATLAB to extract out only the heterozygous SNPs (see above for what those look like) and write them in VCF format. Step 2 lifts the VCFs over from NCBI build 36/dbsnp build 100 to current builds (hg19).

MATLAB[edit]

The MATLAB step reads in the var-*.tsv files and extracts the heterozygous SNP data, outputting as a VCF. It then takes the new VCFs (one for each chromosome/cell line) and merges them across cell lines, using only heterozygous SNPs called by at least 2 cell lines.

  1. Read in reference data and merge using 2/3 voting method (If site was called heterozygous in 2/3 cell lines)
    1. M code

LiftOver[edit]

The next step is to lift over the data from NCBI 36/dbsnp build 100 to the build used by the other data sets (hg19).

  1. Convert positions using liftOver; this converts the output vcf files to bed files, then performs the script liftOver script (available from UCSC Genome Browser)
    1. Shell script
  2. Liftover RSIDs; this is for updating dnsnp references
    1. Shell script
    2. Link to LiftRSNumber.py
  3. Reorder; I discovered later that this left a small number of SNPs out of order in the vcf file; this step is the last step I have in the current pipeline
  4. Shell Script

The final files are labeled merged.chr*.hg19.liftover.reorder.vcf, and are in the directory: /oasis/tscc/scratch/djacobse/haplotypeSeqData/vcfs/merged_wgs_vcfs