Skip to end of metadata
Go to start of metadata

You are viewing an old version of this page. View the current version.

Compare with Current View Page History

« Previous Version 15 Next »

Objectives

In this exercise, you will use a recently developed R package called RIPSeeker to identify binding regions in a small CLIP-seq dataset.

Primary Data

Ameyar-zazoua M, Rachez C, Souidi M, et al. Argonaute proteins couple chromatin silencing to alternative splicing. Nat Struct Mol Biol. 2012;19(10):998-1004.

RIPSeeker

RIPSeeker is a Bioconductor R package, and we are going to run the job using the "Rscript" utility using the script and commands file that exist in the directory you just staged.  Check out the "ripseeker_script.R" file.

source("http://bioconductor.org/biocLite.R")
biocLite("RIPSeeker")
library(RIPSeeker)
bamFiles = list.files(".", "\\.bam$", recursive=TRUE, full.names=TRUE)
ripGal=combineAlignGals(bamFiles[2], genomeBuild="hg19")
ctlGal=combineAlignGals(bamFiles[1], genomeBuild="hg19")

cNAME="input"
binSize=200
genomeBuild = "hg19"
outDir <- file.path(getwd())
reverseComplement = FALSE
uniqueHit = FALSE
assignMultihits = TRUE
rerunWithDisambiguatedMultihits = TRUE

seekOut.AGO2 <- ripSeek(bamPath = bamFiles, cNAME = cNAME, reverseComplement = FALSE, genomeBuild = "hg19", \
uniqueHit = TRUE, assignMultihits = TRUE, rerunWithDisambiguatedMultihits = FALSE, binSize=binSize, outDir=outDir)

Now, before running anything, let's go through this section by section.  First, we need to get the RIPSeeker packages from Bioconductor and load them into R.

> source("http://bioconductor.org/biocLite.R")
> biocLite("RIPSeeker")
> library(RIPSeeker)

Then, we need to read in the BAM files we're going to use with their paths and use the RIPSeeker (or rather GenomicRanges) function "combineAlignGals" to generate a "gapped alignment" object in R using hg19 genomic and contig coordinates.

> bamFiles = list.files(".", "\\.bam$", recursive=TRUE, full.names=TRUE)
> ripGal=combineAlignGals(bamFiles[2], genomeBuild="hg19")
> ctlGal=combineAlignGals(bamFiles[1], genomeBuild="hg19")

Now, we're ready to execute the primary RIPSeeker function, appropriately called ripSeek().  Since it takes so many options, we're going to declare a bunch of variables before executing the command, and then have the command refer back to those already-declared variables.  This is what it would look like:

> binSize=NULL
> cNAME="input"  #Which files are input
> genomeBuild = "hg19"   #Genome to use
> outDir <- file.path(getwd())  #Output file prefix (here, the current directory)
> reverseComplement = FALSE   #Do not reverse complement BAM entries
> uniqueHit = FALSE   #If TRUE, use only unique hits to train initial model
> assignMultihits = TRUE  #If TRUE and uniqueHit is TRUE, assign multi-hits to regions based on initial model
> rerunWithDisambiguatedMultihits = TRUE   #If TRUE and prior two options are TRUE, re-train model after multi-hit assignment
 
> seekOut.AGO2 <- ripSeek(bamPath = bamFiles, cNAME = cNAME, reverseComplement = FALSE, genomeBuild = "hg19", \
uniqueHit = TRUE, assignMultihits = TRUE, rerunWithDisambiguatedMultihits = FALSE, binSize=binSize, outDir=outDir)

The way we're going to run this at TACC is by having the actual R commands in a script called ripseeker_script.R, which is then called by a commands file using Rscript, which is a utility to run R scripts from the command line (or, in this case, a batch job) when it can't be interactive.

$BI/bin/launcher_creator.py -n ripseeker -t 01:00:00 -j ripseeker.cmds -q development -a CCBB
qsub ripseeker.sge

The output should deposit two files in your current directory:

RIPregions.gff3
seekOutput.RData

The important file is RIPregions.gff3, which is a coordinate file of the regions predicted to be enriched in the IP.  Because of the extremely slimmed-down dataset we used, there is only one peak, but this pipeline (with more reads that are actually real) would produce a file that is more interesting.

Alternatives

We could have added these lines to the ripseeker_script.R file:

biocLite("ChIPpeakAnno")
library(ChIPpeakAnno)
 
> biomart="ensembl"   #Source of genome annotations
> biomaRt_dataset="hsapiens_gene_ensembl"   #Specific dataset from biomaRt source
> goAnno="org.Hs.eg.db"   #Annotation database

Which would automatically annotate the resulting peak coordinates with their gene names.  We didn't do this because the package that is used, ChIPpeakAnno, has quite a few dependencies and takes a while to download (maybe ~15 minutes) and run.

  • No labels