How to get help

ICETEA questions could be posted to the Bioconductor support site, which serves as a searchable knowledge base of questions and answers.

Posting a question and tagging with “icetea” will automatically send an alert to the package authors to respond on the support site. See the first question in the list of Frequently Asked Questions (FAQ) for information about how to construct an informative post.

How to cite ICETEA

To cite icetea package in your analysis type citation("icetea") from within your R session.

De-multiplexing the fastq

Experiments like MAPCap and RAMPAGE (Batut et al. 2013) produce multiplexed fastq files with sample indicies and PCR barcodes attached to the fastq sequence. The tool demultiplexFASTQ de-multiplexes the samples using the sample barcode, producing fastq files corresponding to each sample. It also trims off these barcodes and attaches them in the header of the fastq files for further processing.

For RAMPAGE (Batut et al. 2013) protocol, the modified header looks like this : <read ID>#<sample barcode>:<pseudo-random barcode> For MAPCap protocol, the modified header looks like this : <read ID>#<sample barcode>:<random barcode>:<replicate barcode> For conventional CAGE (Kodzius et al. 2006) protocols, we assume that standard Illumina de-multiplexing has already been performed and this step is not required.

Function returns a modified CapSet object that contains location of demultiplexed files along with processing statistics.

# demultiplex fastq and trim the barcodes
dir.create("splitting")
cs <- demultiplexFASTQ(cs, max_mismatch = 1, outdir = "splitting", ncores = 10)
#> de-multiplexing the FASTQ file

It takes about 6 min (340 sec) to trim and de-multiplex 1M PE reads into 12 sampels (12 pair of fastqs), using 1 thread. This can be done under 45 seconds if 10 threads are used.

Optionally, we can skip the above process and perform post-mapping de-multiplexing on the BAM files (see below).

Mapping the fastqs

The demultiplexed fastqs can now be mapped using the mapCaps function. This function is a wrapper over the subjunc function from Rsubread package (Liao, Smyth, and Shi 2013). It additionally performs sorting and collects mapping statistics of the mapped files, stored in the modified CapSet object.

In order to run the function we first create a subread index of our genome.

dir.create("genome_index")
library(Rsubread)
buildindex(basename = "genome_index/dm6", reference = "/path/to/dm6/genome.fa")

We can now perform the mapping.

# provide location of a subread index file
subread_idx <- "genome_index/dm6"
# map fastq
cs <- mapCaps(cs, subread_idx, outdir = "mapping", ncores = 20, logfile = "mapping/subread_mapping.log")
# you can save the CapSet object for later use
save(cs, file = "myCapSet.Rdata")

Note: The package Rsubread is not available for windows. Windows users would need to map their demultiplexed files using another tool.

Note: Optionally, we can perform mapping directly on the multiplexed files and perform the de-multiplexing afterwards (see below). This option is only recommended where the FASTQ files are not very large (for example, where overall sequencing depth is low, below 10 Mil).

(optional) Post mapping de-multiplexing

Demultiplxeing can be performed on the BAM files after mapping, using the splitBAM_byIndex and (for MAPCap) splitBAM_byRepindex functions. These functions assume that the reads have been trimmed and the sample barcodes have been attached to the header in the way as described above.

splitBAM_byIndex(bamFile, index_list, outfile_list, max_mismatch = 1, ncores = 10)

Filtering PCR duplicates

Experiments like MAPCap and RAMPAGE provide us a way to remove sequencing reads which are PCR duplicates from the mapped data. Random barcodes present in the read sequence are used for this purpose. In MAPCap, pre-designed random barcodes present in the oligos serve as UMIs, while in RAMPAGE, the sequences used as RT-PCR primers are treated as pseudo-random barcodes. PCR duplicates are recognized as reads that map to the same 5’ position and contain the same random barcode.

The function filterDuplicates removes these PCR duplicate sequences (keeping only one copy in these cases), and creates de-duplicated BAM files. It returns the modified CapSet object with de-duplication statistics. filterDuplicates works in a strand-specific way. By default, it filters the R2 read from the paired-end BAM files, along with secondary alignments. Both reads can be allowed to kept using the option keepPairs = TRUE. Note that in this case, each read in a pair is evaluated seperately and in case of fragments where only one mate is flagged as duplicated, the other mate would be kept. In order to keep only fragments with both kept mates, the output BAM files should be filtered for proper pairs using samtools. Note that the run-time for paired-end evaluation is significantly longer, and only the first mate is considered for TSS detection (later).

# load a previously saved CapSet object (or create new one)
cs <- exampleCSobject()
#> Checking de-multiplexed R1 reads
#> Checking de-multiplexed R2 reads
#> Checking mapped file
#> Checking de-duplicated file
# filter PCR duplicates and save output in a new directory
dir.create("removedup")
cs <- filterDuplicates(cs, outdir = "removedup")
#> Removing PCR duplicates : /tmp/Rtmp3j0kSQ/Rinst1a4c550fb36e/icetea/extdata/bam/embryo1.bam
#> Removing PCR duplicates : /tmp/Rtmp3j0kSQ/Rinst1a4c550fb36e/icetea/extdata/bam/embryo2.bam
#> Removing PCR duplicates : /tmp/Rtmp3j0kSQ/Rinst1a4c550fb36e/icetea/extdata/bam/embryo3.bam
#> Removing PCR duplicates : /tmp/Rtmp3j0kSQ/Rinst1a4c550fb36e/icetea/extdata/bam/embryo4.bam

Tagging of read IDs using the random barcodes and subsequent removal of PCR duplicates can also be performed using the command line tools, like UMItools.

Detection of TSS

icetea implements a new method of detection of transcription start sites, which is adopted from the recently described methods for differential transcription factor binding analysis in ChIP-Seq data (Lun and Smyth 2015). Genome is divided into small, sliding windows and the TSSs are detected as the windows that show a user-defined fold-enrichment over a local background region. Multiple consicutively enriched windows are then merged to detect broader TSSs. The method works well with replicates, resulting in one TSS file per group.

This method is implemented in the function detectTSS, which returns a modified CapSet object with TSS detection statistics.

# detect TSS
dir.create("tssCalling")
cs <- detectTSS(cs, groups = c("wt", "wt", "mut", "mut"),
                 outfile_prefix = "tssCalling/testTSS", restrictChr = "X", ncores = 1)
#> Counting reads within detected TSS
#> Writing filtering information as .Rdata
# export the detected TSS bed files
exportTSS(cs, merged = TRUE, outfile_prefix = "tssCalling/testTSS")
#> Writing merged .bed files

Using spike-In controls

For the differential TSS expression analysis, the function fitDiffTSS shown above can utilize one of the available internal normalization methods. In some cases however, external (spike-in) normalizations are preferred. icetea provides a way to perform spike-in normalization during differential expression analysis.

Spike-In samples could be processed in the same way as the normal samples, using the CapSet object as described above. Normalization factors can be obtained from the spike-in reads using getNormFactors and provided to the fitDiffTSS function during detection of differential TSS between samples.

## get gene counts for spike-in RNA mapped to human genome
library("TxDb.Hsapiens.UCSC.hg38.knownGene")
normfacs <- getNormFactors(cs_spike, features = genes(TxDb.Hsapiens.UCSC.hg38.knownGene))

csfit <- fitDiffTSS(cs, groups = rep(c("wt","mut"), each = 2),
                    normalization = NULL, normFactors = normfacs,
                    outplots = NULL, plotRefSample = "embryo1", ncores = 1)

Getting the gene counts

Given the annotations, TSS counts for each gene can be summarized to gene counts. The following function sums up the TSS counts per gene from a given txdb object and returns gene counts.

# get transcripts
dm6trans <- transcriptsBy(TxDb.Dmelanogaster.UCSC.dm6.ensGene, "gene")
# get gene counts, counting reads around 500 bp of the TSS
gcounts <- getGeneCounts(cs, dm6trans)

Annotating the TSS distribution

Detected or differentially expressed TSS, which are exported as bed files, can be quickly annotated using the function annotateTSS. annotateTSS reports the numbers/proportions of detected TSS falling into different genomic features (provided by a TxDb object) as a data.frame and optionally as a plot. In order to break ties between overlapping features, we can provide a vector of feature names in the decreasing order of preference.

# save the output as data.frame (outdf) + plot on screen
outdf <- annotateTSS(tssBED = tssfile,
            txdb = TxDb.Dmelanogaster.UCSC.dm6.ensGene,
            plotValue = "number",
            outFile = NULL)
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned many:1 mapping between keys and columns
#> 'select()' returned 1:1 mapping between keys and columns
#> 'select()' returned many:1 mapping between keys and columns
#> 'select()' returned many:1 mapping between keys and columns
#> 'select()' returned 1:1 mapping between keys and columns