My first ever minimal working ChIP-seq pipeline using snakemake

Why using snakemake

Snakemake is a python3 based pipeline building tool (a python variant of GNU make) specialized for bioinformatics. I put my notes managing different versions of python here. You can write any python codes inside the Snakefile. Using snakemake is to simplify the tedious pre-processing work for large genomic data sets and for the sake of reproducibility. There are many other tools you can find here for this purpose.

Key features of snakemake

• Snakemake automatically creates missing directories.
• wildcards and Input function

To access wildcards in a shell command: {wildcards.sample}

{wildcards} is greedy (.+): {sample}.fastq could be matching sampleA.fastq if there is no sub-folder anymore, but evenwhateverfolder/sampleA.fastq can be matched as well.

One needs to think snakemake in a bottom-up way: snakemake will first look for the output files, and substitue the {wildcards} with the file names, and look for which rule can be used to creat the output, and then look for input files that are defined by the {wildcards}.

Read the following

flexible bioinformatics pipelines with snakemake
Build bioinformatics pipelines with Snakemake
snakemake ChIP-seq pipeline example
submit all the jobs immediately
snakemake-parallel-bwa
RNA-seq snakemake example
functions as inputs and derived parameters
snakemake FAQ
snakemake tutorial from the developer

examples

https://github.com/slowkow/snakefiles/blob/master/bsub.py
https://github.com/broadinstitute/viral-ngs/tree/master/pipes

A working snakemake pipeline for ChIP-seq

The folder structure is like this:

├── README.md
├── Snakemake
├── config.yaml
└── rawfastqs
    ├── sampleA
    │   ├── sampleA_L001.fastq.gz
    │   ├── sampleA_L002.fastq.gz
    │   └── sampleA_L003.fastq.gz
    ├── sampleB
    │   ├── sampleB_L001.fastq.gz
    │   ├── sampleB_L002.fastq.gz
    │   └── sampleB_L003.fastq.gz
    ├── sampleG1
    │   ├── sampleG1_L001.fastq.gz
    │   ├── sampleG1_L002.fastq.gz
    │   └── sampleG1_L003.fastq.gz
    └── sampleG2
        ├── sampleG2_L001.fastq.gz
        ├── sampleG2_L002.fastq.gz
        └── sampleG2_L003.fastq.gz

There is a folder named rawfastqs containing all the raw fastqs. each sample subfolder contains multiple fastq files from different lanes.

In this example, I have two control (Input) samples and two corresponding case(IP) samples.

CONTROLS = ["sampleG1","sampleG2"]
CASES = ["sampleA", "sampleB"]

putting them in a list inside the Snakefile. If there are many more samples, need to generate it with pythonprogrammatically.

## dry run
snakemake -np

## work flow diagram
snakemake --forceall --dag | dot -Tpng | display

To Do:

• Make the pipeline more flexiable. e.g. specify the folder name containing raw fastqs, now it is hard coded.
• write a wrapper script for submitting jobs in moab. Figuring out dependencies and --immediate-submit

merge the sub-cytobands to major cytoband and get the coordinates

In my pervious blog post, I wrote How to get the cooridnates of the cytobands

Now, how about I want to get coordinates of the bigger band merging all the sub-bands.
ex. chr1 p36.33 p36.32, p36.31, p36.23, p36.22, p36.13, p36.12, p36.11....p31.1 to p3 chr1 0 84900000
This is a task that requires some efforts for text processing. First look at what the data look like:

head -30 cytoBand.txt  
chr1    0   2300000 p36.33  gneg
chr1    2300000 5400000 p36.32  gpos25
chr1    5400000 7200000 p36.31  gneg
chr1    7200000 9200000 p36.23  gpos25
chr1    9200000 12700000    p36.22  gneg
chr1    12700000    16200000    p36.21  gpos50
chr1    16200000    20400000    p36.13  gneg
chr1    20400000    23900000    p36.12  gpos25
chr1    23900000    28000000    p36.11  gneg
chr1    28000000    30200000    p35.3   gpos25
chr1    30200000    32400000    p35.2   gneg
chr1    32400000    34600000    p35.1   gpos25
chr1    34600000    40100000    p34.3   gneg
chr1    40100000    44100000    p34.2   gpos25
chr1    44100000    46800000    p34.1   gneg
chr1    46800000    50700000    p33 gpos75
chr1    50700000    56100000    p32.3   gneg
chr1    56100000    59000000    p32.2   gpos50
chr1    59000000    61300000    p32.1   gneg
chr1    61300000    68900000    p31.3   gpos50
chr1    68900000    69700000    p31.2   gneg
chr1    69700000    84900000    p31.1   gpos100
chr1    84900000    88400000    p22.3   gneg
chr1    88400000    92000000    p22.2   gpos75
chr1    92000000    94700000    p22.1   gneg
chr1    94700000    99700000    p21.3   gpos75
chr1    99700000    102200000   p21.2   gneg
chr1    102200000   107200000   p21.1   gpos100
chr1    107200000   111800000   p13.3   gneg
chr1    111800000   116100000   p13.2   gpos50

The coordinates are sorted in a way that larger cytoband is also sorted: p3, p2, p1, q1, q2... so we only need to keep the major chromosome band, and keep expanding the end coordinates until it hits a new major band. At the same time, track which chromosome is processed and which major band is processed.

A python script using a list of lists data structure can serve this purpose:

import re
with open("/Users/Tammy/annotations/human/hg19_UCSC_genome/cytoBand.txt", "r") as f:
    # A list of lists
    # { ['chr1', 'p1','107200000', '125000000'], ['chr1',  'p2','84900000', '107200000'] ... ['chr2', 'p1' , '47800000', '93300000'],...} }
    genome_list = []
    # a set of chromosomes seen so far
    chrSet = set()
    # a set of chromosome arms seen so far
    armSet = set()
    # loop over each line
    for line in f:
        # split each line to a list
        lineSplit = line.strip().split()
        chr = lineSplit[0]
        start = lineSplit[1]
        end = lineSplit[2]
        band = lineSplit[3]
        ## this regex capture the p1 of p11, p11.1 or p12...
        arm = re.search('([pq]\d).+', band).group(1)

        # if this is the first time see this chr,
        # empty the armSet
        if chr not in chrSet:
            chrSet.add(chr)
            armSet = set()
            armSet.add(arm)
            arm_list = [chr, start, end, arm]
            genome_list.append(arm_list)
        else:
            if arm not in armSet:
                armSet.add(arm)
                arm_list = [chr, start, end, arm]
                ## append this to the genome_list
                genome_list.append(arm_list)
            else:
                ## keep the start of the arm as previous 
                new_start = genome_list[-1][1]
                ## change the previous end of the arm to current end
                new_end = end
                arm_list = [chr, new_start, new_end, arm]
                # mutate the last entry of arm_list
                genome_list[-1] = arm_list

ofile = open("/Users/Tammy/annotations/human/hg19_UCSC_genome/chrom_arm_list.txt", "w")
for arm_list in genome_list:
    ofile.write(arm_list[0] +"\t" + arm_list[1] + "\t" + arm_list[2] + "\t" + arm_list[3] + "\n")
    ## always remember to close the file!
ofile.close()

Let's look at the results.

$ head chrom_arm_list.txt 
chr1    0   84900000    p3
chr1    84900000    107200000   p2
chr1    107200000   125000000   p1
chr1    125000000   142600000   q1
chr1    142600000   185800000   q2
chr1    185800000   214500000   q3
chr1    214500000   249250621   q4
chr10   0   40200000    p1
chr10   40200000    52900000    q1
chr10   52900000    135534747   q2

There are many ways to achieve the same purpose. Just for some advanture, I used dictionary of dictionaries to record the data

import re
with open("/Users/Tammy/annotations/human/hg19_UCSC_genome/cytoBand.txt", "r") as f:
    # A dictionary of dictionaries
    # { chr1: {'p1': ['107200000', '125000000'], 'p2': ['84900000', '107200000'] ...}, chr2:{ 'p1': ['47800000', '93300000'],...} }
    genome_dict = {}
    # a set of chromosomes seen so far
    chrSet = set()
    for line in f:
        lineSplit = line.strip().split()
        chr = lineSplit[0]
        start = lineSplit[1]
        end = lineSplit[2]
        band = lineSplit[3]
        ## this regex capture the p1 of p11, p11.1 or p12
        arm = re.search('([pq]\d).+', band).group(1)

        ## if it is the first time see this chromosome
        if chr not in chrSet:
            chrSet.add(chr)
            # initiate an empty dictionary
            arm_dict = dict()
            arm_dict[arm] = [start, end]
            genome_dict[chr] = arm_dict
        else:
            # if this chromosome arm seen the first time
            if not arm_dict.get(arm):
                arm_dict[arm] = [start, end]
            else:
                new_start = arm_dict[arm][0]
                new_end = end
                arm_dict[arm] = [new_start, new_end]
            genome_dict[chr] = arm_dict

ofile = open("/Users/Tammy/annotations/human/hg19_UCSC_genome/chrom_arm_dict.txt", "w")
for key in genome_dict.keys():
    for key2 in genome_dict[key].keys():
        ofile.write(key + "\t" + genome_dict[key][key2][0] + "\t" + genome_dict[key][key2][1] + "\t" + key2 +"\n")
ofile.close()

Note that, There is no order of dictionary, so the output may not be sorted. You can sort by chr and start:
cat chrom_arm_dict.txt | sort -k1,1V -k2,2n

V is only for GNU sort which will sort chromsome in alpha-numeric order

If I compare the two results from different ways, they are identical:
cmp <(sort -k1,1V -k2,2n chrom_arm_dict.txt) <(sort -k1,1V -k2,2n chrom_arm_list.txt)

Test significance of overlapping genes from 2 gene sets

The question is very straightforward:
I have three microarray or RNA-seq data sets (control, treatment1, treatment2, each should have several replicates, otherwise you can not get a pvalue for differential expression!) , I then get two sets of differentially up-regulated genes ( treatment1 vs control, treatment2 vs control).  There are some overlapping genes which are both up-regulated by treatment1 and treatment2. How significant the overlapping is?

The easiest way is to use an online program http://nemates.org/MA/progs/overlap_stats.html
The program calculates the p value by using Hypergeometric_distribution.

see below for a python and an R script for this problem.

read  posts here http://stats.stackexchange.com/questions/16247/calculating-the-probability-of-gene-list-overlap-between-an-rna-seq-and-a-chip-c
http://blog.nextgenetics.net/?e=94
https://www.biostars.org/p/90662/

list1=3000, list2=400, Total gene number (population)=15000,
overlapping between list1 and list2 =100
From the stackexchange link:
"The p-value you want is the probability of getting 100 or more white balls in a sample of size 400 from an urn with 3000 white balls and 12000 black balls. Here are four ways to calculate it.

sum(dhyper(100:400, 3000, 12000, 400))
1 - sum(dhyper(0:99, 3000, 12000, 400))
phyper(99, 3000, 12000, 400, lower.tail=FALSE)
1-phyper(99, 3000, 12000, 400)

These give 0.0078.

dhyper(x, m, n, k) gives the probability of drawing exactly x. In the first line, we sum up the probabilities for 100 – 400; in the second line, we take 1 minus the sum of the probabilities of 0 – 99.

phyper(x, m, n, k) gives the probability of getting x or fewer, so phyper(x, m, n, k) is the same as sum(dhyper(0:x, m, n, k)).

The lower.tail=FALSE is a bit confusing. phyper(x, m, n, k, lower.tail=FALSE) is the same as 1-phyper(x, m, n, k), and so is the probability of x+1 or more "

If I run with the python script, I got the same p-value:

tommy@tommy-ThinkPad-T420 ~/scripts_general_use
$ ./hyper_geometric.py 15000 3000 400 100

total number in population: 15000
total number with condition in population: 3000
number in subset: 400
number with condition in subset: 100

p-value <= 100: 0.996049535673

p-value >= 100: 0.00784035207977

The bottom line is that we need more statistics!!

mapping gene ids with mygene

Mapping gene ids is one of the routine jobs for bioinformatics. I was aware of several ways to do it including Biomart.

Update on 10/30/14, a mygene bioconductor package is online http://bioconductor.org/packages/release/bioc/html/mygene.html

Recently I got to know mygene, a python wrapper for the mygene.info services to map gene ids.
I found it very handy to convert gene ids.  see a gist below.

To use it,  cat input.txt | python geneSymbol2Entrez.py > output.txt
or python geneSymbol2Entrez.py input.txt > output.txt  where input.txt contains one gene name in each line. pretty neat!

use python to change the header of a fasta file based on a dictionary in another file

I saw on SEQanwser, someone asked this question here, and I used python to resolve this problem.
I have not written in python for so long, I've been using R for a long time. This provide me a good practice opportunity.

The basic problem is to format some files, and I know it can be done by other tools like sed and awk ( see an awk solution below), but python just gives you much more flexibility in terms of what you can do and gives you more readability.
myfasta file

>comp1558_c0_seq1 len=204 path=[4976:0-203]
ATTGTACTTAATCTA.....
>comp8142_c0_seq1 len=222 path=[8835:0-221]
GTACAATCACGGCACATT......
>comp8570_c0_seq1 len=232 path=[3344:0-232]
GCATCGCTTATCGGTTTA......

annotation file:

comp1558_c0_seq1 repressor protein
comp8142_c0_seq1 aspartate aminotransferase
comp8357_c0_seq1 l-xylulose reductase
comp8387_c0_seq1 succinyl- synthetase alpha
comp8570_c0_seq1 fatty acid synthase beta subunit dehydratase

see the code below:




output:
>comp1558_c0_seq1 repressor protein
ATTGTACTTAATCTA.....
>comp8142_c0_seq1 aspartate aminotransferase
GTACAATCACGGCACATT......
>comp8570_c0_seq1 fatty acid synthase beta subunit dehydratase
GCATCGCTTATCGGTTTA......

if the dictionary file has only two records for each line, this awk should work:
http://stackoverflow.com/questions/11678939/replace-text-based-on-a-dictionary

Usage: awk -f foo.awk dict.dat user.dat
http://www.gnu.org/software/gawk/manual/html_node/String-Functions.html
http://www.gnu.org/software/gawk/manual/html_node/Arrays.html

NR == FNR {
  rep[$1] = $2
  next
} 

{
    for (key in rep) {
      gsub(key, rep[key])
    }
    print
}

interactive visualization packages for R and python

I've been collecting potential packages for (interactive) visualization of Genomic data.
Here are the ones so far I've collected:

python:
Vispy http://vispy.org/gallery.html
pretty plotlib https://github.com/olgabot/prettyplotlib
D3py https://github.com/mikedewar/d3py
Yhat ggplot for python http://blog.yhathq.com/posts/ggplot-for-python.html
galry https://github.com/rossant/galry
Glue http://www.glueviz.org/en/latest/
Seaborn https://github.com/mwaskom/seaborn
Bokeh https://github.com/ContinuumIO/Bokeh

R:
Click me https://github.com/nachocab/clickme
plotly https://plot.ly/

others:
http://selection.datavisualization.ch/
Gpviz http://icbi.at/software/gpviz/gpviz.shtml
Jheatmap http://jheatmap.github.io/jheatmap/
Caleydo http://caleydo.github.io/caleydo-doc/3.1/index.html#!index.md

There are too many out there, I will need to focus on  one or two.

python code for getting the reverse complement DNA strand

Following pythonforbiologist http://pythonforbiologists.com/index.php/business-card-solutions/, I wrote a very simple python script to get the reverse complement DNA strand:


It is a very simple script but several things I want to mention:
1. what is the fast string concatenation method in python?
in the first function, we tend to write something like:

out_seq = ' '
for base in reversed(seq):
    out_seq += seq_dict[base]

but it turned out to be very slow. Read this post for more details  http://www.skymind.com/~ocrow/python_string/
the list comprehension is the fastest and most elegant way to concatenate a string.


2. for more complicated IUPAC ambiguity codes, read this post:
http://pythonforbiologists.com/index.php/how-to-count-non-dna-bases-in-a-sequence-using-python/

Nucleotide Code:  Base:
----------------  -----
A.................Adenine
C.................Cytosine
G.................Guanine
T (or U)..........Thymine (or Uracil)
R.................A or G
Y.................C or T
S.................G or C
W.................A or T
K.................G or T
M.................A or C
B.................C or G or T
D.................A or G or T
H.................A or C or T
V.................A or C or G
N.................any base
. or -............gap

group the same ID to the same line by python, use awk to do SQL like groupby and sum

I was on Seqanswer forum and found someone asked a question:
http://seqanswers.com/forums/showthread.php?p=119263&posted=1#post119263
"Hi guys,
i hv got a GO file for my differentially expressed genes file, it goes like:

FBgn00001 GO:0016301 [Name:****(annotation)]
FBgn00002 GO:0016301 [Name:****(annotation)]
FBgn00003 GO:0016301 [Name:****(annotation)]
FBgn00004 GO:0003700 [Name:****(annotation)]
FBgn00004 GO:0009651 [Name:****(annotation)]
FBgn00004 GO:0006355 [Name:****(annotation)]
FBgn00005 GO:0009556 [Name:****(annotation)]
FBgn00005 GO:0005515 [Name:****(annotation)]
FBgn00005 GO:0080019 [Name:****(annotation)]
FBgn00005 GO:0016563 [Name:****(annotation)]
FBgn00005 GO:0016627 [Name:****(annotation)]
FBgn00006 GO:0003700 [Name:****(annotation)]
FBgn00006 GO:0010018 [Name:****(annotation)]

now i want to use WEGO ,so i need to convert it like:

FBgn00001 GO:0016301
FBgn00002 GO:0016301
FBgn00003 GO:0016301
FBgn00004 GO:0003700 GO:0009651 GO:0006355
FBgn00005 GO:0009556 GO:0005515 GO:0080019 GO:0016563 GO:0016627
FBgn00006 GO:0003700 GO:0010018

I think this could be solved using a perl script. I am not able to do this since i am a beginner. Can someone help me out? A simple perl script is good enough for me^^ "


The question is on very basic text manipulation, and I quickly wrote a python script for that:

import csv
reader = csv.reader(open("GO.txt","r"), delimiter="\t")
new={}
for row in reader:
    if row[0] not in new.keys():
        new[row[0]] = [row[1]]
    else:
        new[row[0]].append(row[1])


with open("wego.txt","w") as f:
    for key, value in sorted(new.items()):
       f.write(key+"\t"+"\t".join(value)+"\n")

update 10/19/13, I used the gist in github to display the code

basically, I store the ID as the key, and a list of GO terms as the value, and then write the dictionary
to a txt file.
I was thinking if there are any easier ways to do it, I mean by awk or sed.
A quick google search
http://stackoverflow.com/questions/10286522/group-by-sum-from-shell/10289422#10289422

the question is a little bit different from the one above:

"I have a large file containing data like this:

a 23
b 8
a 22
b 1

I want to be able to get this:

a 45
b 9

"

awk array is very handy for that purpose see a post here http://www.thegeekstuff.com/2010/03/awk-arrays-explained-with-5-practical-examples/:
let's first creat the file:
tommy@tommy-ThinkPad-T420:~$ cat > foo.txt
a 23
b 8
a 22
b 1
#control+D to indicate cat this is the end of the file.


tommy@tommy-ThinkPad-T420:~$ cat foo.txt | awk '{new[$1]+=$2}END{for (key in new) print key, new[key]}'
a 45
b 9

#or you can redirect the result to a new file


tommy@tommy-ThinkPad-T420:~$ cat foo.txt | awk '{new[$1]+=$2}END{for (key in new) print key, new[key]}' > new.txt

I do not know whether the array value can hold multiple strings like a list in python.
-------------------------------------------
updated on 10/21/13,  an awk one liner can solve this problem

tommy@tommy-ThinkPad-T420:~$ cat foo.txt | awk '{ if (new[$1]) new[$1]=new[$1]"\t"$2; else new[$1]=$2;} END { for (i in new) print i, new[i]}' OFS="\t"
a 23 22
b 8 1

Amazing awk!

Decorators in python

To help you understand Decorators in python, read the posts below:

http://simeonfranklin.com/blog/2012/jul/1/python-decorators-in-12-steps/

http://www.artima.com/weblogs/viewpost.jsp?thread=240808

and here:
http://pythonarticles.com/decorators_part1.html
http://pythonarticles.com/decorators_part2.html

Enjoy!

nice demonstration of *args and **kwargs in python function def

see a post here
http://freepythontips.wordpress.com/2013/08/04/args-and-kwargs-in-python-explained/

python for biologist

I just found this website:
http://pythonforbiologists.com/index.php/welcome-to-python-for-biologists/

it teaches you how to start programming with python for biologists.  I even bought the pdf ebook for $39. I had a look at the book, and I think it is worthy. I went through first 4 chapters really fast, they are really basic but teach biologists in a really friendly way. I am even thinking to hold a bioinformatics training course in the future :)  sounds ambitious, hmm..

Advanced topics can be found in the website:

Here is a link to a .zip file containing solutions to some of the exercises.

Section 1 : Large programs and complex data structures

Section 2 : Retrieving and sending information over the web

Section 3 : Large datasets and treelike data structures

Section 4 : Classes, modules and web applications

I've read the practical computing for biologists by Haddock Dunn http://www.amazon.com/Practical-Computing-Biologists-Steven-Haddock/dp/0878933913
it is a very good introductory book. Besides python, it also teaches you how to use linux command lines and some database basics. By the way, I wrote my first regular expression after reading the book.

well, I feel I just could get some practical problems ( text manipulation, Genomic Interval calculation using pybedtools, NGS by HTSeq) done by python though I began to play around with python in 04/24/2012. It takes time, but the investment definitely is rewarding!!