SAMtools mpileup

The SAMtools mpileup utility provides a summary of the coverage of mapped reads on a reference sequence at a single base pair resolution. In addition, the output from mpileup can be piped to BCFtools to call genomic variants. I'm currently working with some Sanger sequenced PCR products, which I would like to call variants on. There are various tools for variant detection on Sanger sequences but I wanted to take this opportunity to check out SAMtools mpileup and BCFtools. In this post, I illustrate the BWA-MEM, SAMtools mpileup, and BCFtools pipeline with a bunch of randomly generated sequences.

Continue reading

VCF concordance

I want to compare the genotype concordance between two VCF files and I came across SnpSift, which seems to calculate the statistics that I want. However, the format of the results from my run differ from the format in the documentation. In this post, I will try to come up with the exact scenarios that fall into each of the summary statistics (to satisfy my curiosity).

Continue reading

Converting PED into VCF

Updated 2015 August 25th: as suggested by Tim, I checked out PLINK 1.9 and found it much simpler to convert PED into VCF. I updated the post with instructions for performing the conversion using PLINK 1.9.

Being late to the game of analysing genomic variants, I only recently discovered that IGV is capable of visualising VCF files; this is great if your variants are in the VCF. However, I have some PLINK files (a PED and MAP file), which I believe are not supported by IGV. After searching on the web, I came across a Python script written by Brad Chapman that seems to be able to convert PED files into VCF files. Since the script has several dependencies, this short post simply documents how and where these are downloaded.

Continue reading

Creating a coverage plot using BEDTools and R

One of my Top 10 posts is on creating a coverage plot using R. For that post I used CAGE data, which is a transcriptomic data set containing transcription start sites, and I used R exclusively for building a "coverage plot." The main issue with that post was that the plots were density plots rather than a real coverage plot. In this post, I'll use BEDTools to calculate the per base coverage of a defined region and produce an actual coverage plot using R.

Continue reading