The Analytic Garden

Ebola is big news. Governors are ignoring public health officials and the White House . News media publish sensationalist headlines . From CNN we have:  "Ebola in the air? A nightmare that could happen". Is that a possibility? Is it likely? We'll get to that in a moment. The question at hand is what, if anything, does Python have to do with Ebola?  The university where I'm teaching has asked us to devote some time to discussing the Ebola epidemic, so I gave some thought about how programming can contribute to dealing with the problem.. I won't go into detailas about the epidemic here. If you want to find facts, I suggest staying away from the news and going to NIH or NIAD . The best popular account of how health officials and scientists are trying to deal with epidemic that I have encountered is from the New Yorker .  As part of the class I teach, I decided to take a brief look at whatever kind of Ebola genomic data was available and see if I could gain an...

One question of interest when dealing with paired-end reads is what was the length of the sequence that the paired mates originally came from. We mapped D. melanogaster DNA,  100 bp paired-end reads, to the genome with BWA . Among other information, BWA returns the observed template length (TLEN) for each read and its mate. With this information, producing a histogram of the fragment lengths is simply a matter of extracting the proper field from the SAM/BAM file. We used Pysam to read the BAM file. While reading we filtered read pairs and rejected any read that had a base call quality or mapping quality less than 30. We also rejected read pairs when one of the pair did not map. For the pairs that passed the filter, we simply count the number with each TLEN. We save the counts in a file for plotting with R. Here's the counting code: def ReadBamFile(bam_file): """ ReadBamFile - read a bam file and count nucleotides by context bam_file - bam file co...

The GC content of a sequence is the percentage of G's and C's in a fragment or an entire genome. Genes are often characterized by a higher GC content than background.Variations in GC content can lead to bias in read coverage in ChIP-seq and RNA-seq experiments. See Benjamini and Speed . Calculating GC content is straight forward. Select a segment of DNA or RNA sequence and count A's, C's etc. and calculate \[gc = \frac{{G + C}}{{A + C + G + T}}\]. Often GC content is calculated for an entire chromosome by calculating the GC content in a series of either overlapping or sliding windows of fixed size.  Our problem is slightly more complicated. We want to calculate the number of reads mapped to various values of GC content ranging from 0 to 100%. The first step is to calculate the GC content for a chromosome using a series of non-overlapping, fixed width windows. We can use Python string functions to do this. def GC_pct(seq_rec, window): """ ...

Probably the most common bioinformatics task is conversion between file formats. Recently, we needed to extract sequence data from a FASTQ file and write the sequences to a FASTA file. If all you want to do is extract the sequences, it's very easy in BioPython . Our task was slightly more difficult, but not by much. We also want to filter reads with low basecall quality. Our usual criteria is to reject any read with one or more bases with a quality score lower than thirty. A second, minor problem is that some of our FASQ files are in zip format, others are gzipped, and some are just text. We will use BioPython's SeqIO module to handle reading and writing the sequences. We use Python's ZipFile module to handle zipped files and the gzip module for gzipped files. The code is relatively simple. We open the FASTQ file, after first checking for the appropriate file type; read it one record at a time; check for low quality bases; if the record passes teh check, write it to ...

Ever since I started dealing with ChIP-Seq and RNA-Seq data, PCR duplicates have been a bit of a mystery to me. Googling for PCR duplicates provides a variety of opinions about how they arise, how to handle them and so forth. The best introduction I have encountered is the CureFFI.org blog. I going to swipe their description of the library prep process because it's clear and concise (they actually got it from here ): Shatter genomic DNA, e.g. with a sonicator or DNAase. Ligate adapters to both ends of the fragments. Not every molecule successfully ligates. PCR amplify the fragments with adapters Spread DNA molecules across flowcells.  The goal is to get exactly one DNA molecule per flowcell lawn of primers.  This depends purely on probability, based on the concentration of DNA. Use bridge PCR to amplify the single molecule on each bead or each lawn so that you can get a strong enough signal (whether light or pH) to detect.  Usually this requires several hundre...

A student in our group asked whether it was possible to open a .bam file located on a remote site and process it locally without having to download it. While it's possible to do that in general with a some programming trickery, unless you have some sort of special access it's usually better to download a large file and deal with it locally. If you have ssh access and reasonable connection speed to the site containing the file, you can mount the directory as a local file system and access files as if they reside on your computer. On Linux, you can use sshfs . On a Mac, there is oxfuse , an sshfs client for Mac OS-X. For Windows, there is win-sshfs . After connecting to a remote site with one of these programs, the remote directory appears as a local path. You pay a performance price for operating on files across the internet, but often convenience overrides the performance hit. Rather than directly accessing the file across the net, let's look at a simpler problem. Find...

In genomics sense refers to the direction of transcription of translation. A strand of DNA is called the sense strand if an RNA version of the same sequence is translated or translatable into protein. In the genome reference sequence genes are labeled as being on the '+' or '-' strand to indicate their sense. We want to count the number of reads mapped to genes and exons in the sense and anti-sense directions. As we did, previously, we will use an interval tree to store the gene and exon locations. Given the location, we search the gene interval tree for an overlap. If found and the direction of the read mapping matches the gene direction, we count it as being in the sense direction. We perform a similar search for the exons. One complication is that because of alternate splicing, a read may be indicated as overlapping several genes. We will need to be careful to count each read only once. def ReadBamFile(chrom, chrom_path, bam_file, gene_tree): ""...