Most RNA-seq is probably done to study gene expression. This objective implies a set of experimental methods that specifically isolate mRNA molecules within a particular expected size range. The relevant steps can include processes like rRNA depletion, size selection, and/or fragmentation. Most of this course, accordingly, focuses on pipelines and tools that assume as input mRNA-targeted data (with all that implies), and will output lists of differentially expressed genes mapped to various conditions. However, RNA has numerous biological functions beyond acting as mRNA, and many of these functions can be studied using many of the same principles. Such studies simply require taking into account the expected differences between the data of interest and gene expression-oriented data. In this section, we will discuss two major categories of 'non-mRNA' RNA-seq data: (1) sequencing to study small RNAs such as microRNAs and (2) sequencing to study RNA-protein interactions.
Many types of small RNA have been characterized, and their biological functions are extremely wide-ranging. The table below describes the different forms and biological functions of small RNAs.
Clearly, there are many biologically important functions executed by small RNA, and they can be studied by sequencing by simply cutting (for example) the 25-50bp range out of a size selection gel followed by otherwise normal library preparation. Otherwise, all these species share certain qualities that allow sequencing data derived from each to be analyzed in a similar fashion. These qualities can include (but are not limited to):
All of these issues can be taken into account effectively, and in some regards can produce results simpler to understand and evaluate relative to standard gene expression data. Our first exercise will focus on one class of small RNAs, microRNAs, and will use principles that generalize to other interesting small RNAs.
Similarly, RNA-protein interactions are required for an equally diverse set of biological functions, and hundreds of RNA-binding proteins have been identified. It is frequently interesting to isolate protein-RNA complexes, remove the protein, and sequence the resulting RNA. The methods involve combine components of RNA-seq, because the underlying molecule is RNA, and chromatin immunoprecipitation (ChIP), because the most common mechanism to isolate a protein-RNA complex is with an antibody raised against a fragment of the protein of interest. Below is a sample protocol flow for a RIP-seq experiment.
For "normal" RIP-seq, one usually expects to recover full RNA molecules regardless of where on an RNA molecule the protein was bound, since all of it is 'pulled down' together. However, such protocols generally do not use any chemical or physical means to covalently attach the RNA to the protein, which allows for the possibility that the RNA and protein complexes disassociate and re-associate from each other during sample preparation (there have been published papers that claim this - see here). Moreover, proteins will often bind to specific RNA sequence motifs or positions, and retrieval of the full RNA molecule provides no information about the specific binding site. To accommodate these concerns, methods have been developed to cross-link protein to RNA in a way that leaves a signature of interaction where the protein and RNA actually come into contact. Below is a table of the three methods that modify the RNA in various ways to enable binding site detection by sequencing.
In our second exercise, we will use a recently developed tool to analyze some sample PARCLIP data to identify specific binding sites of a protein across the entire human transcriptome.
For standard RIP-seq, many of the methods already covered in this class are useful since one can expect to recover a full RNA molecule, and IP and Mock/No Antibody/IgG samples can be thought of as "conditions" to be compared by differential expression analysis. If the only comparison is between IP and Input, than the tools you have already learned about can be used to quantitate expression for each transcript, and fold changes can be subsequently calculated. However, more specific tools do exist, particularly for CLIP-seq and its variants. Below is a table from a semi-recent paper that summarizes some of the most widely used tools in RIP and CLIP experiments.
Some of these tools, like Cuffdiff (and similar tools like edgeR and DESeq) can be used as you would for normal differential expression analysis in standard RIP-seq experiments. These programs tend to be available directly at TACC. Others, like MACS, are available at TACC but are not really designed for use with RNA-seq data. Finally, programs like RIPseeker and PARalyzer are much less widely-used (since they are much more recent), but are designed for extremely specific experimental structures. PARalyzer, for example, is explicitly and only to be used with PAR-CLIP data.
In the exercises that follow, we will use samtools to generate miRNA profiles (Exercise #1), parse Cuffdiff results to evaluate mRNA enrichment in a 'normal' RIP-seq experiment (Exercise #2), and implement PARalyzer to analyze a down-sampled toy PARCLIP dataset.
In this exercise, we will analyze a sample microRNA-seq dataset derived from H1 human embryonic stem cells that was generated for the ENCODE project and made publicly available a few years ago. These are 1x36 Illumina reads derived from all cellular RNA that is less than 200bp. Our end goal will be to obtain a microRNA profile, or counts of how many reads are derived from each microRNA.
Recall that, because these RNAs are very short, they may align multiple times throughout the genome. Moreover, our goal (as is frequently the case) is to quantitate all known small RNAs of a given class, rather than discover new members. Thus, it makes sense to align our sequences against a database of miRNA sequences (or snRNA, or tRNA, or...) where identical sequences are collapsed. As we will see, this facilitates down stream analysis and is also significantly faster since the genomic search space is dramatically reduced.
To obtain a FASTA file with all human miRNA sequences, execute these commands:
mkdir -p $SCRATCH/my_rnaseq_course/day_4b/mirbase cd $SCRATCH/my_rnaseq_course/day_4b/mirbase cp /corral-repl/utexas/BioITeam/rnaseq_course_2015/day_4b/hairpin_cDNA_hsa.fa . less hairpin_cDNA_hsa.fa |
If you run commands like wc -l or less on the reference FASTA, you will see how reduced the sequence space is. Additionally, each "contig" is a feature, meaning that rather than count reads in a genome (using something like HTSeq), we could simply count reads per contig, which is substantially faster. You can do this with ANY reference FASTA data, since most aligners can accept any FASTA as a reference index. To build the index:
module load perl module load bowtie/2.2.0 bowtie2-build hairpin_cDNA_hsa.fa hairpin_cDNA_hsa.fa |
That was fast! For comparison, doing the same command with the hg19 reference FASTA would take several hours (and should never be done on a login node). If you ls, you will see:
hairpin_cDNA_hsa.fa hairpin_cDNA_hsa.fa.1.bt2 hairpin_cDNA_hsa.fa.2.bt2 hairpin_cDNA_hsa.fa.3.bt2 hairpin_cDNA_hsa.fa.4.bt2 hairpin_cDNA_hsa.fa.rev.1.bt2 hairpin_cDNA_hsa.fa.rev.2.bt2 |
These files, together, constitute the bowtie2 reference index.
We will first copy them over from the BioITeam area on Corral, stage them in a directory in your scratch area, and look at them a little bit. The commands to do that would look something like this:
cd $SCRATCH/my_rnaseq_course/day_4b cp /corral-repl/utexas/BioITeam/rnaseq_course_2015/day_4b/human_mirnaseq.fastq.gz . less human_mirnaseq.fastq.gz |
A sample miRNA FASTQ entry, using less, might look like this:
@TUPAC_0037_FC62EE7AAXX:2:1:2000:1139#0/1 TAGCAGCACGTCAGTATTGNCGTAAAAAAAAAAAAG +TUPAC_0037_FC62EE7AAXX:2:1:2000:1139#0/1 ffafffffff\U_La[[W[B^a^abfffcccccccc |
The third line has the name attached after the "+", which is an artifact of a storage method that we won't go into here. However, everything else is basically the same - read name, followed by sequence, strand, and quality scores. However, note the string of A's towards the end. This is because, as for many very short RNAs, our read extends past the actual RNA fragment. In this case, the 'adapter' sequence is obvious - it is just a poly-A string. However, what if it wasn't? Indeed, working with publicly available small RNA data, you will often not know what the adapter is (it may not be obvious), or you might not even know if for data coming from your lab (if we're being honest).
So, we would like to use an alignment strategy that can intelligently ignore the parts that won't align to a reference (the 'adapter') and align correctly the parts that align well. This is called a 'local' alignment, in contrast to a 'global' alignment, which would count the 'mismatches' in the adapter against the alignment score. Fortunately, you have already used a local-alignment-capable aligner in this class. Tophat2 runs on the Bowtie2 alignment engine, which (if used directly, i.e. not with Tophat2), can perform local alignment.
To run the alignment, we execute a command that is very similar to BWA or Tophat2, but with different syntax:
bowtie2 --local -L 16 -x mirbase/hairpin_cDNA_hsa.fa -U human_mirnaseq.fastq.gz -S human_mirnaseq.sam |
Parameters are:
|
You can create and submit a cmds file like you have already learned how to do - this alignment would only take around 30 minutes or less. However, if TACC is being very slow (which is likely), than you can get what the results will look like by copying the result to your area with the following commands
cd $SCRATCH/my_rnaseq_course/day_4b cp /corral-repl/utexas/BioITeam/rnaseq_course_2015/day_4b/human_mirnaseq.sam . |
Recall that, because each contig in our reference is a feature, we do not need a true annotation to quantify our miRNAs. Instead, we just need to count how many reads aligned to each contig, and sort them accordingly. Fortunately, samtools has a simple utility called idxstats that reports exactly this. The following commands will produce this information by converting the SAM file to BAM, sorting and indexing it, then running idxstats:
module load samtools/0.1.19 cd $SCRATCH/my_rnaseq_course/day_4b samtools view -bS human_mirnaseq.sam > human_mirnaseq.bam samtools sort human_mirnaseq.bam human_mirnaseq.sorted samtools index human_mirnaseq.sorted.bam samtools idxstats human_mirnaseq.sorted.bam |