The Analytic Garden

The other day I showed how to filter a bam file to remove low quality reads . As a side effect of this, we also removed all reads that did not map to the particular chromosome that the bam file referenced. This greatly reduced the size of the bam file. In addition, it makes the next task much easier. We want to know about the read coverage in each of the chromosomes to visualize how the coverage might differ among runs and sources of RNA and DNA. Coverage is the number of reads overlapping a particular genomic position. Counting coverage is pretty easy with the PySam module. Here's one way to do it: def ReadBamFile(chrom, bam_file, genome_file):     """     ReadBamFile - read a bam file and count nucleotides by context     chrom - the chromosome name as aa string, e.g. chr2L     bam_file - bam file containing reads     genome_file - file containing the genome in fasta format     returns: a nump...

Recently, we wanted to remove low quality and duplicate reads from .bam alignment files to make viewing them in the UCSC genome browser easier. Our definition of low quality is stringent. We consider a read as low quality if it has any base with a base call quality below 30. To filter low quality reads and duplicates, we use a combination of samtools and the python PySam library. samtools has an option to remove PCR and optical duplicates and we use PySam to filter low quality reads. PySam allows us to read and write bam files. To filter out low quality reads, we loop through the input bam file, check the read quality to see if any base call falls below 30. If all of the base call qualities are 30 or above, we write the read to a new file. Here's one way to filter: def FilterReads(chrom, in_file, out_file):     def read_ok(read):         """         read_ok - reject reads with a low quality (be...

Sequence data is often stored in Fasta format . Fasta files are text files, so to access them, you can use standard I/O routines to read and write fasta files. However, libraries like BioPerl and BioPython provide a convenient way to manipulate these files. Fasta files often contain multiple sequences. Here's a simple example of using BioPython to read a Fasta file and return a list of SeqRecord objects. SeqRecord objects contain sequence data, the sequence ID, possibly annotation and other information. It's a uniform object for storing sequence data. def ReadFasta(file):     """     ReadFasta - read a sequence from a fasta formatted file.     file - fasta file name     returns: a list of BioPython sequence record     """     seq_list = []     handle = open(file, 'rU')     for record in  SeqIO.parse(handle, 'fasta'):        ...

One of the most basic things you want to know when aligning sequence reads to a genome is how many reads cover a given position. Here we'll take an alignment from reads produced by an Illumina HiSeq 2500 sequence and aligned with BWA .  To count, we'll make use of PySam 's pileup routine. The process is pretty simple. We have pileup iterate through each position in the genome that is covered by a read. Pileup returns a list of read overlapping the position. We through away any read that contains a base call quality below 30 and increment a count array for each read that overlaps a given position. Here's the relevant code    for pileup_col in f.pileup(id, max_depth=12000):         for pileup_read in pileup_col.pileups:             if not pileup_read.indel or  pileup_read.is_del:                 aln =...

I put up a custom track on the UCSC genome browser for the BWA alignment of a small portion of the reads from INTACT DNA run YAS11_CGATGT_L006_R1_001 for chr2L. I wanted to let everyone take a look at it and see if you want more tracks. Here's how to view the track.  Point your browser to  https://genome.ucsc.edu/cgi-bin/hgCustom Select the insect clade, the D. melanogaster genome and the Apr 2006 assembly Copy the following text and paste it into the text box: track type=bam name="INTACT DNA YAS11" bigDataUrl=http://ccmbweb.ccv.brown.edu/bam/YAS11_CGATGT_L006_R1_001/chr2L.extract-sorted.bam You should see something like this Click the submit button. You should see something like  Click go to genome browser The genome browser may not take you to chr2L or it may take you to a region without reads. Paste chr2L:826,001-851,000 into the location box. You should see something like the following, but the details will vary ...

I described a method for  counting mismatches previously. This time, I'll count something different in a slightly different manner, but the principles will be the same. In this case, we're looking for mismatches to the official genome in sequence reads, but we want to know how the mismatch rate varies along the read. The position along the read is called the cycle , because it indicates the machine cycle that sequenced the particular base. The basic method is similar to the one used in the previous post, but this time we'll use perl instead of python. Perl receives a lot of bad press these days, having gone from the being the Swiss-army chainsaw of the web to having its web processing thunder stolen by Php and its general programming ubiquity being usurped by python. Still, I like perl and use it frequently. Perhaps, this is because I have been using it for so long that I know many of its idiosyncrasies and can either workaround them or use them to my advantage. Despite...

Counting is hard. It seems like it shouldn't be. After all, we all learned to count as part of our first experiences in school, or maybe even before. In computational biology, we often have to count to determine basic properties of a collection of sequence or to estimate probabilities and expected values. The problem we frequently run into is how do we determine what to count and how do we go about counting it. Sequencing Errors One of the questions we have been interested in answering is; what are the basic error rates in the base calls coming from the genome center's Illumina sequencing machine. A brief overview of the sequencing and alignment processes can be downloaded here . If you looked at https://www.youtube.com/watch?v=l99aKKHcxC4 you can see how reads are built one nucleotide at a time. The base calling process determines the nucleotide type by image processing of  emitted fluorescence. For each nucleotide, intensity values are calculated for each type: A, C, G, a...