ChIP_Project

ChIP-Seq Pipeline

This is an example of how to construct pipelines using ngslib. The final goal is: providing a configuration file, run the ChIP-Seq basic analysis automatically.

Demo configuration file: HDAC.cfg

• genome: Bowtie indexed genome file Generated by Bowtie_build:

>bowtie_build mm9.fa mm9

In configuration file:

@genome    ~/Data/mm9/mm9

• gsizes: genome size file
  Generated by UCSC toolkit:

>faSize -detailed mm9.fa >mm9.sizes

Format:

chr1    197195432
chr2    181748087
chr3    159599783
chr4    155630120
chr5    152537259
chr6    149517037
...

In configuration file:

@gsize    ~/Data/mm9/mm9.sizes

• paired: paired end reads or not. Choice of True or False
  In configuration file:

@paired    False

• aligner: bowtie(2) or bwa, and the general parameters used
  In configuration file:

@aligner    bowtie  -m 1 -v 2 -p 6

• peakcalling: peakcalling tools, macs14 or macs2
  In configuration file:

@peakcalling    macs14 -f BAM -g mm

• bam2wig: Convert Bam file to BigWiggle file. This is a wrapper of wBamToWig.py in ngslib package.
  In configuration file: Normalized to 10M reads. Extend reads to 200bp. By default it extend it to 150bp.

@bam2wig    -n 10 -e 200

• data: Format: SRA link or file, file label, control label 1, [control label 2, ...]
  SRA link example:

/sra/sra-instant/reads/ByExp/sra/SRX/SRX122/SRX122624/SRR414932/SRR414932.sra

In configuration file:

SRR351397.sra    HDAC1_ChIP_ZT22    Control_ChIP_ZT22
SRR351398.sra    HDAC2_ChIP_ZT22    Control_ChIP_ZT22
SRR351399.sra    HDAC3_ChIP_ZT22    Control_ChIP_ZT22    IgG_ChIP_ZT22
SRR351400.sra    Control_ChIP_ZT22
SRR351401.sra    IgG_ChIP_ZT22

Real example:

###################################################
# bowtie index genome prefix
@genome ~/Data/mm9
# genome size file in format: chrN\\tlength
@gsize  ~/Data/mm9/mm9.sizes 

# paired end or not
@paired False 

# Aligner: bowtie (default) or bowtie2
# When paired end: @fastq-dump  paired
# example parameter for bowtie: -m 1 -v 2 -p 6
# example parameter for bwotie2: --fast -k 1
@aligner    bowtie  -m 1 -v 2 -p 6 

# macs version: macs14 (default) or macs2
# example parameter -f BAM -g mm
@peakcalling    macs14 -f BAM -g mm 

# bam2wig parameter: default is "-n 10 -e 200"
@bam2wig    -n 10 -e 200

# data
SRR351397.sra    HDAC1_ChIP_ZT22    Control_ChIP_ZT22
SRR351398.sra    HDAC2_ChIP_ZT22    Control_ChIP_ZT22
SRR351399.sra    HDAC3_ChIP_ZT22    Control_ChIP_ZT22    IgG_ChIP_ZT22
SRR351400.sra    Control_ChIP_ZT22
SRR351401.sra    IgG_ChIP_ZT22
###################################################

Run ChIP-Seq pipeline:

ngslib.Pipeline.ChIP_Project('HDAC.cfg')

Wrappers of existing tools

• Example of wrapping Bowtie Aligner:

class Pipeline(object):
    def bowtie(genome, fqfile,fqfile2=None, samfile=None, paras="-m 1 -v 2 -p 6",overwrite=False):
        '''
        Bowtie aligner.
        Shell example:
            bowtie ~/Data/mm9/mm9 -m 1 -v 2 -p 6 -S input.fastq output.sam
        Usage:
            Pipeline.bowtie("~/Data/mm9/mm9","input.fastq", None, "output.sam")
        '''
        sys.stderr.write("Starting bowtie ...\n")
        if not Utils.cmd_exists('bowtie'):
            raise ValueError('ERROR: bowtie cannot be found.')
        genome = os.path.expanduser(genome)
        if samfile is None:
            samfile = os.path.splitext(fqfile)[0]+".sam"
        else:
            samfile = os.path.expanduser(samfile)
        # check if samfile exists
        if os.path.isfile(samfile):
            if overwrite:
                os.remove(samfile)
            else:
                sys.stderr.write("Skipped: output file: {0} exists.\n".format(samfile))
                return
        # check if fastq file exists
        fqfile = os.path.expanduser(fqfile)
        Utils.mustexist(fqfile)
        if fqfile2 is not None:
            fqfile2 = os.path.expanduser(fqfile2)
            Utils.mustexist(fqfile2)
        # check if genome exists
        Utils.mustexist(genome+".1.ebwt")
        # run bowtie
        if fqfile2:
            cmd = "bowtie {0} {1} -1 {2} -2 {3} -S {4}".format(genome, paras, fqfile, fqfile2, samfile)
        else:
            cmd = "bowtie {0} {1} {2} -S {3}".format(genome, paras, fqfile,samfile)
        sys.stderr.write("Running command: {0}\n".format(cmd))
        p = Popen(cmd.split(), stdin=None, stdout=PIPE, stderr=PIPE)
        pstdout, pstderr = p.communicate()
        # Writing log file
        logfile = os.path.splitext(fqfile)[0]+"_bowtie.log"
        sys.stderr.write("Writing log information into {0}.\n".format(logfile))
        with open(logfile,'w') as ofh:
            print >> ofh, pstderr
    bowtie=staticmethod(bowtie)

Call Bowtie in shell:

bowtie ~/Data/mm9/mm9 -m 1 -v 2 -p 6 -S input.fastq output.sam

The equivalent code in Python:

Pipeline.bowtie("~/Data/mm9/mm9","input.fastq", None, "output.sam","-m 1 -v 2 -p 6")