refTrans

conf.py

@@ -45,7 +45,7 @@
 
 # General information about the project.
 project = u'khmer-protocols'
-copyright = u'2013, C. Titus Brown et al.'
+copyright = u'2015, C. Titus Brown et al.'
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
@@ -54,7 +54,7 @@
 # The short X.Y version.
 version = '0.8'
 # The full version, including alpha/beta/rc tags.
-release = '0.8.4'
+release = '0.8.5'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

index.txt

@@ -51,12 +51,23 @@ This is a protocol for assembling low- and medium-diversity metagenomes.
 Marine sediment and soil data sets may not be assemblable in the cloud
 just yet.
 
+Reference based RNAseq assembly: the refTrans Protocol
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+:doc:`refTrans/index`
+
+This is a protocol for reference based transcriptome assembly. 
+Materials are based on two workshops given at MSU by Titus Brown and Matt Scholz on Dec 4/5 ad Dec 10/11. 
+Both workshops were sponsored by `iCER <http://icer.msu.edu/>`__ 
+and made use of the `MSU High Performance Compute Center <http://hpcc.msu.edu/>`__ .
+
+
 Additional information
 ----------------------
 
 Need help? Either post comments on the bottom of each page, OR
 `sign up for the mailing list <http://lists.idyll.org/listinfo/protocols>`__.
-
+ 
 Have you used these protocols in a scientific publication?  We'll have
 citation instructions up soon.
 

refTrans/index.txt

@@ -0,0 +1,31 @@
+===================================================
+Reference based RNAseq assembly (refTrans) protocol
+===================================================
+
+:author: Tamer A. Mansour and Titus Brown
+
+This is a protocol for reference based transcriptome assembly. 
+Materials are based on and extending two workshops given at MSU by Titus Brown and Matt Scholz on Dec 4/5 ad Dec 10/11. 
+Both workshops were sponsored by `iCER <http://icer.msu.edu/>`__ 
+and made use of the `MSU High Performance Compute Center <http://hpcc.msu.edu/>`__ .
+
+.. figure:: files/model-rnaseq-pipeline.png
+
+Tutorials:
+
+.. toctree::
+   :maxdepth: 2
+
+   m-dataNsoftware_HPC
+   m-quality
+   m-tophat
+   m-count
+   m-data-analysis
+   m-func-analysis
+   m-advice
+
+Reference material
+------------------
+:doc:`../docs/command-line`
+
+.. :doc:`../amazon/index`

refTrans/m-HPC-cuffdiff.txt

@@ -0,0 +1,62 @@
+Submitting cuffdiff job to the MSU HPC queue
+============================================
+
+.. shell start
+
+
+We should start by specifying the resources required for our job and load the needed softwares
+::
+
+   echo "#PBS -l walltime=04:00:00,nodes=1:ppn=8,mem=12Gb" > run_cuffdiff.sh
+   echo "module load cufflinks/2.1.1" >> run_cuffdiff.sh
+
+Define the groups you want to compare (According to our sample naming conventions, 
+it should be exactly matching the stating name of your samples)
+::
+
+   echo "groups=(lung salivary)" >> run_cuffdiff.sh
+
+We will prepare the replicates as comma separated lists to pass them to the shell scripts
+::
+
+   gp_replicates=()
+   for gp in "${groups[@]}"; do
+        replicates=()
+        for rep in $gp*_tophat; do
+            replicates+=("$workingPath/$rep/accepted_hits.bam")
+        done
+        gp_replicates+=($(echo ${replicates[*]} | tr ' ' ','))
+   done
+
+   echo "
+   cuffdiff \
+    --output-dir "$workingPath/cuffdiff_out" \
+    --num-threads 8 \
+    --multi-read-correct \
+    --frag-bias-correct ${homo_sapiens_genome} \
+    --labels ${groups[0]},${groups[1]} \
+    $GTF_file  ${gp_replicates[0]}  ${gp_replicates[1]}" >> run_cuffdiff.sh
+
+
+you can see your script::
+
+   less run_cuffdiff.sh
+
+(use 'q' to exit). 
+
+
+Now submit your job
+::
+
+   qsub run_cuffdiff.sh
+
+
+.. note:: You will need to adjust the required wall time, memory and number of nodes(processor) required for your job
+
+You can use 'qstat | grep <username>' to check on your jobs' status.
+
+.. shell stop
+
+----
+
+Back: :doc:`m-cuffdiff`

refTrans/m-HPC-tophat.txt

@@ -0,0 +1,44 @@
+Submitting TopHat jobs to the MSU HPC queue
+===========================================
+
+.. shell start
+
+In :doc:`m-more-tophat`, we showed you a *shell script*, which was a way
+of telling the computer to do multiple things in a row.  We discussed
+several advantages to scripting --
+
+1. It automates long-running processes;
+2. It tracks what you did, and you can edit it (to modify analyses) and
+   also copy it (to do a family of analyses);
+3. You can also provide it as part of your Methods in your paper;
+
+There is also a fourth advantage, or really a necessity, to scripting:
+it's how you use the MSU HPC to run your computation.  Briefly, to run
+a job on the HPC, you create a shell script and then run 'qsub'.
+
+This code will automate making the scripts for each sample files and run the qsub for these scripts
+::
+
+   mkdir qsub_scripts
+   for f in *R1.pe*; do
+        input_R1="$f"
+        input_R2=$(echo "$f" | sed s/R1/R2/)
+        input_R=$(echo "$f" | sed s/R1.pe/R.se/)
+        output=$(basename "$f" | cut -f 1,2 -d "_")_tophat
+        script="./qsub_scripts/"$output"_script.sh"
+        echo "#PBS -l walltime=00:30:00,nodes=1:ppn=4,mem=4Gb" >> $script
+        echo "module load TopHat2/2.0.12" >> $script
+        echo "cd \$PBS_O_WORKDIR" >> $script
+        echo "tophat -p 4 -G "$GTF_file" --transcriptome-index "$workingPath/trans_Index/transcriptome" -o $output $Bowtie2Index $input_R1 $input_R2,$input_R" >> $script
+        qsub $script
+   done
+
+.. note:: You will need to adjust the required wall time, memory and number of nodes(processor) required for your job
+
+You can use 'qstat | grep <username>' to check on your jobs' status.
+
+.. shell stop
+
+----
+
+Back: :doc:`m-tophat`

refTrans/m-advice.txt

@@ -0,0 +1,71 @@
+Miscellaneous advice
+====================
+
+Sequencing depth and number of samples
+--------------------------------------
+
+`Hart et al. (2013)
+<http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3842884/>`__ provides a
+nice description and a set of tools for estimating your needed
+sequencing depth and number of samples.  They provide an `Excel based
+calculator
+<http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3842884/bin/supp_data.zip>`__
+for calculating number of samples.  Their numbers are surprisingly
+large to me ;).
+
+In a proposal for an exploratory effort to discover differentially
+expressed genes, I would suggest 3-5 biological replicates with 30-50
+million reads each.  More reads is usually cheaper than more replicates,
+so 50-100m reads may give you more power to resolve smaller fold changes.
+
+Developing your own pipeline
+----------------------------
+
+Even if all you plan to do is change the filenames you're operating on,
+you'll need to develop your own analysis pipeline.  Here are some tips.
+
+1. Start with someone else's approach; don't design your own.  There
+   are lots of partly done examples that you can find on the Web,
+   including in this tutorial.
+
+2. Generate a data subset (the first few 100k reads, for example).
+
+2. Run commands interactively on an HPC dev node until you get all of
+   the commands basically working; track all of your commands in a
+   Word document or some such.
+
+3. Once you have a set of commands that seems to work on small data,
+   write a script.  Run the script on the small data again; make sure
+   that works.
+
+4. Turn it into a qsub script (making sure you're in the right 
+   directory, have the modules loaded, etc.)
+
+5. Make sure the qsub script works on your same small data.
+
+6. Scale up to a big test data set.
+
+7. Once that's all working, SAVE THE SCRIPT SOMEWHERE.  Then,
+   edit it to work on all your data sets (you may want to make subsets
+   again, as much as possible).
+
+8. Provide your scripts and raw counts files as part of any publication
+   or thesis, perhaps via `figshare <http://figshare.com>`__.
+
+
+Informational resources
+----------------------
+
+`UT (Austin) Sequencing Core prices <https://wikis.utexas.edu/display/GSAF/Pricing>`__ - costs and yields for sequencing.
+
+`ANGUS - summer NGS course <http://angus.readthedocs.org/en/2014>`__ - lots of resources and materials and book reference
+
+`Data Carpentry <http://www.datacarpentry.org>`__ - intro to R, etc.
+
+`Software Carpentry <http://software-carpentry.org>`__ - more scripting, Python, etc.
+
+
+Places to share data, scripts, and results files
+------------------------------------------------
+
+`Figshare <http://figshare.com>`__.

refTrans/m-count.txt

@@ -0,0 +1,72 @@
+Counting the reads that map to each feature
+===========================================
+
+After mapping of reads back to the transcriptome. Now it is time to figure out 
+the expression level of each expression unit. There is two main flavours of this process:
+
+.. toctree::
+   :maxdepth: 2
+
+   m-htseq
+   m-cufflinks
+
+
+----
+
+:doc:`./m-htseq`
+
+Simply count all the reads that map **uniquelly** to all exons of the gene. There are several 
+software packages that can do this e.g. 
+`htseq-count <http://www-huber.embl.de/users/anders/HTSeq/doc/count.html>`__.
+and `Featurecount <http://bioinf.wehi.edu.au/featureCounts/>`__.
+Both programs require the user to provide a gene model (Gene annotation file; GTF).
+Annotaions file of model oraganisms are avaliable on all gene banks (e.g. 
+`NCBI <http://www.ncbi.nlm.nih.gov/guide/howto/dwn-genome/>`__, 
+`UCSC <http://hgdownload.cse.ucsc.edu/downloads.html>`__, and
+`Ensemble <http://www.ensembl.org/info/data/ftp/index.html>`__). 
+Also they are avaliable through other genomic software projects e.g. 
+`Illumina iGenomes project <http://cufflinks.cbcb.umd.edu/igenomes.html>`__. 
+
+----
+
+:doc:`./m-cufflinks`
+
+Mapping the reads to transcripts (gene isoforms) rather than the whole genes will enable 
+the identification of changes in alterative splicaing machinary. Alternative splicing is believed to be 
+one of the most powerful regulatory mechanisms in severel biological conditions including 
+`cancer <http://www.ncbi.nlm.nih.gov/pubmed/24336324>`__, 
+`development <http://www.ncbi.nlm.nih.gov/pubmed/24752171>`__,
+`immmunology <http://www.ncbi.nlm.nih.gov/pubmed/23550649>`__, 
+`fertility <http://www.ncbi.nlm.nih.gov/pubmed/25474150>`__, and many others.
+
+`Cufflinks <http://cole-trapnell-lab.github.io/cufflinks/>`__ is the canonical software 
+package to do that. But there are other softwares for analysis of differential 
+splicing events (exon skipping, intron retention, etc.) 
+e.g. `MATS <http://rnaseq-mats.sourceforge.net/>`__, and 
+`MISO <http://genes.mit.edu/burgelab/miso/>`__.
+
+Cufflinks also 
+pools the results of all transcripts per gene to give an answer similar to the one obtained from 
+htseq-like packages. One comparative feature that should be considered is the tendency of Cufflinks to 
+utilize ambiguously aligned reads (versus htseq wich use anly the uniquelly mapped reads). 
+This feature might explain a lot of the differences in the output of the 2 softwares.
+
+Cufflinks can do the mapping based on the avaliable annoataion files which will 
+make Cufflinks operate on the known transcripts only. However Cufflinks can make 
+its own gene model. This feature has two inportant applications:
+* Prediction of denovo transcripts of known genes which is a valid possibility even in well annoatated oraganisms.  
+* Prediction of novel genes. This feature is very helpful for semi-model organisms that does not have a well annoatated transcriptome.
+On the other hand, mapping to transcripts requires the implementation of error prone modeling 
+approaches and thus requires rigorous laboratory confirmation. 
+
+
+
+----
+
+Links:
+
+* `Trapnell et al. 2013 Differential analysis of gene regulation at transcript resolution with RNA-seq <http://www.nature.com/nbt/journal/v31/n1/full/nbt.2450.html>`__
+
+----
+
+Next: :doc:`m-m-htseq`

refTrans/m-cuffdiff.txt

@@ -0,0 +1,56 @@
+Differential expression by Cuffdiff
+===================================
+
+.. shell start
+
+We'll be using `Cufdiff
+<http://www.bioconductor.org/packages/release/bioc/html/edgeR.html>`__
+to do the differential expression analysis based on the cufflinks output.
+
+define the groups you want to compare (According to our sample naming conventions, 
+it should be exactly matching the stating name of your samples)
+::
+
+   groups=(lung salivary)
+
+Cufflinks recive replicates as comma separated lists for each condition.    
+You can write those lists yourself but we can also create these lists 
+and store them in an array.
+::
+
+   gp_replicates=()
+   for gp in "${groups[@]}"; do
+        replicates=()
+        for rep in $gp*_tophat; do
+            replicates+=("$workingPath/$rep/accepted_hits.bam")
+        done
+        gp_replicates+=($(echo ${replicates[*]} | tr ' ' ','))
+   done
+
+You can see how these comma separated lists looks like by::  
+
+   declare -p gp_replicates
+   #printf "${gp_replicates[*]}"
+
+Now we can run the Cuffdiff
+
+::
+
+   cuffdiff \
+    --output-dir "$workingPath/cuffdiff_out" \
+    --num-threads 2 \
+    --multi-read-correct \
+    --frag-bias-correct ${homo_sapiens_genome} \
+    --labels ${groups[0]},${groups[1]} \
+    $GTF_file  ${gp_replicates[0]}  ${gp_replicates[1]}
+
+On HPCC, most probably you need to schadule this job 
+
+.. shell stop
+
+----
+
+Next: :doc:`m-func-analysis`
+
+
+FOLDERS=`ls -dm ./ | tr -d ' '`