Content Comparison

Table of Contents

Overview and Objectives

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.

Studying Small RNAs

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.

Image Removed

Yao Y, Sun Q. Exploration of small non coding RNAs in wheat (Triticum aestivum L.). Plant Mol Biol. 2012;80(1):67-73.

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):

...

Table of Contents

Overview and Objectives

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.

Studying Small RNAs

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 some small and/or non-coding RNAs (though small RNAs are, almost by definition, non-coding).

Image Added

Yao Y, Sun Q. Exploration of small non coding RNAs in wheat (Triticum aestivum L.). Plant Mol Biol. 2012;80(1):67-73.

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):

• Single-end sequencing: if an RNA-species is (for example) 16-25 bp, than paired-end sequencing of any kind provides little (though some) additional data compared to a single end run provided the reads are long enough
• Increased adapter contamination: as implied above, adapter sequence is almost always included in your reads requiring either pre-processing to remove such sequence or alignment adjustment to account for it
• Low-complexity libraries: there are often far fewer members of a category of small RNA in a genome than there are reads in the data, meaning the exact same sequence will occur in multiple reads
• Extensive genomic duplication: there are often many copies of the same sequence of a given small RNA in a genome, meaning most genomic alignments will contain numerous "multi-mappers"

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.

Zhao J, Ohsumi TK, Kung JT, et al. Genome-wide identification of polycomb-associated RNAs by RIP-seq. Mol Cell. 2010;40(6):939-53.

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.

König J, Zarnack K, Luscombe NM, Ule J. Protein-RNA interactions: new genomic technologies and perspectives. Nat Rev Genet. 2011;13(2):77-83.

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.

Li Y, Zhao DY, Greenblatt JF, Zhang Z. RIPSeeker: a statistical package for identifying protein-associated transcripts from RIP-seq experiments. Nucleic Acids Res. 2013;41(8):e94.

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 real (but down-sampled toy PARCLIP dataset) PARCLIP data (Exercise #2).

Exercise #1: miRNA Sequencing and Profiling (miRNA-seq)

...

module load perl
module load bowtie/2.2.0
bowtie2-build hairpin_cd $SCRATCH/my_rnaseq_course/day_4b/mirbase
bowtie2-build hairpin_cDNA_hsa.fa hairpin_cDNA_hsa.fa

...

To run the alignment, we execute a command that is very similar to BWA or Tophat2, but with different syntax:

cd $SCRATCH/my_rnaseq_course/day_4b
bowtie2 --local -N 1 -L 16 -x mirbase/hairpin_cDNA_hsa.fa -U human_mirnaseq.fastq.gz -S human_mirnaseq.sam

...

Photo-activatable ribonucleoside cross-linking and immunoprecipitation followed by sequencing (PARCLIP-seq) is one of several emerging technologies for studying protein-RNA interactions with dramatically increased nucleotide resolution compared to traditional RIP-seq methods.  Earlier, we discussed the differences between RIP, CLIP, and PARCLIP.  In this exercise, we will take some public PARCLIP data and analyze it using SAMtools and PARalyzer, which is a tool developed to identify the specific mutational signature indicating real protein-RNA interaction.

Data Staging and PARalyzer Usage

PARalyzer can be downloaded from (and  We are using heavily reduced/downsampled data from a paper by Kishore et al (found here) that targeted human Argonaute2 (Ago2).

Data Staging and PARalyzer Usage

PARalyzer can be downloaded from (and its documentation found) here.  It was created by one of the first groups to every successfully publish a PARCLIP dataset.  Since it is not in the TACC module system, you can copy it and all necessary associated files to your scratch area by executing the following commands:

cd $SCRATCH/my_rnaseq_course/day_4b
cp -r /corral-repl/utexas/BioITeam/rnaseq_course_2015/day_4b/PARalyzer_v1_5 .
lscd -la PARalyzer_v1_5
ls -la

As you will see, this directory contains many files, including the PARalyzer executable itself.  It also contains the SAM file that we will analyze.  A BAM file is not provided, since PARalyzer (this is fairly unusual) requires that input alignments be in SAM format, without headers.  The other two files that are worth mentioning are "hg19.2bit", which is a special binary form of the hg19 human genome build, and "sample.ini," which contains all specifications for PARalyzer.  In fact, to run PARalyzer, we provide it only with a memory allotment and the name of the sample.ini file.  Go ahead and see what is in sample.ini using less:

less PARalyzer_v1_5/sample.ini
 
BANDWIDTH=3
CONVERSION=T>C
MINIMUM_READ_COUNT_PER_GROUP=1
MINIMUM_READ_COUNT_PER_CLUSTER=5
MINIMUM_READ_COUNT_FOR_KDE=5
MINIMUM_CLUSTER_SIZE=5
MINIMUM_CONVERSION_LOCATIONS_FOR_CLUSTER=1
MINIMUM_CONVERSION_COUNT_FOR_CLUSTER=1
MINIMUM_READ_COUNT_FOR_CLUSTER_INCLUSION=5
MINIMUM_READ_LENGTH=1
MAXIMUM_NUMBER_OF_NON_CONVERSION_MISMATCHES=0
BOWTIE_FILE=./_____.sam
GENOME_2BIT_FILE=./hg19.2bit
OUTPUT_CLUSTERS_FILE=./paralyzer_output
EXTEND_BY_READ

All of these are parameters used by PARalyzer, and we don't have time to get into each one.  If you head to the PARalyzer website, they have a README file that contains full descriptions.  The option we care about, for our purposes, is "BOWTIE_FILE" which includes the path to our input SAM file.

Running PARalyzer

Results and Visualization

  Go ahead and take a look at the SAM file using less.  In particular, below is a sample read that contains the characteristic mutational signal that (in principle) indicates close proximity to a protein during crosslinking:

SRR189785.13517867      0       chr1    565558  255     24M     *       0       0       CATACTCCTCAACTACCCACATAG        ?=7>AC>B>@CA<BB<BB@>94AA        XA:i:1  MD:Z:12T11      NM:i:1

The above read demonstrates the characteristic T>C mutation in the center of the read that (taken together with several such events), implies real cross-linking.  In contrast to pure read counting, identification of clusters with many similar mutations is the objective of our analysis.

The other two files in the PARalyzer directory that are worth mentioning are "hg19.2bit", which is a special binary form of the hg19 human genome build, and "sample.ini," which contains all specifications for PARalyzer.  In fact, to run PARalyzer, we provide it only with a memory allotment and the name of the sample.ini file.  Go ahead and see what is in sample.ini using less or cat:

less sample.ini
 
#BANDWIDTH=3
#CONVERSION=T>C
#MINIMUM_READ_COUNT_PER_GROUP=10
#MINIMUM_READ_COUNT_PER_CLUSTER=1
#MINIMUM_READ_COUNT_FOR_KDE=1
#MINIMUM_CLUSTER_SIZE=1
#MINIMUM_CONVERSION_LOCATIONS_FOR_CLUSTER=2
#MINIMUM_CONVERSION_COUNT_FOR_CLUSTER=1
#MINIMUM_READ_COUNT_FOR_CLUSTER_INCLUSION=1
#MINIMUM_READ_LENGTH=1
#MAXIMUM_NUMBER_OF_NON_CONVERSION_MISMATCHES=5
#SAM_FILE=./Kishore_Ago2_parclip.sam
#GENOME_2BIT_FILE=./hg19.2bit
#OUTPUT_CLUSTERS_FILE=./Kishore_Ago2_clusters.csv
#EXTEND_BY_READ

All of these are parameters used by PARalyzer, and we don't have time to get into each one.  If you head to the PARalyzer website, they have a README file that contains full descriptions.  The option we care about, for our purposes, is "SAM_FILE" which includes the path to our input SAM file, and "OUTPUT_CLUSTERS_FILE" which specifies the comma-separated output file.

Running PARalyzer

To bring up the PARalyzer help page (which is very minimal), execute these commands:

cd $SCRATCH/PARalyzer_v1_5
PARalyzer 4g

Since the bulk of the tool is written in Java, the '4g' specifies how much memory we allocate to the Java virtual machine.  To actually perform the analysis, we would run the command (either at the command line or through a batch submission file):

cd $SCRATCH/PARalyzer_v1_5
PARalyzer 4g sample.ini

Running this command would produce the output file Kishore_Ago2_clusters.csv, but would take quite a bit time since the tool does not scale particularly well to increasing read depths.  Consequently, I have already prepared the output file, and named it Kishore_Ago2_clusters_done.csv.  Feel free to try to generate the output at some point, if you like, but in the interests of time, we will proceed to parsing the PARalyzer output to identify interesting interactions.

Results Parsing

Use head to look at the Kishore_Ago2_clusters_done.csv file using the following commands:

cd $SCRATCH/PARalyzer_v1_5
head Kishore_Ago2_clusters_done.csv
 
#Chromosome,Strand,ClusterStart,ClusterEnd,ClusterID,ClusterSequence,ReadCount,ModeLocation,ModeScore,ConversionLocationCount,ConversionEventCount,NonConversionEventCount
#chr1,+,565558,565581,G1.1,CATACTCCTCAATTACCCACATAG,19,565571,0.6521311984865596,2,4,110
#chr1,+,565635,565662,G2.1,ACTATTTATATTATCCTAACTACTACCG,36,565662,0.6357747408016309,2,2,352
#chr1,+,565663,565686,G2.2,CATTCCTACTACTCAACTTAAACT,85,565664,0.6581599027923404,7,10,630
#chr1,+,566884,566910,G3.1,ACCTAACCATCTTCTCCTTACACCTAG,22,566910,0.6731219485026998,3,4,172
#chr1,+,567716,567744,G4.1,TAAATCTAACTTTCTTCCCACAACACTTT,88,567726,0.6173693006052031,6,8,941
#chr1,+,16266654,16266679,G8.1,GTTTTTTTTAATCTGTGCCAAAAATG,41,16266668,0.8083176199348235,4,39,425
#chr1,+,32637046,32637078,G12.1,TTATTTAATTTTTTTTTCTTTTGCACATTTTTC,14,32637065,0.728509998913673,4,11,279
#chr1,+,36645990,36646010,G15.1,TTTTTAATCTGCACCTTATAG,18,36646000,0.8427833343143215,2,19,159
#chr1,+,41220043,41220064,G18.1,TGTAAACATCCTTGACTGGAAG,24,41220050,0.7405424084870008,3,12,131

The first four fields provide the genomic coordinates of each read cluster, while the other fields each provide various pieces of information about the cluster.  The most important field is likely "ModeScore," which represents the "strength" of the signal from that cluster.  Also, the "sequence" is fairly important as well, since PARalyzer has a clear bias towards read groups with poly-T chains (which is unsurprising, given the signal we are looking for).  These can be accommodated with pre-processing or trimming for poly-A or poly-T sequences.

Now, we want to view the highest confidence clusters, and see which transcripts, if any, they interact with.  To do so, we can pipe together a set of commands that will give us what we want.  For example:

cd $SCRATCH/PARalyzer_v1_5
sort -r -t , -k 9 -n Kishore_Ago2_clusters_done.csv | head

#chr17_gl000205_random,-,86334,86355,G751.1,TGCGCAGAACCTTCCCATCCTT,34,86355,0.9996957487817257,2,6,157
#chr14,+,99289965,99289986,G537.1,AAAAGCCGGGTTGAGAGGGTGA,15,99289986,0.9960425240641467,2,16,27
#chr17,+,76220266,76220302,G703.1,AATAGCACAAACTACAATTAAAACTAAGCACAAAGCC,33,76220266,0.9954430928200928,4,32,96
#chr3,-,49397535,49397566,G1141.1,CCTTTTTCATTTATCTATAATTTACCTAAGAT,10,49397535,0.9877804735048394,2,9,149
#chrM,+,4869,4900,G1739.1,ACAAAAACTAGCCCCCATCTCAATCATATACC,12,4869,0.9781133037398211,3,3,51
#chr19,+,33668019,33668058,G812.1,CGAATCCCAGCGGTGCCTCAACCGAGCGTCCAAGCTCTTT,11,33668030,0.9776451604002823,2,3,51
#chr7,+,10493027,10493065,G1469.1,CACAGGGGAAGGAAATAACATTGCACTTTATAAACACTG,52,10493027,0.9735726980290786,5,48,356
#chr20,+,3898183,3898212,G980.1,GGTCAAGCAGCATTGTACAGGGCTATGAAA,504,3898183,0.9717404318250604,6,290,2229
#chr4,-,109818065,109818085,G1236.1,TGTGGAATAGTTTAAAACTGT,32,109818085,0.9668473986530141,2,39,196
#chr12,+,31549462,31549500,G380.1,CACAGGGGAAGGAAATAACATTGCACTTTATAAACACTG,71,31549462,0.965539145175358,6,68,478

Above, we see the top ten clusters as identified by PARalyzer.  This is not the only useful ranking metric, since one might want to rank clusters by read count, or length, or whatever - but it gives some useful starting information.  Let's look at the mutational signature of the reads that are within a top cluster in more detail using SAMtools with the following commands:

module load samtools/0.1.19
cd $SCRATCH/PARalyzer_v1_5
samtools view Kishore_Ago2_parclip.bam chr17:76220260-76220302 | cut -f 13 | less

Clearly, there are MANY cases of T>C mutations that (in reality - you can't see this from the above commands alone) are all at the same nucleotide.  If one were to view the reads in IGV or another browser, one would expect to see something like this:

Image Added

Figure from Hafner M, Lianoglou S, Tuschl T, Betel D. Genome-wide identification of miRNA targets by PAR-CLIP. Methods. 2012;58(2):94-105.

Data from Lipchina I, Elkabetz Y, Hafner M, et al. Genome-wide identification of microRNA targets in human ES cells reveals a role for miR-302 in modulating BMP response. Genes Dev. 2011;25(20):2173-86.

In the figure, red represents reference bases and yellow represents T>C mutation events.  miRNA binding sites are shown because this PAR-CLIP experiment targeted Ago2, so the regions retrieved were primarily miRNA binding sites.