The End of 2017

In the end of last year, I wrote a post summarizing the passing 2016. Now, it is time to write the same for 2017! How time flies!

Last year, I wrote:

For the coming 2017, I should be:
1. busy with Phoebe.
2. writing 1-2 papers.
3. writing a book chapter on ChIP-seq for the biostar handbook. It should come out in the mid of 2017.
3. writing a small R package for practice.
4. learning a piano song.

Looking back, it seems I only accomplished 1 and 3 :)  I do have two first-author papers in writing-stage, but I have not finished them yet. I wish I  could get them out in 2018.

The book chapter on ChIP-seq was published here. If you want a PDF of my chapter, you can download from https://drive.google.com/open?id=1UxB0uhsoWlPvymukP3em8v4vBDI-CKwK

I still have not got a chance to write an R package, which is on my list for long. The coming 2018 is a good time for me to get my hands wet. Our epigenomic project was selected by the  Data Science Road-Trip program !! I received the confirmation in the end of 2017. I look forward to learn more R and machine learning for 2 weeks. And the plan is to turn the work into an R package. Many thanks to @JasonWilliamsNY, I saw this opportunity from his tweet. The development of R package is becoming easier with hadley wickham's work usethis and many others.

I failed #4 totally...I did not get any time to practice the piano. My wife is occupied by Phoebe and does not have time to teach me either.

I have some other achievements that I think I need to celebrate though:)

I am proud to say that I am an official instructor for Data Carpentry! I went to the instructor training in UC Davis hosted by Dr.Titus Brown. I am excited to be part of this welcoming community.

I am also excited to be a  Genomics Advisory committee member of Data Carpentry and a maintainer of the Wrangling Genomics course materials.

I got to practice some teaching skills learned from the instructor training in a ChIP-seq lab, which is part of the GS01 1143 Introduction to Bioinformatics course.

Together with Dr.Kunal Rai, I authored a book chapter on Computational Analysis of Epigenetic Modifications in Melanoma. It should be out in 2018.

The other thing I want to mention is that I wrote several NGS processing workflows using Snakemake, a python extension. It was my first time to write something seriously using python and I like it a lot.

The most complex workflow I have written is the mutation calling pipeline (find it here).
I follow GATK best practices and incorporate two somatic callers: mutect and lancet. In my current lab of Dr.Jianjun Zhang and Dr.Andrew Futreal, I deal with multi-region or longitudinal tumor samples from the same patient. In the pipeline, I implemented copy-number analysis and pyclone clonal architecture analysis. The Snakefile is over 1000 lines, I need to think about modularizing it.

Of course, my ChIP-seq snakemake processing pipeline is used to uniformly process thousands of data sets generated in Dr.Kunal Rai's lab. I am happy that many wet lab members are learning computation and are using it by themselves.

In addition to DNA-seq and ChIP-seq. I developed several other workflows:
RNA-seq
ATAC-seq
RRBSeq

I will need to better document all the workflows.

For the coming 2018:

1. I will need to get out at least 2 papers. I enjoyed the public service such as involving Data Carpentry, but I know if I do not have papers, my academic future is doomed.

2. That does not mean I will spend less time on teaching. In fact, I plan to bring at least one Data Carpentry workshop down to Genomic Medicine department in MD Anderson Cancer Center.

3. Finish that R package of course!

4.  I am expecting my second child Noah in April 2018 ( We survived Hurricane Harvey 2017!!).  I know it will be even busier with 2 kids :) I love my first kid Phoebe, she is now 17 months. The joy she has brought to me is not exchangeable with anything else. I endeavor to be a better researcher, but first I need to be a better husband and father.

Looking forward to a brand new 2018! Good luck everyone.

linux commands basics

I am attending the NGS course at MSU. This is a great course with great instructors and friendly colleagues.
I highly recommend this course to everyone. http://bioinformatics.msu.edu/ngs-summer-course-2014

This morning, we learned SNP calling by samtools and sam file specification (I will write another blog for the SNP calling) .in the night , TA Elijah gave an awesome introduction to linux commands.
personally, I think this should be taught in the first day of the course. ( I am already pretty familiar with basic linux commands, but it does cause a lot of frustrations for beginners).

I took the notes, and put the commands that taught in a gist, see below and enjoy linux commands!


linux basics by Tommy Tang is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Several NGS bioinformatics training materials

1.  Next Generation sequencing wiki book http://en.wikibooks.org/wiki/Next_Generation_Sequencing_(NGS)

2. UCDavis training course https://training.bioinformatics.ucdavis.edu/documentation/

3. MSU training course by Titus Brown http://ged.msu.edu/angus/index.html

4. GOBLET http://www.mygoblet.org/training-portal

5.  more http://ged.msu.edu/angus/bioinformatics-courses.html

6. a list maintained by Stephen Turner in UVA http://stephenturner.us/edu.html

7. GVL practial protocols for ChIP-seq, RNA-seq, variant calling etc https://genome.edu.au/wiki/GVL

These will keep me busy for a while...I feel I do not have enough time to learn :)

ChIP-seq analysis programs

This is the first blog in the year of 2014.  I am graduating in the coming August, there are so many things to do and I do not have much time to blog. However, I do not want it to be inactive.

I have extensive experience with ChIP-seq data analysis, and I want to list all the programs that I usually use.

With raw Fastq data, one needs to align it to the reference genome by an aligner. I usually use Bowtie http://bowtie-bio.sourceforge.net/index.shtml, others like BWA (developed by Heng Li, a legendary bioinformatician), MAQ, SOAP etc are also very popular (http://massgenomics.org/short-read-aligners)

Many times, I download data from SRA, sra toolkit http://eutils.ncbi.nih.gov/Traces/sra/?view=software is used to convert the sra format to fastq, and then one can align it with bowtie.

After mapping with Bowtie, one get a sam file http://genome.sph.umich.edu/wiki/SAM. samtools (also developed by Heng Li) is used to convert it to bam file, a binary form of sam file.

MACS https://github.com/taoliu/MACS/ is the most widely used peak calling program developed by Tao Li, previously in shirley Liu's lab at Harvard. MACS2 now can detect differentially changed peaks (see the link above). It also generates a bedgraph file to be visualized by IGV http://www.broadinstitute.org/igv/.

For motif enrichment analysis:
Homer  http://homer.salk.edu/homer/ngs/index.html
pscan-ChIP http://159.149.160.51/pscan_chip_dev/
MEME http://meme.nbcr.net/meme/
GREAT for GO analysis http://bejerano.stanford.edu/great/public/html/
Homer, CEAS http://liulab.dfci.harvard.edu/CEAS/ or PAVIS http://manticore.niehs.nih.gov:8080/pavis/ to annotate peaks

you might be interested in my previous blogs also:
http://crazyhottommy.blogspot.com/2013/12/chip-seq-peaks-overlapping-significance.html
http://crazyhottommy.blogspot.com/2013/10/compare-chip-seq-data-for-different.html
http://crazyhottommy.blogspot.com/2013/08/how-to-make-heatmap-based-on-chip-seq.html
http://crazyhottommy.blogspot.com/2013/04/how-to-make-tss-plot-using-rna-seq-and.html

2014 will be an exciting year for me! Follow me and many others I follow on twitter https://twitter.com/tangming2005.

ChIP-exo data analysis

Our neighboring lab just generated some ChIP-exo data, (if you do not know the technique, look at this paper http://www.ncbi.nlm.nih.gov/pubmed/23026909) and the company did some analysis but not what they want. The bedgraph files generated by the company were separated for plus strand and minus strand. They want just peak files like in the ChIP-seq experiment.

I was aware of several software can be used to deal with this type of data:
Tao Liu's MACS https://github.com/taoliu/MACS/
this is probably the most widely used method for ChIP-seq peak calling, and it is also suitable for ChIP-exo data
https://github.com/taoliu/MACS/issues/15

Peakzilla https://github.com/steinmann/peakzilla
A new tool for ChIP-exo

"Peakzilla identifies sites of enrichment and transcription factor binding sites from transcription factor ChIP-seq and ChIP-exo experiments at hight accuracy and resolution. It is designed to perform equally well for data from any species. All necessary parameters are estimated from the data. Peakzilla is suitable for both single and paired end data from any sequencing platform."

A quick google I found several others:
MACE http://chipexo.sourceforge.net/
and GEM http://www.psrg.csail.mit.edu/gem/

Since MACS (version 1.4) was pre-installed in the HPC at UFL, I decided to give it a try

[mtang@dev1 mm10]$ module load macs

macs -t ChIP.bam -c control.bam -f BAM -g mm -n ChIP-exo -B

it took some time to finish the process (~30mins).  with -B flag, I want to output the bedgraph files for each chromosomes. If you specify -S, it will give you a single bedgraph file for the whole genome.

Look at the model built by MACS

You can have a quick look at the data by loading the bedgraph files into IGV. I just pick VEGFa to check

The peaks look very specific and sharp.

Then, I can find all the genes that contain a peak nearby.
Homer annotatePeaks http://biowhat.ucsd.edu/homer/ngs/annotation.html
if you use R, ChIPpeakAnno http://www.ncbi.nlm.nih.gov/pubmed/20459804
cistrome can do it very easily http://cistrome.org

BETA-minus: Targets prediction with binding only Predict the factors (TFs or CRs) direct target genes by only binding data
CEAS http://liulab.dfci.harvard.edu/CEAS/  also from Liu's lab
the easiest way PAVIS http://manticore.niehs.nih.gov:8080/pavis/

Two pipelines for NGS analysis

bcbio-nextgen from Brad Chapman

https://github.com/chapmanb/bcbio-nextgen
Documentation https://bcbio-nextgen.readthedocs.org/en/latest/index.html

Metaseq from Ryan Dale

https://github.com/daler/metaseq

Documentation http://pythonhosted.org/metaseq/

DESeq2 is now in Bioconductor

http://seqanswers.com/forums/showthread.php?t=28350
http://bioconductor.org/packages/2.12/bioc/html/DESeq2.html

Thank you Simon and your team!

install bioawk in unbuntu

Bioawk is written by Heng Li, and it can handle formats like:

bed:
     1:chrom 2:start 3:end 4:name 5:score 6:strand 7:thickstart 8:thickend 9:rgb 10:blockcount 11:blocksizes 12:blockstarts
sam:
    1:qname 2:flag 3:rname 4:pos 5:mapq 6:cigar 7:rnext 8:pnext 9:tlen 10:seq 11:qual
vcf:
    1:chrom 2:pos 3:id 4:ref 5:alt 6:qual 7:filter 8:info
gff:
    1:seqname 2:source 3:feature 4:start 5:end 6:score 7:filter 8:strand 9:group 10:attribute
fastx: 1:name 2:seq 3:qual 4:comment

 I followed the tutorial here https://github.com/vsbuffalo/bioawk-tutorial

my first try did not work

tommy@tommy-ThinkPad-T420:~$ git clone git://github.com/lh3/bioawk.git && cd bioawk && make && mv awk bioawk && sudo cp bioawk /usr/local/bin/
Cloning into 'bioawk'...
remote: Counting objects: 163, done.
remote: Compressing objects: 100% (89/89), done.
remote: Total 163 (delta 95), reused 136 (delta 74)
Receiving objects: 100% (163/163), 112.32 KiB, done.
Resolving deltas: 100% (95/95), done.
yacc -d awkgram.y
make: yacc: Command not found
make: *** [ytab.o] Error 127

It looks like I do not have  yacc or bison (the GNU
equivalent) installed.

tommy@tommy-ThinkPad-T420:~/bioawk$ sudo synaptic

search yacc, and install the bison.

after that, it worked.

tommy@tommy-ThinkPad-T420:~/bioawk$ bioawk
usage: bioawk [-F fs] [-v var=value] [-c fmt] [-H] [-f progfile | 'prog'] [file ...]

a quick tutorial for git:
https://github.com/vsbuffalo/git-demo

a python script for extracting fastq sequences

a simple python script to handle fastq file.



fast fastq file parser, extract ~10,000 reads with their sequences and quality 
from a fastq file (~70milion reads) based on names  of the reads. HengLi in Princeton https://github.com/lh3/seqtk
has a wrapper for this kind of task.
the following code is from http://www.biostars.org/p/10353/  
this script demonstrates the usage of set(hash-able), a data structure that is much faster than list
when you have two files, put the information from one file into a container, loop over the other file.
usage: cat file.fastq | python parse_fastq.py id_file.txt > selected.fastq


#! /usr/bin/env python
import sys
# get filename from parameter
idfile = sys.argv[1]

# load ids in a set with  a set comprehension 

ids = set( x.strip() for x in open(idfile) )

# read the fastq file
handle = sys.stdin

while ids:
    #parse fastq
    idline = handle.readline()
    seq   = handle.readline()
    spacer = handle.readline()
    quals = handle.readline()

    id_name = idline[:-1] # except the newline \n
    if id_name in ids:
        #print fastq
        sys.stdout.write( '%s%s%s%s%' % ( idline, seq, spacer,\ quals) )
        ids.remove(id_name)



see link here https://github.com/lh3/seqtk

The NGS WikiBook: a dynamic collaborative online training effort with long-term sustainability

The NGS wiki book http://en.wikibooks.org/wiki/Next_Generation_Sequencing_(NGS)

http://intl-bib.oxfordjournals.org/content/early/2013/06/21/bib.bbt045.full

for more details, refer to the link above.
Nine simple rules to begin with NGS analysis

RULE 1: DO NOT FEAR THE COMMAND LINE

RULE 2: KNOW THE CONVENTIONS

RULE 3: READ INTRODUCTORY REVIEWS

RULE 4: START WITH QUALITY CHECKING

RULE 5: PLAN FOR MISTAKES AND DOCUMENT WORKFLOW

RULE 6: ALWAYS GET INFORMED AND GET HELP IF STUCK

RULE 7: USE AN EFFICIENT INTEGRATIVE APPROACH

RULE 8: AVOID REINVENTING THE WHEEL  ( I mentioned this several times in my blog :) )

RULE 9: EDUCATION IS IMPORTANT

VCF file annotation and manipulation tools

VCF (variant calling format) file, as specified in 1000 Genome http://www.1000genomes.org/wiki/Analysis/Variant%20Call%20Format/vcf-variant-call-format-version-41
is a common file to handle the variant information from the NGS sequencing data.

many tools have been developed to manipulate, extract information with this format.
some popular ones are:
tabix http://samtools.sourceforge.net/tabix.shtml for fast indexing and extracting certain regions of the whole genome.
VCFtools http://vcftools.sourceforge.net/ for many more manipulations including annotating, merging, concatenating, comparing etc.
bedtools http://code.google.com/p/bedtools/ for genomic region overlapping calculation. etc
I have some experience with the above three tools.

by the way, I have not used the GATK tools http://www.broadinstitute.org/gatk/from the broad Institute yet, I am sure I will have a try sometime later.

after google, I found several more:
variationtoolkit http://code.google.com/p/variationtoolkit/
vcflib https://github.com/ekg/vcflib
varianttools http://varianttools.sourceforge.net/
variantannotation http://www.bioconductor.org/packages/2.12/bioc/html/VariantAnnotation.html
Plinkseq http://atgu.mgh.harvard.edu/plinkseq/

Taser is the one I just came across http://www.zhanxw.com/taser/
It is based on R and can extract variant info by gene names fairly easily


Again, there are so many tools out there. Depending on your needs, choose the right one for you.

count how many mapped reads in a bam file

I wanted to count how many mapped reads are in a ChIP-seq bam file downloaded from UCSC ECODE project. http://genome.ucsc.edu/ENCODE/downloads.html

A quick google search:
http://left.subtree.org/
and http://seqanswers.com/forums/showthread.php?t=891

The output from short read aligners like Bowtie and BWA is commonly stored in SAM/BAM format. When presented with one of these files a common first task is to calculate the total number of alignments (reads) captured in the file. In this post I show some examples for finding the total number of reads using samtools and directly from Java code. For the examples below, I use the HG00173.chrom11 BAM file from the 1000 genomes project which can be downloaded here.

First, we look at using the samtools command directly. One way to get the total number of alignments is to simply dump the entire SAM file and tell samtools to count instead of print (-c option):

1 2	$ samtools view -c HG00173.chrom11.ILLUMINA.bwa.FIN.low_coverage.20111114.bam 5218322

If we’re only interested in counting the total number of mapped reads we can add the -F 4 flag. Alternativley, we can count only the unmapped reads with -f 4:

1 2 3 4 5 6 7

# Mapped reads only $ samtools view -c -F 4 HG00173.chrom11.ILLUMINA.bwa.FIN.low_coverage.20111114.bam 5068340 # Unmapped reads only $ samtools view -c -f 4 HG00173.chrom11.ILLUMINA.bwa.FIN.low_coverage.20111114.bam 149982

To understand how this works we first need to inspect the SAM format. The SAM format includes a bitwise FLAG field described here. The -f/-F options to the samtools command allow us to query based on the presense/absence of bits in the FLAG field. So -f 4 only output alignments that are unmapped (flag 0×0004 is set) and -F 4 only output alignments that are not unmapped (i.e. flag 0×0004 is not set), hence these would only include mapped alignments.

An example for paired end reads you could do the following. To count the number of reads having both itself and it’s mate mapped:

1 2	$ samtools view -c -f 1 -F 12 HG00173.chrom11.ILLUMINA.bwa.FIN.low_coverage.20111114.bam 4906035

The -f 1 switch only includes reads that are paired in sequencing and -F 12 only includes reads that are notunmapped (flag 0×0004 is not set) and where the mate is not unmapped (flag 0×0008 is not set). Here we add 0x0004 + 0x0008 = 12 and use the -F (bits not set), meaning you want to include all reads where neither flag 0×0004 or 0×0008 is set. For help understanding the values for the SAM FLAG field there’s a handy web tool here.

There’s also a nice command included in samtools called flagstat which computes various summary statistics. However, I wasn’t able to find much documentation describing the output and it’s not mentioned anywhere in the man page. This post examines the C code for the flagstat command which provides some insight into the output.

1 2 3 4 5 6 7 8 9 10 11 12

$ samtools flagstat HG00173.chrom11.ILLUMINA.bwa.FIN.low_coverage.20111114.bam 5218322 + 0 in total (QC-passed reads + QC-failed reads) 273531 + 0 duplicates 5068340 + 0 mapped (97.13%:-nan%) 5205999 + 0 paired in sequencing 2603248 + 0 read1 2602751 + 0 read2 4881994 + 0 properly paired (93.78%:-nan%) 4906035 + 0 with itself and mate mapped 149982 + 0 singletons (2.88%:-nan%) 19869 + 0 with mate mapped to a different chr 15271 + 0 with mate mapped to a different chr (mapQ>=5)

So, samtools flagstat  is what I need.


tommy@tommy-ThinkPad-T420:~/Desktop$ time samtools flagstat Mcf7Ctcf.sorted.bam 
20769652 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
20769652 + 0 mapped (100.00%:-nan%)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (-nan%:-nan%)
0 + 0 with itself and mate mapped
0 + 0 singletons (-nan%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

real 0m12.472s
user 0m11.237s
sys 0m0.384s

it is pretty fast.  with 12 seconds it finished counting a bam file with 20 million reads.

I used to use HTSeq python package http://www-huber.embl.de/users/anders/HTSeq/doc/overview.html
 to count the mapped reads, but it is much slower (several mins to count the same file) than samtools.


def aligned_counts(ifile1):
    '''count how many alignments are aligned back to genome, ifile1 is a sorted bam file'''
    import HTSeq
    sortedbamfile= HTSeq.BAM_Reader(ifile1)
    aligned_counts=0
    unaligned_counts=0
    for almnt in sortedbamfile:
        if almnt.aligned:
            aligned_counts+= 1
        else:
            unaligned_counts+=1
    print "number of aligned tags of %s is %d " % (ifile1, aligned_counts)
    print "number of unaligned tags of %s is %d "% (ifile1, unaligned_counts)
    return aligned_counts

Again, I am happy with the existing tools. Do not re-invent the wheels !