In this tutorial, we will analyse differential gene/isoform expression using CuffDiff and MISO. You can run the pipeline from scratch, which will start from mapping reads with TopHat2 or you can skip to DGE analysis using our pre-mapped reads. Please keep in mind that mapping could take more than four hours. Two samples from D. melanogaster are used in this tutorial for demonstration purposes.
All programs are required to run the analysis from scratch. If you only want to analyse DGE with pre-mapped reads, only Cufflinks and MISO are needed.
Installation:
# install bowtie2
cd ~
wget -O bowtie2-2.1.0-linux-x86_64.zip http://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.1.0/bowtie2-2.1.0-linux-x86_64.zip/download
unzip bowtie2-2.1.0-linux-x86_64.zip
cp bowtie2 bowtie2-align bowtie2-build bowtie2-inspect /usr/local/bin/
# install tophat2
cd ~
wget http://tophat.cbcb.umd.edu/downloads/tophat-2.0.8b.Linux_x86_64.tar.gz
tar xvfz tophat-2.0.8b.Linux_x86_64.tar.gz
ln -s ~/tophat-2.0.8b.Linux_x86_64/tophat2 /usr/local/bin
# install cufflinks
cd ~
wget http://cufflinks.cbcb.umd.edu/downloads/cufflinks-2.1.1.Linux_x86_64.tar.gz
tar xvfz cufflinks-2.1.1.Linux_x86_64.tar.gz
cp cufflinks cuffdiff cuffcompare /usr/local/bin
# install MISO with easy_install
cd ~
wget http://pypi.python.org/packages/source/m/misopy/misopy-0.4.9.tar.gz
tar xvfz misopy-0.4.9.tar.gz
cd misopy-0.4.9/
python setup.py install
Map paired-end reads to Drosophila genome with TopHat2:
cd /mnt
# build bowtie index
bowtie2-build /data/dmel-all-chromosome-r5.51.fasta dmel_genome_bowtie2
# map reads to the genome
# this step could take more than four hours
tophat -r 300 -p 2 -o /mnt/tophat/vg1_005 dmel_genome_bowtie2 \
/mnt/OREf_SAMm_vg1_CTTGTA_L005_R1_001.fastq /mnt/OREf_SAMm_vg1_CTTGTA_L005_R2_001.fastq
tophat -r 300 -p 2 -o /mnt/tophat/w_006 dmel_genome_bowtie2 \
/data/OREf_SAMm_w_GTCCGC_L006_R1_001.fastq /data/OREf_SAMm_w_GTCCGC_L006_R2_001.fastq
Analyse differential expression with CuffDiff:
# prepare the annotation with cuffcompare
# cuffcmp.combined.gtf will be created
cd /mnt
cuffcompare -s /data/dmel-all-chromosome-r5.51.fasta -CG -r \
/data/Drosophila_melanogaster.BDGP5.71.gtf /data/Drosophila_melanogaster.BDGP5.71.gtf
# run cuffdiff, this could take more than an hour
cd /mnt
cuffdiff -p 2 -o cufflinks/comparison cuffcmp.combined.gtf \
/mnt/tophat/vg1_005/accepted_hits.bam tophat/w_006/accepted_hits.bam
First we need to sort and index BAM files:
cd /mnt
# sort two BAM files from two samples
# this could take a while
samtools sort /mnt/tophat/vg1_005/accepted_hits.bam accepted_hits.sorted
samtools sort /mnt/tophat/w_006/accepted_hits.bam accepted_hits.sorted
# index BAM files
samtools index /mnt/tophat/vg1_005/accepted_hits.sorted.bam
samtools index /mnt/tophat/w_006/accepted_hits.sorted.bam
# move BAM files to a new directory
mkdir -p /mnt/miso-data/bam-data
mv /mnt/tophat/vg1_005/accepted_hits.sorted.bam /mnt/miso-data/bam-data/vg1_005.sorted.bam
mv /mnt/tophat/vg1_005/accepted_hits.sorted.bam.bai /mnt/miso-data/bam-data/vg1_005.sorted.bam.bai
mv /mnt/tophat/w_006/accepted_hits.sorted.bam /mnt/miso-data/bam-data/w_006.sorted.bam
mv /mnt/tophat/w_006/accepted_hits.sorted.bam.bai /mnt/miso-data/bam-data/w_006.sorted.bam.bai
Then, we need to build alternative splicing index from alternative splicing annotation. Annotations for fly, mouse and human can be downloaded from MISO website.
Note that regular gene annotations are not supported by MISO:
# build alternative event index
# modENCODE_SE_4.gff is a small subset of skipped exons
# You will need to build the index for all annotation files if you
# want to do a complete analysis.
index_gff.py --index /data/modENCODE/modENCODE_SE_4.gff /mnt/miso-data/altevents/4_SE
For demonstration purposes, we will use the default settings provided by MISO. MISO provides a guide to estimate the insert size and standard deviation on its webpage. However, in this tutorial we will use the insert size and the standard deviation estimated by TopHat:
# calculate psi value (percent spliced-in),
run_events_analysis.py --compute-genes-psi /mnt/miso-data/altevents/4_SE \
/mnt/miso-data/bam-data/vg1_005/accepted_hits.sorted.bam --output-dir \
/mnt/miso-data/4_SE/vg1_005 --read-len 100 --paired-end 250 99 --event-type=SE
run_events_analysis.py --compute-genes-psi /mnt/miso-data/altevents/4_SE \
/mnt/miso-data/bam-data/w_006/accepted_hits.sorted.bam --output-dir \
/mnt/miso-data/4_SE/w_006 --read-len 100 --paired-end 250 99 --event-type=SE
# summarize the results
run_miso.py --summarize-samples /mnt/miso-data/4_SE/vg1_005/ /mnt/miso-data/4_SE/w_006/
# compare samples
run_miso.py --compare-samples /mnt/miso-data/4_SE/w_006/ \
/mnt/miso-data/4_SE/vg1_005/ /mnt/miso-data/4_SE/comparisons
MISO provides a script that filters out results based on some parameters. For example, in this tutorial, we filter out all events that have bayes-factor less than 1.
Note, bayes-factor > 10 is considered significant. However, none of events in our small dataset has bayes-factor greater than 1.0:
# filter results
filter_events.py --filter /mnt/miso-data/4_SE/comparisons/w_006_vs_vg1_005/bayes-factors/w_006_vs_vg1_005.miso_bf \
--num-inc 1 --num-exc 1 --num-sum-inc-exc 10 --delta-psi 0.20 --bayes-factor 1 \
--output-dir /mnt/miso-data/4_SE/comparisons/filtered/
Sashimi plot requires some settings written in a file, e.g. sashimi_plot_settings.txt. Among other things, the number of mapped reads is required for each dataset in order to calculate expression in RPKM unit.
To get the number of reads, run:
# You do not need to run this.
samtools view -c -f 1 -F 12 /mnt/miso-data/bam-data/w_006.sorted.bam
# total number of mapped reads = 59,324,738
samtools view -c -f 1 -F 12 /mnt/miso-data/bam-data/vg1_005.sorted.bam
# total number of mapped reads = 106,437,936
cd /mnt/miso-data/
# plot the graph from two alternative events
# both graphs will be in sashimi_plot directory
plot.py --plot-event "chr4:848485:848533:-@chr4:844127:844244:-@chr4:823024:823953:-" \
altevents/4_SE/ sashimi_plot_settings.txt --output-dir sashimi_plot
plot.py --plot-event "chr4:497918:498709:-@chr4:495586:495730:-@chr4:495079:495394:-" \
altevents/4_SE/ sashimi_plot_settings.txt --output-dir sashimi_plot