7. Expression analysis (with RSEM) ================================== .. shell start .. :: echo 7-expression-analysis START `date` >> /root/times.out .. :: set -x set -e rm -fr /root/rsem* rm -fr /root/bowtie* rm -fr /mnt/rsem In addition to screed, khmer, and eel-pond, you'll also need to install bowtie (see :doc:`3-big-assembly`). .. note:: You can grab the partitioned and renamed data for nematostella here:: cd /mnt curl -O http://athyra.idyll.org/~t/trinity-nematostella.renamed.fa.gz gunzip -c trinity-nematostella.renamed.fa.gz > nematostella.fa Installing rsem --------------- .. :: echo 7-expression-analysis install `date` >> /root/times.out We'll be using the `RSEM package `__ to do some expression analysis, and `EBSeq `__ to do differential expression. To install these packages, do: :: cd /root curl -O http://deweylab.biostat.wisc.edu/rsem/src/rsem-1.2.8.tar.gz tar xzf rsem-1.2.8.tar.gz cd rsem-1.2.8 make cd EBSeq make And now add this directory into your PATH, which is where Unix looks for things to run: :: echo 'export PATH=$PATH:/root/rsem-1.2.8' >> ~/.bashrc export PATH=$PATH:/root/rsem-1.2.8 Installing bowtie ----------------- If you didn't install bowtie on this machine already (e.g. as part of :doc:`3-big-assembly`), RSEM needs it; do:: cd /root curl -O -L http://sourceforge.net/projects/bowtie-bio/files/bowtie/0.12.7/bowtie-0.12.7-linux-x86_64.zip unzip bowtie-0.12.7-linux-x86_64.zip cd bowtie-0.12.7 cp bowtie bowtie-build bowtie-inspect /usr/local/bin Prepare the reference --------------------- .. :: echo 7-expression-analysis makeref `date` >> /root/times.out Go to a working directory on /mnt: :: cd /mnt mkdir rsem cd rsem Link in the nematostella file: :: ln -fs ../nematostella.fa . Make a transcript-to-gene-map file: :: python /usr/local/share/eel-pond/make-transcript-to-gene-map-file.py nematostella.fa nematostella.fa.tr_to_genes and ask RSEM to prepare the reference against which to map the reads: :: rsem-prepare-reference --transcript-to-gene-map nematostella.fa.tr_to_genes nematostella.fa nema (Here, the 'nema' at the end is what to call the reference; the other two are just file names.) This last step will take about half an hour or more. Find and list the reads ----------------------- Find the QC reads, and link them in; e.g. if using the Nematostella reads, make a volume from snap-126cc847, mount it as /data, and do: :: ln -fs /data/*.pe.qc.fq.gz . Now, make a list of the data files: :: ls -1 *.pe.qc.fq.gz > list.txt Note, the order of the files in this list is going to determine the order in the final RSEM output matrix. You might consider rearranging it so that your controls are first, etc. Run RSEM -------- .. :: echo 7-expression-analysis rsem `date` >> /root/times.out Now, for each one of the files in 'list.txt', run RSEM. This will take a long time for lots of data, so definitely run this step in screen! : :: n=1 for filename in $(cat list.txt) do echo mapping $filename gunzip -c $filename > ${n}.fq /usr/local/share/khmer/scripts/split-paired-reads.py ${n}.fq rsem-calculate-expression --paired-end ${n}.fq.1 ${n}.fq.2 nema -p 4 ${n}.fq rm ${n}.fq ${n}.fq.[12] ${n}.fq.transcript.bam ${n}.fq.transcript.sorted.bam n=$(($n + 1)) done Gather results: :: rsem-generate-data-matrix [0-9].fq.genes.results 10.fq.genes.results > 0-vs-6-hour.matrix ...and voila, you now have a file, '0-vs-6-hour.matrix', which is a tab-separated file (that Excel can load) containing a matrix of gene expression levels in FPKM (rows) vs condition (columns). The '1' condition will be the first file in list.txt, the '2' condition will be the second file, etc. If you want the conditions in a specific order, you can specify the files in the order you want -- e.g. : :: rsem-generate-data-matrix 1.fq.genes.results 3.fq.genes.results > results.matrix .. note:: Our current protocol only supports pairwise differential expression analysis, i.e. comparing two conditions, which is why we only create the 0-vs-6 hour matrix, above. .. :: echo 7-expression-analysis DONE `date` >> /root/times.out .. shell stop Next: :doc:`8-differential-expression`