Post-Processing Bismark Bisulphite Sequencing Data

This is Part 1 in a series on Bisulphite Sequencing.
You can skip to Part 2 (Using the Binomial Distribution in Bisulphite Sequencing).

In this post I will present a couple of scripts that are useful for processing the data produced by the Bismark bisulphite sequencing aligner into a form that is a little easier to work with.

When you run Bismark it will produce output in either its own "vanilla" format (version <= 0.5x) or in SAM format (version >= 0.6x).

Either way, the approach at this point is to remove duplicate reads that might have occurred as PCR duplicates using deduplicate_bismark_alignment_output.  I usually also remove any fragments which look like they haven't been bisulphite converted (multiple apparently methylated CHH or CHG within a single fragment).

For Bismark vanilla output you can use something like this:

awk 'FNR>1{
      countl = split($8,a,"[HX]");
      countr = split($11,a,"[HX]");
      if (countl-1 <= 3 && countr-1 <= 3)
           print;
     }' $file > $file.unconvertedreadsremoved

(Where the 8th and 11th columns contain the methylation call string in Bismark format).

Next up is to extract the methylation calls for each base using the methylation_extractor script (included with Bismark) using the --no_overlap and --comprehensive options.

This will give you a file that looks a little like this:

FCD0KMMACXX:3:1101:17551:2000#CTTCCTCC/1    +    chr3    116402006    Z
FCD0KMMACXX:3:1101:17551:2000#CTTCCTCC/1    +    chr3    116401980    Z
FCD0KMMACXX:3:1101:17551:2000#CTTCCTCC/1    +    chr3    116401960    Z
FCD0KMMACXX:3:1101:17865:1999#CTTCCTCC/1    -    chrX    115417650    z
FCD0KMMACXX:3:1101:18710:1999#CTTCCTCC/1    -    chr15    51474193    z
FCD0KMMACXX:3:1101:18710:1999#CTTCCTCC/1    +    chr15    51474352    Z
FCD0KMMACXX:3:1101:18599:2000#CTTCCTCC/1    +    chr3    101217644    Z

First we make it a little easier by sorting by chromosome and base position such that all the methylation calls for the same base are together as well as removing the Bismark header line.

awk 'FNR>1' CpG_context.txt | \
sort -k3,3 -k4,4n -S10G > CpG_context.sorted.txt

Again, as in the previous post, you can control sorts memory usage using the -S flag (set to 10G above). We can then sum up (aggregate) the file on a per CpG basis.

import getopt
import sys
import csv

def parseMethLine(line):
    (readid,meth,chrm,coord,methcall) = line
    return (True if meth == '+' else False, chrm, coord)

if __name__ == '__main__':
    try:
        opts, args = getopt.getopt(sys.argv[1:], "",
                                  ["bismark=", "aggregate="])
    except getopt.GetoptError, err:
        # print help information and exit
        # will print something like "option -a not recognized"
        print str(err) 
        sys.exit(2)
    
    infile = None
    aggregatefile= None
    
    for o, a in opts:
        if o=="--bismark":
            infile = a
            print "Bismark", a
        elif o=="--aggregate":
            aggregatefile = a
            print "Aggregate", a
    
    assert infile != None
    assert aggregatefile != None
    
    bismark = csv.reader(open(infile, "r"), delimiter="\t")
    
    aggregate = csv.writer(open(aggregatefile,"w"),delimiter='\t')      
    
    currentchrm = None
    currentpos = None
    
    methylated = 0
    unmethylated = 0
    
    for line in bismark:
        (methcall, chrm, coord) = parseMethLine(line)
        
        # meth calls in 1 base, ucsc / bedgraph in 0 base
        coord = int(coord) -1 
        
        # init
        if currentpos == None or chrm == None:
            currentchrm = chrm
            currentpos = coord
        elif currentpos != coord or currentchrm != chrm:
            # a new position or chromosome encountered
            # output our counts until now
            aggregate.writerow([currentchrm,currentpos,
                                methylated,unmethylated])
            
            # reset counts for next position
            currentchrm = chrm
            currentpos = coord
            methylated = 0
            unmethylated = 0
        
        if methcall:
            methylated += 1
        else:
            unmethylated += 1
    
    # output last line
    aggregate.writerow([currentchrm,currentpos,
                        methylated,unmethylated])

Executed as:

python aggregateBismark.py \
--bismark "CpG_context.sorted.txt" \
--aggregate "CpG_context.aggregate"

This will produce a file that's a little more useful.  Each CpG is represented by a single line with two data columns indicating the number of methylated and unmethylated calls made for that base (3rd and 4th columns respectively).

In other words we will have a file that looks like this:

chr1    1343    1    0
chr1    1625    4    0
chr1    1626    1    0
chr1    1655    1    0
chr1    1689    2    1
chr1    1690    2    0

Since we probably have more than one sample that we will want to be able to compare with each other we will want to combine two (or more) of the aggregate files above.  The script below will combine any two files with chromosome and coordinate as the first columns which are sorted the same way.

import getopt
import sys
import csv

if __name__ == '__main__':
    try:
        opts, args = getopt.getopt(sys.argv[1:], "",
                                  ["one=", "two=", "combine="])
    except getopt.GetoptError, err:
        # print help information and exit
        # will print something like "option -a not recognized"
        print str(err) 
        sys.exit(2)
    
    oneloc = None
    twoloc = None
    combineloc = None 
    
    for o, a in opts:
        if o=="--one":
            oneloc = a
            print "One", a
        elif o=="--two":
            twoloc = a
            print "Two", a
        elif o=="--combine":
            combineloc = a
            print "Combine", a
    
    assert oneloc != None
    assert twoloc != None
    assert combineloc != None
    
    left = csv.reader(open(oneloc, "r"), delimiter="\t")
    right = csv.reader(open(twoloc, "r"), delimiter="\t")
    
    combine = csv.writer(open(combineloc,"w"),delimiter='\t')      
    
    def unpackLine(chrm,coord,*data):
        return chrm.strip(), int(coord), [int(d) for d in data]
    
    def extractLine(line):
        chrm,coord,data = unpackLine(*line)
        if not chrm.startswith("chr"):
            chrm = "chr"+chrm
        return (chrm,coord,data)
    
    # get the first data points
    
    lchr,lcoord,ldata = extractLine(left.next())
    rchr,rcoord,rdata = extractLine(right.next())
    
    # find out how many data columns in each object
    # first line will tell us
    leftcols = len(ldata)
    rightcols = len(rdata)
    
    leftFinished = False
    rightFinished = False
    
    
    def spoolBoth():
        global lchr,lcoord,ldata,leftFinished
        global rchr,rcoord,rdata,rightFinished
        
        row = [lchr,lcoord]
        row.extend(ldata)
        row.extend(rdata)        
        combine.writerow(row)
        
        try:    
            lchr,lcoord,ldata = extractLine(left.next())
        except StopIteration:
            leftFinished = True    
        try:            
            rchr,rcoord,rdata = extractLine(right.next())
        except StopIteration:
            rightFinished = True 
        
    def spoolLeft():
        global lchr,lcoord,ldata,leftFinished
        
        rdata = [0 for i in range(0,rightcols)]
        
        row = [lchr,lcoord]
        row.extend(ldata)
        row.extend(rdata)        
        combine.writerow(row)
        
        try:    
            lchr,lcoord,ldata = extractLine(left.next())
        except StopIteration:
            leftFinished = True          
    
    def spoolRight():
        global rchr,rcoord,rdata,rightFinished
        
        ldata = [0 for i in range(0,leftcols)]
        
        row = [rchr,rcoord]
        row.extend(ldata)
        row.extend(rdata)
        combine.writerow(row)
        
        try:            
            rchr,rcoord,rdata = extractLine(right.next())
        except StopIteration:
            rightFinished = True    
    
    while not leftFinished or not rightFinished:
   
        if leftFinished:
            spoolRight()
        elif rightFinished:
            spoolLeft()
        else:
            if lchr == rchr and lcoord == rcoord:
                # we have a match, print it and spool both
                spoolBoth()
            elif lchr < rchr:
                spoolLeft()
            elif lchr > rchr:
                spoolRight()   
            # we are on the same chromosome
            # but we have skipped one (or more) position
            elif lcoord < rcoord:
                spoolLeft()
            else:
                spoolRight()

Execute this as follows:

python combineAggregate.py \
--one "Sample1.CpG_context.aggregate" \
--two "Sample2.CpG_context.aggregate" \
--combine "Samples.1.and.2.CpG_context.aggregate"

Ultimately you'll be left with a file where each CpG (or CHH, CHG) is a single line. A pair of columns represents the number of methylated and unmethylated calls in one sample and multiple pairs indicate multiple samples. Giving you something that looks like this.

chr1    1343    1    0    2    0
chr1    1625    4    0    5    3
chr1    1626    1    0    1    0
chr1    1655    1    0    1    5
chr1    1689    2    1    2    1
chr1    1690    2    0    3    1

I find this file format much easier to deal with for further downstream analysis. Once you've combined Sample 1 and Sample 2 you can then take the output and combine that with Sample 3 and so on.

False Discovery Rates and Large Files

Welcome to my first, of what will hopefully be many, blog post.  I'm currently a Computational Biologist at the University of Glasgow based within the Beatson Institute for Cancer Research.

While I don't officially represent either institution this blog will aim to explain some of the things that I do and to share information or code which might be useful to others working within the field.

To start with I'll present a small problem I encountered a few months ago but which I think demonstrates how powerful UNIX tools can be for tackling bioinformatics problems.

Task: I have a large file where I wish to calculate the FDR of a single column.

Problem: When the file is loaded into R in order to use the p.adjust() function the R process uses more memory than I want it to.  A file which is 6.5GB on disk; with a mixture of strings, integers and floats; peaks at about 15GB to read in using R's read table function.

I've got a pretty beefy machine for doing bioinformatic analysis on but memory is still a finite resource and I'm often working on more than one project at a time.  As such I quite often want to control the amount of memory that I use for various pipelines.

Solution: Reimplement the FDR algorithm using GNU tools.

#!/bin/sh

if [ $# -ne 2 ]
then
 echo "usage: PtoFDR file column"
 exit 1
fi

export LC_ALL=C
sort -S 10G -k"$2","$2"gr $1 | \
 awk -v col=$2 -v numrows=`wc -l $1 |  awk '{print $1}'` '
  function min(a,b)
  {
    if (a <= b)
      return a
    else
      return b
  }

  BEGIN {  cummin = 1.0; OFS="\t"; }
  {
    cummin = min(cummin,$col*(numrows/(numrows - NR + 1)));
    $col = cummin;
    print
  }'

This is just an reimplementation of the Benjamini-Hochberg procedure found in p.adjust().  This is a step-up procedure so starts from small test statistics (large p-values) and iterates through them in increasing test statistic order (decreasing p-values) -- hence the reverse sort.

You can control the amount of memory available by adjusting the -S parameter given to sort (set to 10G above).  This will limit the main memory buffer size (thus forcing sort to use temporary files on disk to do the sort).  The awk command processes the sorted file line by line so should use a minimal amount of memory.

As it's just a simple shell script it should neatly slot into many bioinformatics pipelines.

Notes:

• This will not preserve the original sort order of your input file.  If you want to do this then add a column with an incrementing counter 1..n to your file which you can then use to resort afterwards.
• The column number parameter is 1-based (i.e. 1,2,3 rather than 0,1,2).
• The P-value is replaced with the FDR(p) in the same column.  This prevents the file format from changing but if you need to keep the original p-value as well then you will need to change the script to add a new column.
• LC_ALL is set to C in the script above which has tremendous performance advantages during the sort but which limits your input files to, essentially, ANSI if you want to avoid codepoint problems.
• One alternative approach if you want to stay within an R environment is to use the mmap package to avoid loading the file into memory.

Update: Tommaso Leonardi contributed the following code which adds the FDR as an additional column instead of replacing the original p-value as my approach did.

#!/bin/sh
 
if [ $# -ne 2 ]
then
echo "usage: PtoFDR file column"
exit 1
fi

export LC_ALL=C
sort -S 10G -k"$2","$2"gr $1 | \
  awk -v col=$2 -v numrows=`wc -l $1 | awk '{print $1}'` '
  function min(a,b)
  {
  if (a <= b)
    return a
  else
    return b
  }

  BEGIN { cummin = 1.0; OFS="\t"; }
  {
    cummin = min(cummin,$col*(numrows/(numrows - NR + 1)));
    for (i = 1; i <= col ; i++){
      printf $i"\t"
    }
    printf cummin"\t";

    for (i = col+1; i <= NF; i++){
      printf $i"\t"
    }
    printf "\n"
  }'

Update 2: I've quickly made the following script which can use either of the two PtoFDR approaches above but will then resort your file back to its original order. It works by adding the row number to the beginning of every line, running an existing PtoFDR script then resorting the file by row number and removing that column. It uses a few more temporary files than is strictly necessary but this guarantees that memory usage won't ever exceed what you set for -S (otherwise if you pipe two sorts together then you can use double the memory).

#!/bin/sh

if [ $# -ne 2 ]
then
 echo "usage: PtoFDR.sorted file column"
 exit 1
fi

export LC_ALL=C
awk 'BEGIN {OFS="\t"; }{print FNR,$0}' $1 > $1.tmp
PtoFDR $1.tmp `calc $2+1` > $1.tmp2
sort -S 10G -k1,1n $1.tmp2 | cut -f2-
rm $1.tmp $1.tmp2