Preprocesses and Aligns Run-On Sequencing (PRO/GRO/ChRO-seq) data from Single-Read or Paired-End Illumina Sequencing
Switch branches/tags
Nothing to show
Clone or download
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Permalink
Failed to load latest commit information.
images first Aug 21, 2018
input_file_exmaples add example files Aug 21, 2018
output_file_exmaples add example files Aug 21, 2018
README.md change to -4DREG Sep 17, 2018
proseq2.0.bsh change to -4DREG Sep 17, 2018

README.md

proseq2.0

Preprocesses and Aligns Run-On Sequencing (PRO/GRO/ChRO-seq) data from Single-Read or Paired-End Illumina Sequencing. It is a new version of ProseqMapper (https://github.com/Danko-Lab/utils/tree/master/proseq).

Overview

Our proseq2.0 pipeline will take single-end or paired-end sequencing reads in fastq.gz format as input. The pipeline will automate three routine pre-processing and alignment options, including

  • pre-processing reads: remove the adapter sequence and quality trim the reads (cutadapt), deduplicate the reads if UMI barcodes are used (prinseq-lite.pl)
  • mapping reads to a reference genome (BWA)
  • converting BAM files into bedGraph and BigWig formats (kentsource). When converting to bedGraph and BigWig, the pipeline only report the 5’ end position of the reads after UMI/adapter removal. For pair-end sequencing, user can choose to report the 5’ end of R1 or R2 reads.

To run our pipeline users must first download the pipeline files and install dependencies indicated in this README.md. In addition, user need to provide a path to a BWA index file and the path to the chromInfo file for the genome of choice. After running this pipeline, users should have processed data files in the specified output directory.

Dependencies:

The pipelines depend on several common bioinformatics tools:

Please make sure you can call the bioinformatics tools from your current working directory.

Usage

Preprocesses and aligns PRO-seq data.

Takes PREFIX.fastq.gz (SE),  PREFIX_R1.fastq.gz, PREFIX_R2.fastq.gz (PE)
or *.fastq.gz in the current working directory as input and writes
BAM and bigWig files as output to the user-assigned output-dir.

Requirements in current working directory:
cutadapt 1.8.3, prinseq-lite.pl 0.20.2, bwa, samtools, bedtools, and bedGraphToBigWig.

bash proseq2.0.bsh [options]

options:

To get help:
-h, --help             Show this brief help menu.

Required options:
-SE, --SEQ=SE          Single-end sequencing.
-PE, --SEQ=PE          Paired-end sequencing.
-i, --bwa-index=PATH   Path to the BWA index of the target genome
                       (i.e., bwa index).
-c, --chrom-info=PATH  Location of the chromInfo table.

I/O options:
-I, --fastq=PREFIX     Prefix for input files.
                       Paired-end files require identical prefix
                       and end with _R1.fastq.gz and _R2.fastq.gz
                       eg: PREFIX_R1.fastq.gz, PREFIX_R2.fastq.gz.
-T, --tmp=PATH         Path to a temporary storage directory.
-O, --output-dir=DIR   Specify a directory to store output in.

Required options for SE
-G, --SE_READ=RNA_5prime Single-end sequencing from 5' end of
                         nascent RNA, like GRO-seq.
-P, --SE_READ=RNA_3prime Single-end sequencing from 3' end of
                         nascent RNA, like PRO-seq.

Options for PE
--RNA5=R1_5prime    Specify the location of the 5' end of RNA
                    [default: R1_5prime].
--RNA3=R2_5prime    Specify the location of the 3' end of RNA
                    [default: R2_5prime].
                    Available options: R1_5prime: the 5' end of R1 reads
                                       R2_5prime: the 5' end of R2 reads
-5, --map5=TRUE     Report the 5' end of RNA [default on, --map5=TRUE].
-3, --map5=FALSE    Report the 3' end of RNA,
                    only available for PE [default off, --map5=TRUE].
-s, --opposite-strand=TRUE
                    Enable this option if the RNA are at the different strand
                    as the reads set at RNA5 [default: disable].

Optional operations:
--ADAPT_SE=TGGAATTCTCGGGTGCCAAGG
                    3' adapter to be removed from the 3' end of SE reads.
                   [default:TGGAATTCTCGGGTGCCAAGG]
--ADAPT1=GATCGTCGGACTGTAGAACTCTGAACG
                    3' adapter to be removed from the 3' end of R2.
                   [default:GATCGTCGGACTGTAGAACTCTGAACG]
--ADAPT2=AGATCGGAAGAGCACACGTCTGAACTC
                    3' adapter to be removed from the 3' end of R1.
                   [default:AGATCGGAAGAGCACACGTCTGAACTC]

--UMI1=0           The length of UMI barcode on the 5' of R1 read.
                   [default: 0]
--UMI2=0           The length of UMI barcode on the 5' of R2 read.
                   [default: 0]
When UMI1 or UMI2 are set > 0, the pipeline will perform PCR deduplicate.

--Force_deduplicate=FALSE
                   When --Force_deduplicate=TRUE, it will force the pipeline to
                   perform PCR deduplicate even there is no UMI barcode
                   (i.e. UMI1=0 and UMI2=0). [default: FALSE]
--ADD_B1=0         The length of additional barcode that will be trimmed
                   on the 5' of R1 read. [default: 0]
--ADD_B2=0         The length of additional barcode that will be trimmed
                   on the 5' of R2 read. [default: 0]
--thread=1         Number of threads can be used [default: 1]

-4DREG             Using the pre-defined parameters to get the most reads
                   for dREG package. Please use this flag to make the bigWig
                   files compatible with dREG algorithm.

Examples

The pipeline requires two parameters for genome information, including BWA index (--bwa-index) and chrom info (--chrom-info).

BWA index should be generated using the bwa index command according to BWA manual at http://bio-bwa.sourceforge.net/bwa.shtml . Please note that the program only take in the prefix when you assign the index, no ".bwt" in the end. See the BWA manual for more details.

Chrom info is a tab-delimited file with two columns. The first column is the chromosome name and the second is the size of the chromosome. Please see see /input_file_exmaples/mm10.chromInfo for example

bwa index Genome.fa.gz                # Don't do this if you already have bwa index
export bwaIndex=Genome.fa.gz          # Or path to your bwa index
export chromInfo=PathToChromInfo

Example 1

PREFIX.fastq.gz were made according to GRO-seq protocol as in https://www.ncbi.nlm.nih.gov/pubmed/19056941 Give UMI1=6, the pipeline will remove PCR duplicates and trim the 6bp UMI barcode.

bash proseq2.0.bsh -i $bwaIndex -c $chromInfo -SE -G -T myOutput1 -O myOutput1 --UMI1=6 -I PREFIX

Example 2

PREFIX.fastq.gz were made according to PRO-seq protocol as in https://www.ncbi.nlm.nih.gov/pubmed/23430654 UMI1 was set to 0 by default. The pipeline will NOT remove PCR duplicates.

bash proseq2.0.bsh -i $bwaIndex -c $chromInfo -SE -P -T myOutput2 -O myOutput2 -I PREFIX

Example 3

PREFIX_R1.fastq.gz and PREFIX_R2.fastq.gz were Paired-End sequenced as in chromatin run-on and sequencing (ChRO-seq) in https://www.biorxiv.org/content/early/2017/09/07/185991

  • Please note that Paired-end files require identical PREFIX and end with _R1.fastq.gz and _R2.fastq.gz.

    Assign the file use -I PREFIX. No _R1.fastq.gz, _R2.fastq.gz, nor *fastq.gz is in the end.

  • UMI1 and UMI2 were set to 0 by default. The pipeline will NOT remove PCR duplicates.

bash proseq2.0.bsh -i $bwaIndex -c $chromInfo -PE -T myOutput3 -O myOutput3 -I PREFIX

Example 4

Same as in Example 3 but with UMI barcode 4bp on the 5’ end and 2bp on the 3’ end of inserts. Pipeline will perform PCR deduplicat.

bash proseq2.0.bsh -i $bwaIndex -c $chromInfo -PE -T myOutput4 -O myOutput4 -I PREFIX --UMI1=4 --UMI2=2

Useful references

Notes for CBSUdanko users:

  1. Setup your environment to use the bioinformatics tools (e.g. prinseq-lite.pl,bedGraphToBigWig,samtools...)
export PATH=$PATH:/programs/prinseq-lite-0.20.2:/home/zw355/lib/ucsc:/programs
  1. Find the BWA index and chromosome table in the server:
export human_genome=/local/storage/data/short_read_index/hg19/bwa.rRNA-0.7.5a-r405/hg19.rRNA
export human_chinfo=/local/storage/data/hg19/hg19.chromInfo

export mouse_genome=/local/storage/data/short_read_index/mm10/bwa.rRNA-0.7.8-r455/mm10.rRNA.fa.gz
export mouse_chinfo=/local/storage/data/mm10/mm10.chromInfo

export dog_genome=/local/storage/data/short_read_index/canFam3/bwa.rRNA-0.7.8-r455/canFam3.rRNA.fa
export dog_chinfo=/local/storage/data/canFam3/canFam3.chromInfo
  1. Using --UMI1=6 to replace -b6 if you have used it in the old version (proseqMapper.bsh).

Notes for dREG users:

In order to make the most compatible with dREG algorithm, please use -4DREG flag when you process the PRO-seq and GRO-seq reads. The dREG package needs enriched reads to detect the transcriptional peaks, we use the "bwa aln" to do mappping and set lower filtering score (0) to get the most reads in this pipeline.

Here is an examples to generate the bigWig for dREG.

bash proseq2.0.bsh -SE -G -4DREG -i $dog_genome -c $dog_chinfo -I ./example1_R1.fastq.gz  -T ./tmpdir -O ./outputdir