The goal of this is to make distribution of your data easy and consistent. The steps are not hard if you are using RStudio. The New Project wizard will do some of the work for you.
Create a directory for the R project (myproject)
Create a subdirectory (myproject/data-raw)
In the subdir, do the data manipulation that needs to be done to get your data in shape (mydata.R)
Use devtools to create the myproject/data directory and create an RDA (data) file.
This will create a file called mydata.rda in the data folder.
in the R/ folder, create a data.R file and document the data set (below)
update the DESCRIPTION file
In Rstudio, run Build and Reload to see if the package builds correctly (control shift B)
You can also run Build Source to create a tar.gz file that you can install (R CMD INSTALL mydata.tar.gz)
Create Documentation
To document our data, use the R/data.R file we created. We use roxygen2. Text formatting reference sheet.
#' @title dmdExonsHg19
#'
#' @description List of DMD exons and their positions
#'
#' @details Precalculated list of DMD exons to skip
#' in order to create an inframe transcript.
#'
#' @format A data frame of 4 columns and 79 rows
#' \code{exon_stop}, \code{skip_to_render_inframe}, \code{starstop}
#' \describe{
#' \item{exon_start}{DMD exon number at start of DMD deletion}
#' \item{exon_stop}{DMD exon number at end of DMD deletion}
#' \item{skip_to_render_inframe}{DMD exon to skip to render in-frame DMD transcript}
#' \item{startstop}{comma sep list of exon start and exon end}
#' }
#'
#' @section: Creator/maintainer
#' My Name (my@email)
"skips"
To create the documentation
Run devtools::document() or Build|Document in RStudio
Preview the documentation ?skips
Repeat til satisfied
Another option is to use RStudio to create a new package.
Just a few notes on the painful installation of EIGENSOFT on my Ubuntu 18 LTS machine. Among the errors I received were the “undefined reference to LAPCKE_dlange”. It turns out that some distributions have LAPACK as one package, which is what EIGENSOFT expects, but my system had it split into different packages.
Not sure which thing actually make the compilation work, but i did
pulling the repo from github
install liblapacke-dev
install libopenblas-dev
install libatlas
running make in src/nicksrc, src/ksrc, src/eigensrc
Finally, this compiled without errors
make LDLIBS="-llapacke " CFLAGS="-I/usr/include/" LDFLAGS="-L/usr/lib/x86_64-linux-gnu/liblapacke"
Currently, the University of California is boycotting Elsevier journals due to the onerous charges. In my field, that includes several journals such as:
AJHG
Cell, Cell Reports, Cell Stem Cell
Gene
Genomics
Lancet
European Journal of Medical Genetics
Neuromuscular Disorders
Molecular Therapy
Trends in Genetics
Here’s a list from 2012 of Elsevier journals under boycott previously.
I am studying alternative journals with similar focus. A few publishers that colleagues have mentioned are:
Biomart is an Ensembl service for querying and retrieving data. Here I’ll use R/Bioconductor to get gene positions for a list of genes.
We basically need a mart object which encapsulates the species, host and database information. Here I’m using the US west coast mirror. This defaults to the current build (GRCh38) of the human genome.
library(biomart) # from bioconductor
library(DT) # for easier viewing
mart = biomaRt::useMart(biomart='ENSEMBL_MART_ENSEMBL',
dataset='hsapiens_gene_ensembl', host='uswest.ensembl.org')
In Biomart parlance, filters are search parameters (like the WHERE clause in SQL) while atttributes are the information you are searching for.
all_filters = listFilters(mart = mart)
all_attributes listAttributes(mart=mart)
# DT just makes it easier to look at a large table
datatable(all_filters)
datatable(all_attributes)
Here, I have a list of gene symbols and I want to query Biomart. This basically says to query the database using the field external_gene_name for any values that match my list of genes. Then return to me the ensembl gene id, gene symbol, chrom, start, stop positions as a data.frame.
Dockerizing your bioconductor code is one way to ensure a reproducible workflow. It also enables easier deployment to cloud environments (e.g. AWS Batch) in case you want to scale your analysis.
In my case, I want to dockerize Bioconductor so that the summarizeOverlaps() function from GenomicAlignments is available. You could start the container, make changes to it and then commit. But the recommended route is to just make a Dockerfile to build a docker container with the libraries you want to install. Here I will install GenomicAlignments and the hg19 human transcriptome annotation. Create a file called Dockerfile and add this
# Docker inheritance
FROM bioconductor/bioconductor_docker:latest
# Install required Bioconductor package
RUN R -e 'BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")'
RUN R -e 'BiocManager::install("GenomicAlignments")'
Now you can build the docker container and tag it. It will look for a file named Dockerfile by default.
docker build -t bioconductor_docker_hg19:latest .
When you type docker images, your new image should show up.
You can then run it locally with something like
docker run -p 8787:8787 -e PASSWORD=bioc \
bioconductor_docker_hg19:latest
The password bit is for Rstudio-server so that when you fire up your browser to http://localhost:8787 you can use the username rstudio (hardcoded into the bioconductor docker) and password you want (bioc). You can change the password to anything you want but the user is always rstudio.
Some useful things to do:
docker images
docker ps -a
docker run --rm -it -v $HOME/tmp:/tmp
docker stop
docker run -it --user rstudio bioconductor/bioconductor_docker:devel bash
You can create your own launch script as a docker entrypoint for even more customization and usability.
With its long single molecule reads, Pacbio IsoSeq offers the tantalizing possibility of reading entire full-length RNAs. The library prep is not for the faint of heart and a recent provider asked me for 5ug of RNA to start – not an easy task! Moreover, a thorough interrogation of a transcriptome would require at least 4 SMRT cells, which are the “chips” that the Sequel instrument currently uses with 1 million wells. The newer version is 8M wells at was released recently.
The raw data came back as unaligned BAM files and (BAM.pbi which presumably is pacbio index) and a bunch of other files. The subreads.bam is the file of interest since that is the “sequenced” data with adapters removed. For a useful glossary of Pacbio terms see here or a bit about version3 file formats. The scraps.bam file is basically whatever remains after removing adapter and barcode sequences.
I used the IsoSeq3 workflow which basically does this:
Build a circular consensus sequence (CCS) from your subreads.bam
Classify the reads. This means finding and removing 5′ and 3′ primers from each read. Limma is the tool that does this. It produces a file named unpolished.flnc.bam.
Cluster the FLNC (full length non-chimeric) reads, which basically groups the reads into clusters of identical/similar reads (your isoforms).
Polish the final reads. Not exactly sure what this does.
This is all nicely summed up here and here. Once you have your polished reads, you can then align them to a reference genome (eg minimap2), see what genes/isoforms are expressed, and visualize in IGV.
SMRTLink
This is the software provided by Pacbio. It’s an end-to-end tool that handles sample setup, run design, data management and analysis. As such its complicated and a pain to install. Base requirements were 32 CPUs and 64GB of RAM. I got it to install on my desktop which had 8 CPUs and 24Gb of RAM but it was not fun using this. The viewer only works with Google Chrome.
I couldn’t get SMRTLink to start using the services-start command, so I ended up doing this:
Then point chrome to http://servername:9090 and you should see
Import your data
Specify SMRT Link Server, and pick Sequel Data (XML).
You need to find the XML file that describes everything. In my case it was m54098_190719_203808.sts.xml. If everything goes well, your data will be imported and some QC data will be available to see.
Overview|Status|Number of Records: number of “reads”
Loading Report|Statistics| Productivity: Number of wells (ZMWs), number of empty wells (0), wells with 1 molecule (1), wells with 2 molecules (2). Obviously a high number of Productivity 1 wells means your chip was well utilized.
Raw Data Report | Polymerase Read Length: how many bases “read”
Raw Data Report | Subread Length: how many bases “read” that we actually care about (RNA without adapters, 1x pass)
Raw Data Report | Insert Length vs Read Length: In my case, all insert reads were about 1000-2000 bases despite longer and longer polymerase reads.
Analysis
You can only analyze data that has been imported into SMRTLink.
Iso-Seq Classify only: classify reads into full-length or non-full length reads
Iso-Seq with Mapping. Runs the first option, then uses GMAP to align the outputted transcripts to a reference genome. You need to provide a reference genome.
On my modest computer (8 core, 24Gb), Iso-Seq took about 12 hours per sample (1M ZMW).
Aligning with Minimap2
From the SMRTLink Iso-Seq output, I retrieved the hq_transcripts.fastq file.
Running Minimap2 is basically the same as running BWA, so just refer to the documentation for minimap2 or pacbio instructions.
From here, you can view the results in IGV, figure out which transcripts were detected, etc. At the moment, quantifying pacbio Iso-Seq is not recommended.
Velvet is a hash-based deBrujin assembler. I occasionally use it to build contigs from short read RNAseq data, although it can be used for long reads, metagenomics, etc.
To install, just download, extract and make.
There are two parts to the program:
1. velveth – prepare the sequences
2. velvetg – build the graph and output the contig.fa and some summary statistics
GMAP is an old (circa 2006) software for long read alignment. Its use case is for mapping RNA reads back to a genome. It has found new life in the world of long read RNAseq such as from Pacbio reads.
Perhaps because of its age and architecture, it has some quirks and dependencies that seem odd to the modern bioinformatician.
Extract the files using
tar zxvf gmap-2018-05-30.tar.gz
change into the source directory and edit the config.site file. config.site
sets some system wide parameters the most important will be the path to your
GMAP database, which is where the index, one you create it, will be placed.
Then compile from the command line
./configure
make
make install
It will help to read the notes README and INSTALL for particular issues with different architectures, operating systems, etc. Because of the many optimizations, some deep diving may be required.
Building an index
There might be some prebuilt indexes, but I couldn’t find them. To build your own you need to specify fasta file(s). One per chromosome or one for the whole genome or contigs or whatever.
# in this case, i called it hg19
gmap_build -d nameOfGenome /path/to/genome_fasta_files
The index files will be placed in whatever you set in config.site or on the gmap_build command line options.
Align a cdna
To align, there are a few options.
map
map and align (show alignment)
align only
batch mode
Since I had several fastqs in one file, I needed to convert FASTQ to FASTA.
Using seqtk
seqtk seq -a input.fastq > output.fasta
And align
# this aligns
gmap -d hg19 myfile.fasta
# this aligns and maps (shows alignment)
gmap -d hg19 -A myfile.fasta
Examining output
The output has genome sequence on the top row and the cDNA on the bottom with the cDNA in its forward direction.
Numbers below introns indicate its length in nucleotides.
GMAP uses these symbols:
| match
> canonical intron
- gap
= non-canonical intron
Figure 4 in the Wu et al paper is helpful in understanding the output. [Wu & Watanabe. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics May 2005]
Its handy to have a webserver where you can put up BAM files for sharing with collaborators. But if you have tons of files or available webserver, you can use Amazon S3 to serve up your files. While S3 is pretty simple, setting it up for serving BAM (and index) files has a few tricky aspects to it.
Setting up an S3 bucket
You can do this pretty simply from the AWS dashboard. Pick a region (close to you), create a folder or two. You can specify an encryption type so that data in S3 is encrypted (I used server side encryption with Amazon S3-Managed Keys).
Remember S3 bucket names are unique! It cannot be the same as any other bucket that has ever existed. Also, it’s a good idea to make the bucket private so that the BAM files are not open to the public.
Copy files to S3
I’ll assume you have AWS CLI for copying files. To list your buckets
aws s3 ls
# list a specific archive
aws s3 ls s3://mybucket
We want to copy our files into S3. We need both the BAM file and the BAM index file.
S3 does not support authentication or any such mechanism. So if your bucket is private, you will need to use some methods to expose your data. There are a few options.
AWS CLI to presign a URL
Using the awscli, you can generate a signed URL. Basically this generates a URL that will work for anyone who has the link. The URL will expire (default 3600 seconds). We actually need to make two signed URLs. One for the BAM and one for the BAM index file.
# this will make one signed URL for the BAM file
aws s3 presign s3://mybucket/path/to/mybam.bam
# this makes a signed URL for the BAM index file
aws s3 presign s3://mybucket/path/to/mybam.bam.bai
# to specify an expiry time
aws s3 presign s3://mybucket/path/to/mybam.bam --expires-in
View in IGV
Then open up IGV, use File, Load from URI and input the signed URL for both the BAM and BAM index. Note that my standalone binary distribution version of IGV did not have two places to input the BAM and Bam index signed URL (there was only one input box). Instead, I had to goto the Broad website and load from Java Web Start which started IGV after downloading from the Broad. This version of IGV (2.4.10 03/20/2018) did have two input boxes for the BAM and BAM index.
Generate presigned URL using Java Eclipse
There is a code example here. A few steps need to be completed:
download the AWS Java SDK
add the SDK jar files to your build path
You must specify these things:
clientRegion (such as us-west-2)
bucketName (mybucket)
objectKey (name of the object stored in S3). Note that if you have a subfolder, that goes into the objectKey name (examples/bamfiles/mybam.bam). This is because S3 doesn’t really have folders, so the object name actually includes the whole path.
