Python List Comprehension

Python list comprehension enables for easy manipulation of a list variable.  For example, list comprehension can be used to convert from a string to an integer or float as well as do some common operation on a list.  Below is an example script that shows the syntax for list comprehension and how to utilize the code in Python.  Hope its helpful.

Run Unix Commands and Other Programs within Python

The python subprocess module is a part of the standard python library that enables a user to run a variety of subprocess within python and collect their output and error pipes. It is a wonderfully handy tool that really expands the possibilities of what python is capable of doing.  The tool is intended to build upon and essentially replace a lot of the functionality of the os package which I have found cumbersome to use.  The main commands of the subprocess module include:

subprocess.call-runs a command and returns output and the exit code
subprocess.check_call-runs a command; if exit code is 0 returns output, otherwise CalledProcessError

Here are some simple example scripts to help get started using the subprocess inside python code:

Creating and Accessing SQL Databases with Python

Python has numerous ways of outputting and storing data.  Recently, I investigated using shelves in Python. This is a way to store indexed data that uses command syntax similar to that of Python dictionaries, but I found it too time consuming to create shelves of large datasets.  Searching for a way to efficiently build databases in Python, I came across the SQL functionality. The library sqlite3 is an indispensable way to create databases in Python.  This package permits Python users to create and query large databases using syntax borrowed from SQL.  SQL stands for Structured Query Language and is used for managing data held in a relational database management system. The sqlite3 is a nonstandard variant of SQL query language that is compliant with the DB-API 2.0 specification. As a quick reference, I thought I would create an example script that could be used to build a SQL database using the Python programming language.  Below is a simple tutorial to follow that hopefully is useful for learning how to use the sqlite3 package.

Quick Primer on Python Shelve

Python shelve is a convenient means of storing a Python data object to disk for later use. The feature behaves similarly to a Python dictionary object and uses a lot of the same syntax.  Each object that is shelved has a key which is associated with the object so it can be quickly and efficiently accessed from disk. A shelved object must be something that can be pickled with the pickle package, essentially making shelves an easy way to organize and store pickled objects. "Pickling" is the process of making a Python object into a byte stream, and the inverse, "unpickling" is where the byte stream is restored back into the original Python object. Pickling is synonymous with serialization, marshalling, and flattening of data. Example Python objects that can be pickled include: integers, lists, tupules, sets, dictionaries, and classes. In addition to saving items to disk, shelves allow for quick access to portions of large data objects and can store them in binary format by specifying a protocol greater than 0.  Here are some simple Python scripts that hopefully serve as a useful tutorial to learn about shelves and their syntax.  The example focuses on saving dictionaries into a shelve, but it can be easily extended into other objects as well.

Easily Parse XML Files in Python

While not the most popular means of storing data, there are instances where data is stored in the XML file format.  The TCGA project is one such group that keeps a record of samples and analysis information in this format.  It can be difficult to extract data from these files without the appropriate tools.  In some cases a simple grep command can be used, but even then the output usually needs to be cleaned.  Luckily, Python has some packages that are aid in parsing XML files so that the desired data can be extracted.  The library I found most useful was ElementTree.  This package combined with the urllib2 package enables one to download an XML file remotely from the internet, parse the XML file, and extract the desired information from the data stored within the XML file.  Below is an example Python script that downloads an XML description file from TCGA (link) and then places each extracted element into a Python dictionary object.  Of course this could be easily modified for your particular application of interest, but at least it provides a simple backbone to easily build other scripts off of.

Below is the expected output:

Generate Random Genomic Positions

Generating random genomic positions or coordinates can be useful in comparing characteristics of a set of genomic loci to that what would be expected from permutations of the underlying genomic distribution.  Below is a Python script to aid in selecting random genomic positions.  The script chooses a chromosome based on probabilities assigned by chromosome length and then chooses a chromosomal position from a uniform distribution of the chromosome's length.  An added gap checking statement is included to ensure the chosen position lies within the accessible genome.  You can choose the number of positions you want, the number of permutations to conduct, the size of the genomic positions, and the genomic build of interest.  A UNIX shell script is included as a wrapper to automatically download needed chromosomal gap and cytoband files as well as run the Python script.  Useage for the UNIX script can be seen by typing ./make_random.sh from the command line after giving the script executable privileges.  An example command would be ./make_random 100 10 1000 hg19.  This command would make 10 .bed files each with 100 random 1Kb genomic regions from the hg19 genome build.  Below are the make_random.sh and make_random.py scripts.

Split a Variable in R into Components

I had a TCGA barcode that I wanted to extract information about the sample type (ie: cut out characters 14-15).  This would have been easy to do in Python (ex: id[13:15] or id.split("-")[3][0:2]), but I wanted to be able to do this inside R.  To do this I found a handy little base function called substr.  This is a function that allows you to take a subset of a string.  Here is the code to extract characters 14-15 from a string:

substr(x=id, start=14, stop=15)

or also from a variable called ids:

substr(x=data$ids, start=14, stop=15)

Sum Overlapping Base Pairs of Features from Chromosomal BED File

I had a .bed file of genomic features on a chromosome that I wanted to figure out the extent of overlap of the features to investigate commonly covered genes as well as positions where features were likely to form.  I wanted to generate a plot similar to a coverage depth plot from next-generation sequencing reads.  I am sure more efficient methods exist, but here is some Python code that takes in a .bed file of features (features.bed) and creates an output file (features.depth) with the feature overlap "depth" every 5,000 base pairs across the areas which contain features in your chromosomal .bed file.

Python Function to Calculate P-value from Z score

Every time I needed to quickly convert a z score to a p-value, I would have to search online for an online calculator or refresh my mind how to do in in R/Excel.  Getting tired of periodically having to do this, I built a simple Python function to do this for me.  Here's the code if you find it helpful:

How to Fully Utilize All Cores of a UNIX Compute Node

Parallelizing tasks can drastically improve computation time of a program.  One way to do this is to ensure your code is utilizing all available cores of a processor.  To do this you need to write your code in such a way that background tasks are being carried out simultaneously.  This is done by inserting the ampersand (&) at the end of a line of code.  If the number of background tasks running equals the number of cores of the computer node, then you are efficiently and fully utilizing the resource.  The final necessary piece of code is to use the wait command.  This tells the computer to wait until all the background tasks are completed before moving on to the next line of code.  The wait command comes in handy to ensure the number of background tasks submitted does not exceed the number of processor cores.  If this happens you will likely overwhelm the processor with too many tasks.  To prevent this the idea is to simultaneously submit a number of tasks equal to the number of cores and then use the wait command to wait for the jobs to finish before submitting more tasks.  Here is some example code I put together to utilize all 8 cores of a node when running 1,000 permutations of a program:


This will run a script with the following commands:

Of course, this could further be parallelized to run permutations simultaneously on different compute nodes as well to further speed up run time. Hope this is a helpful example to help you fully utilize compute nodes and speed up your processing time.

How to Infer Ancestry from SNP Genotypes

Self-reported ancestry is poor metric to use when attempting to statistically adjust for the effects of ancestry since several individuals falsely report their ancestry or are simply unaware of their true ancestry.  Worse yet, sometimes you don't even have information collected on an individual's ancestry.  As many of you know, if you have SNP genotyping data you can rather precisely estimate the ancestry of an individual.  Classically, to do this one needed to combine genotypes from their study sample with genotypes from a reference panel (eg: HapMap or 1000 Genomes), find the intersection of SNPs in each dataset, and then run a clustering program to see which samples clustered with the reference ancestral populations.  Not a ton of work, but a minor annoyance at best.  Luckily, a relatively new program was just released that, in essence, does a lot of this ground work for you.  It is called SNPWEIGHTS and can be downloaded here.  Essentially, the program takes SNP genotypes as input, finds the intersection of the sample genotypes with reference genotypes, weights them based on pre-configured parameters to construct the first couple of principle components (aka. eigenvectors) and then calculates an individual's percentage ancestry for each ancestral population.  Here is how to run the program.

First, make sure you have Python installed on your system and that your genotyping data is in EIGENSTRAT format.  A brief tutorial to convert to EIGENSTRAT format using the convertf tool is here.

Next, download the SNPWEIGHTS package here and a reference panel. I usually use the European, West African and East Asian ancestral populations, but there are other options on the SNPWEIGHTS webpage as well.

Then, create a parameter file with directories of input files, your input population (designated "AA", "CO", and "EA"), and a output file.  An example of one is below:


Finally, run the program using the command inferancestry.py --par par.SNPWEIGHTS.  For the program to run correctly, make sure the inferancestry.info and snpwt.co files are in the same directory as your inferancestry.py file.

For more details, see the SNPWEIGHTS paper or the README file included in the SNPWEIGHTS zip folder.  For code on generating eigenvector plots with overlaid ancestry percentages and triangle plots according to percent ancestry, see this post.

Create FASTA sequences for Phased Haplotypes

Here is some Python code I put together to convert a .haps file (and associated .sample file) into a .fasta file with an entry for each haploytpe sequence.  Haplotypes are designated >ID_A and >ID_B for each ID in the .sample file.   The program can easily be modified to accept a list of SNPs or IDs that you would like to extract from the .haps file.  Also, this program removes indels that may be present in the .haps file to avoid alignment issues.  This program was useful to feed haplotype input into phylogenetic tree programs, such as MEGA.  Just run the program by typing python make_fasta.py data.haps in at the commnad prompt and you will get a data.fasta file as output.  Hope it is useful.

Calculate Minor Allele Frequencies from VCF File Variants

Today I needed to calculate minor allele frequencies (MAFs) for sequence variants called in a .vcf file.  I couldn't find any programs that would do this for me, so I wrote a quick script to do it in Python.


This can be run in python from the command prompt by typing:


where project.vcf is the vcf file you want to calculate MAFs for.  It will return a project.txt file that contains the calculated MAF values.  This script will only work for SNPs and does not work on insertions and deletions.

Alternatively, if Python scares you there is a bit of a round about way that will do this for you too.  First use vcftools to convert your .vcf file into a Plink compatible .ped and .map file.


Then, open Plink and run the --freq option on the newly created .ped file.

**UPDATE**
Today I found an updated way to use Vcftools to directly calculate the MAF values for you.  It just takes the simple command --freq.  Here is some example code:

How To Import Command Line Arguments into a Python Script

It is often handy to be able to feed an argument from the UNIX command line to a Python script.  This is very simple to do with the sys package installed.  Below is an example of feeding three arguments (chr, start, end) from the command line into Python to be used as variables in a script.

In the script.py file:

At the command line:

Transpose List of Lists in Python

Today in Python I wanted to be able to transpose a list of lists that I created. So, for example, I wanted to be able to transpose the list

l=[[1,2,3],[4,5,6],[7,8,9]]

to this

t=[[1,4,7],[2,5,8],[3,6,9]]

I found an easy command to do this, where l is the list above and t is the transposed list.

Import List of Filenames from a Directory into Python

The python glob module is useful for importing a list of all files of a desired type from a UNIX directory.  For example, if I wanted to get a list of all .txt files in my usr/home/projects directory, I would use the following code in python.

Transform fastPHASE Output to Haploview Input

After successfully phasing with fastPHASE, I wanted to view the resulting haplotypes in Haploview.  Since I couldn't find any scripts online to convert the hapguess_switch.out file from fastPHASE into the .haps file that Haploview requires, I wrote a brief Python script to do the heavy lifting.


This script is run from the command line by typing python make_haploview_input.py hapguess_switch.out and will produce a hapguess_switch.haps file.  Because of the use of generators, this will only work on Python 2.6 or higher.  You will also need to create a .info file to match the .haps file.  The .info file has the marker loci in the same order as the .haps file and consists of two columns; the first being the marker name and the second being the marker position.