From 1c24f28616ee3bce65a1110eb8b0ccf5026a5ba1 Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Fri, 30 Jan 2015 22:27:05 -0600 Subject: [PATCH 01/14] initalizing my first commitment --- refTrans/firstCommitment.txt | Bin 0 -> 132 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 refTrans/firstCommitment.txt diff --git a/refTrans/firstCommitment.txt b/refTrans/firstCommitment.txt new file mode 100644 index 0000000000000000000000000000000000000000..59d9296100432b9afcc9613e8a4e7ab6afa4a5b1 GIT binary patch literal 132 zcmXAh!3}^Q6hnJ9aR=sb4i`X#ghasr7nfHba%lR}{@?d0@(NmIl~lFJKt`8}vKn&e oAPIQE$G!hO(R#9tp9<}8qZwNciY)b&6LPe@vFzR4M>K2F2Sx=NAOHXW literal 0 HcmV?d00001 From 16fd6f4a505baa49ea90ed51a494018b2e976d79 Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Sun, 1 Feb 2015 11:08:12 -0600 Subject: [PATCH 02/14] change copyrights and version --- conf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/conf.py b/conf.py index 5c57d90..e3996a2 100644 --- a/conf.py +++ b/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. From 8c391bbdce447fe0c9ebaafe8b4193496be38425 Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Sun, 1 Feb 2015 11:20:56 -0600 Subject: [PATCH 03/14] copy the rst files of the msu workshop for reference based RNAseq assembly and editting the indicies files --- index.txt | 13 +- refTrans/GETTING-STARTED.txt | 11 ++ refTrans/index.txt | 21 +++ refTrans/m-advice.txt | 93 ++++++++++ refTrans/m-data-analysis.txt | 108 ++++++++++++ refTrans/m-hpc-qsub.txt | 41 +++++ refTrans/m-htseq.txt | 42 +++++ refTrans/m-more-tophat.txt | 67 ++++++++ refTrans/m-quality.txt | 205 ++++++++++++++++++++++ refTrans/m-tophat.txt | 66 ++++++++ refTrans/model-orgs.txt | 21 +++ refTrans/more-resources.txt | 18 ++ refTrans/s-building-a-reference.txt | 121 +++++++++++++ refTrans/s-data-analysis.txt | 60 +++++++ refTrans/s-even-more-tophat.txt | 65 +++++++ refTrans/s-more-tophat.txt | 48 ++++++ refTrans/s-quality.txt | 252 ++++++++++++++++++++++++++++ refTrans/s-tophat-round2.txt | 80 +++++++++ refTrans/semi-orgs.txt | 23 +++ 19 files changed, 1354 insertions(+), 1 deletion(-) create mode 100644 refTrans/GETTING-STARTED.txt create mode 100644 refTrans/index.txt create mode 100644 refTrans/m-advice.txt create mode 100644 refTrans/m-data-analysis.txt create mode 100644 refTrans/m-hpc-qsub.txt create mode 100644 refTrans/m-htseq.txt create mode 100644 refTrans/m-more-tophat.txt create mode 100644 refTrans/m-quality.txt create mode 100644 refTrans/m-tophat.txt create mode 100644 refTrans/model-orgs.txt create mode 100644 refTrans/more-resources.txt create mode 100644 refTrans/s-building-a-reference.txt create mode 100644 refTrans/s-data-analysis.txt create mode 100644 refTrans/s-even-more-tophat.txt create mode 100644 refTrans/s-more-tophat.txt create mode 100644 refTrans/s-quality.txt create mode 100644 refTrans/s-tophat-round2.txt create mode 100644 refTrans/semi-orgs.txt diff --git a/index.txt b/index.txt index 1c4ad8e..bc54d8e 100644 --- a/index.txt +++ b/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 `__ +and made use of the `MSU High Performance Compute Center `__ . + + Additional information ---------------------- Need help? Either post comments on the bottom of each page, OR `sign up for the mailing list `__. - + Have you used these protocols in a scientific publication? We'll have citation instructions up soon. diff --git a/refTrans/GETTING-STARTED.txt b/refTrans/GETTING-STARTED.txt new file mode 100644 index 0000000..70c85bc --- /dev/null +++ b/refTrans/GETTING-STARTED.txt @@ -0,0 +1,11 @@ +1. Edit conf.py and replace 'labibi' in the 'project' line. + +2. Either replace or eliminate the Google Analytics ID, the disqus name, + and the github information. :: + + html_context = { + "google_analytics_id" : 'UA-36028965-1', + "disqus_shortname" : 'labibi', + "github_base_account" : 'ctb', + "github_project" : 'labibi', + } diff --git a/refTrans/index.txt b/refTrans/index.txt new file mode 100644 index 0000000..a203b42 --- /dev/null +++ b/refTrans/index.txt @@ -0,0 +1,21 @@ +=================================================== +Reference based RNAseq assembly (refTrans) protocol +=================================================== + +:author: Tamer Mansour and Titus Brown + +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 `__ +and made use of the `MSU High Performance Compute Center `__ . + +.. toctree:: + :maxdepth: 2 + + semi-orgs + model-orgs + +Reference material +------------------ +:doc:`../docs/command-line` +.. :doc:`../amazon/index` \ No newline at end of file diff --git a/refTrans/m-advice.txt b/refTrans/m-advice.txt new file mode 100644 index 0000000..3bf6f32 --- /dev/null +++ b/refTrans/m-advice.txt @@ -0,0 +1,93 @@ +Miscellaneous advice +==================== + +Sequencing depth and number of samples +-------------------------------------- + +`Hart et al. (2013) +`__ 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 +`__ +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. + +Downloading your data +--------------------- + +If you do your sequencing at the MSU Core Facility, you'll get an e-mail +from them when you're samples are ready. The e-mail will give you an +FTP site, a username, and a password, as well as a URL. You can +use these to download your data. For example, if you get:: + + hostname: titan.bch.msu.edu + username: rnaseqmodel + password: QecheJa6 + + URI: ftp://rnaseqmodel:QecheJa6@titan.bch.msu.edu + +you can go to ftp://rnaseqmodel:QecheJa6@titan.bch.msu.edu in your +Web browser; that is, it lets you combine your username and password +to open that link. + +In this case, you will see a 'testdata' directory. If you click on that, +you'll see a bunch of fastq.gz files. These are the files that you +want to get onto the HPC. + +To download these files onto the HPC, log into the HPC, go to the +directory on the HPC you want to put the files in, and run a 'wget' -- +for example, on the HPC:: + + mkdir ~/testdata + cd ~/testdata + + wget -r -np -nH ftp://rnaseqmodel:QecheJa6@titan.bch.msu.edu/testdata/ + +This will download _all_ of the files in that directory. You can also do +them one at a time, e.g. to get 'Ath_Mut_1_R1.fastq.gz', you would do + + wget ftp://rnaseqmodel:QecheJa6@titan.bch.msu.edu/testdata/Ath_Mut_1_R1.fastq.gz + +Tada! + +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 `__. + +Next: :doc:`more-resources` diff --git a/refTrans/m-data-analysis.txt b/refTrans/m-data-analysis.txt new file mode 100644 index 0000000..c249c47 --- /dev/null +++ b/refTrans/m-data-analysis.txt @@ -0,0 +1,108 @@ +Data analysis & differential expression +======================================= + +We'll be using `edgeR +`__ +to do the basic differential expression analysis of our counts. + +To run edgeR, you need to write a data loading and manipulation script +in R. In this case, I've provided one -- `lung_saliva.R +`__. +This script will load in two samples with two replicates, execute an +MA plot, and provide a spreadsheet with differential expression +information in it. To download it, `click here +`__. + +Links: + +* `False Discovery Rate `__ +* `Learn R with Swirl `__ +* `Data Carpentry `__ + +Running edgeR on a data subset +------------------------------ + +To run the script on the HPC, download the script and data files for +the 100k read subsets. At the command line, do:: + + cd ~/rnaseq + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/lung_saliva.R + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/lung_repl1_counts.txt + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/lung_repl2_counts.txt + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/salivary_repl1_counts.txt + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/salivary_repl2_counts.txt + +Note: the ``saliva_repl1_counts.txt`` and ``lung_repl1_counts.txt`` +are the files produced by :doc:`m-tophat` and :doc:`m-more-tophat`, +respectively. I've also run them on the replicate data sets, which +produced the *repl2_counts.txt files. + +Next, to run the R script, do:: + + module load R + Rscript lung_saliva.R + +This will produce two files, `edgeR-MA-plot.pdf +`__ +and `edgeR-lung-vs-salivary.csv +`__; +they will be in your ``rnaseq`` folder in your home directory +on the HPC. The CSV file can be opened directly in Excel; you can +also look at it `here +`__. +It consists of five columns: gene name, log fold change, P-value, and +FDR-adjusted P-value. + +Links: + +* `edgeR tutorial from UT Austin `__ + +Functional and network analysis on differentially expressed genes +----------------------------------------------------------------- + +There are a number of sites that let you analyze gene lists, including +`DAVID `__. I have gotten a number of +recommendations for `PANTHER `__. PANTHER +lets you upload gene lists and explore their functional categories +interactively or in a gene list/annotation format. + +More specifically, you can explore your RNAseq data with + +* functional classifications, in pie chart or in list; +* tests for statistical overrepresentation; +* tests for statistical enrichment based on associated fold change. + +You need to crop and transform the data a little bit before using the +functional classification. The steps are: + +1. Download `the CSV file `__. If you've produced your own, copy it over from the HPC. + +2. Open it in Excel. + +3. Choose an FDR cutoff (suggest FDR < 0.05 or lower) and delete all the rows after that (or, copy the rows into a new spreadsheet -- might be quicker). + +4. Save as a "Tab-delimited text." (Note, on Mac OS X, you may need to save this as "Windows formatted text" instead.) + +This is now a file you can upload to PANTHER. + +(`You can download my copy of this here `__) + +To do the analysis, go to http://www.pantherdb.org/. In box 1, select "Choose file" +and find the CSV file you want to upload to PANTHER. Nothing else in box +one should be changed. + +In box 2, select "Homo sapiens." + +In box 3, select either "Functional analysis classification viewed in +gene list" or "Functional analysis classification viewed in pie +chart." + +Click submit. + +Now you can explore the results! + +Links: + +* `PANTHER database `__ + +Next: :doc:`m-hpc-qsub` diff --git a/refTrans/m-hpc-qsub.txt b/refTrans/m-hpc-qsub.txt new file mode 100644 index 0000000..4fd0d51 --- /dev/null +++ b/refTrans/m-hpc-qsub.txt @@ -0,0 +1,41 @@ +Submitting jobs to the MSU HPC queue +==================================== + +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'. + +Here are four scripts (plus a common library script) that you could use +to run some analyses on the HPC. To run, do something like this:: + + module load powertools + getexample RNAseq-model + + mkdir ~/rnaseq + mkdir ~/rnaseq/script + cd ~/rnaseq/script + ln -fs ~/RNAseq-model/data/*.fastq . + + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/env.sh + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/process-1.sh + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/process-2.sh + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/process-3.sh + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/process-4.sh + + qsub process-1.sh + qsub process-2.sh + qsub process-3.sh + qsub process-4.sh + +You can use 'qstat | grep ' to check on your jobs' status. + +Next: :doc:`m-advice` diff --git a/refTrans/m-htseq.txt b/refTrans/m-htseq.txt new file mode 100644 index 0000000..cc2f732 --- /dev/null +++ b/refTrans/m-htseq.txt @@ -0,0 +1,42 @@ +Counting the reads that map to each gene +======================================== + +Now that we know which reads go with which gene, we'll use +`htseq-count `__. + +First, load the PySAM and HTSeq software packages:: + + module load PySAM/0.6 + module load HTSeq/0.6.1 + +And next, run HTSeq:: + + htseq-count --format=bam --stranded=no --order=pos \ + tophat_salivary_repl1/accepted_hits.bam \ + ~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Annotation/Genes/genes.gtf > salivary_repl1_counts.txt + +When this is done, type:: + + less salivary_repl1_counts.txt + +(again, use 'q' to exit). These are your gene counts. + +Note, these are *raw* gene counts - the number of reads that map to +each feature (gene, in this case). They are not normalized by length +of gene. According to `this post on seqanswers +`__, both +DEseq and edgeR want exactly this kind of information! + +Questions: + +* what are the 'ENS*...' names? +* what do these parameters mean? +* what parameters does HTSeq take? +* why are we using so many programs? + +Links: + +* `HTSeq `__ +* `htseq-count documentation `__ + +Next: :doc:`m-more-tophat` diff --git a/refTrans/m-more-tophat.txt b/refTrans/m-more-tophat.txt new file mode 100644 index 0000000..ea518af --- /dev/null +++ b/refTrans/m-more-tophat.txt @@ -0,0 +1,67 @@ +Processing another sample with TopHat and HTSeq +=============================================== + +Let's script this! + +Type:: + + cd ~/rnaseq/ + module load Trimmomatic/0.32 + cat > lung-tophat.sh < lung_repl1_R1.fq.gz + head -400000 ~/RNAseq-model/data/ERR315326_2.fastq | gzip > lung_repl1_R2.fq.gz + + # run Trimmomatic. Here, the inputs are 'lung_repl1_R1.fq.gz' and + # 'lung_repl1_R2.fq.gz', and the outputs are 'lung_repl1_R1.qc.fq.gz' + # and 'lung_repl1_R2.qc.fq.gz'. + java -jar \$TRIM/trimmomatic PE lung_repl1_R1.fq.gz lung_repl1_R2.fq.gz \ + lung_repl1_R1.qc.fq.gz s1_se lung_repl1_R2.qc.fq.gz s2_se \ + ILLUMINACLIP:\$TRIM/adapters/TruSeq2-PE.fa:2:40:15 \ + LEADING:2 TRAILING:2 \ + SLIDINGWINDOW:4:2 \ + MINLEN:25 + + # now run Tophat! + # The inputs are the outputs of the previous Trimmomatic step. + # The outputs are going to be under the 'tophat_lung_repl1' directory. + tophat -p 4 \ + -G ~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Annotation/Genes/genes.gtf \ + --transcriptome-index=$HOME/RNAseq-model/transcriptome \ + -o tophat_lung_repl1 \ + ~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Sequence/Bowtie2Index/genome \ + lung_repl1_R1.qc.fq.gz lung_repl1_R2.qc.fq.gz + + # count the hits by gene -- 'tophat_lung_repl1' is the main output, + # from Tophat. + htseq-count --format=bam --stranded=no --order=pos tophat_lung_repl1/accepted_hits.bam \ + ~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Annotation/Genes/genes.gtf > lung_repl1_counts.txt + EOF + +(Be sure to press the Enter or Return key after pasting this in!) This is +`called a 'heredoc' `__ +and it gives a way to write a shell script via copy-paste ;). + +Next, type:: + + bash lung-tophat.sh + +This will run all of the commands in the file 'lung-tophat.sh'. + +You can use the 'nano' editor to modify this file -- type:: + + nano lung-tophat.sh + +Next: :doc:`m-data-analysis` diff --git a/refTrans/m-quality.txt b/refTrans/m-quality.txt new file mode 100644 index 0000000..d9b7aa7 --- /dev/null +++ b/refTrans/m-quality.txt @@ -0,0 +1,205 @@ +Short read quality and trimming +=============================== + +Log into the HPC with SSH; use your MSU NetID and log into the machine +'gateway.hpcc.msu.edu'. There copy/paste:: + + module load powertools + getexample RNAseq-model + +and you should see something about linking. + +0. Getting the data +------------------- + +We'll be using a few RNAseq data sets from Fagerberg et al., `Analysis +of the Human Tissue-specific Expression by Genome-wide Integration of +Transcriptomics and Antibody-based Proteomics +`__. + +You can get this data from the European Nucleotide Archive under +ERP003613 -- go to http://www.ebi.ac.uk/ena/data/view/ERP003613 +to see all the samples. + +I picked two samples: `salivary gland +`__ and `lung +`__. Note that each +sample has two replicates, and each replicate has two files. + +**Don't download them**, but if you were downloading these yourself, +you would want the "Fastq files (ftp)", both File 1 and File 2. (They +take a few hours to download!) + +We've already loaded the data onto the MSU HPC, and you've loaded +it with 'module load powertools'. + +To log into a compute node, type:: + + ~/RNAseq-model/login.sh + +If this doesn't work, do:: + + ssh dev-intel14-phi + +Now do:: + + ls -la ~/RNAseq-model/data/ + +You'll see something like :: + + -r--r--r-- 1 mscholz common-data 3714262571 Dec 4 08:44 ERR315325_1.fastq + -r--r--r-- 1 mscholz common-data 3714262571 Dec 4 08:44 ERR315325_2.fastq + -r--r--r-- 1 mscholz common-data 2365633645 Dec 4 08:44 ERR315326_1.fastq + -r--r--r-- 1 mscholz common-data 2365633645 Dec 4 08:44 ERR315326_2.fastq + +which tells you that this file is 900,000,000 bytes or about 900 MB. +Quite large! + +1. Copying in some data to work with. +------------------------------------- + +First, make a directory:: + + mkdir ~/rnaseq + cd ~/rnaseq + +Copy in a subset of the data (100,000 reads):: + + head -400000 ~/RNAseq-model/data/ERR315325_1.fastq | gzip > salivary_repl1_R1.fq.gz + head -400000 ~/RNAseq-model/data/ERR315325_2.fastq | gzip > salivary_repl1_R2.fq.gz + +These are FASTQ files -- let's take a look:: + + less salivary_repl1_R1.fq.gz + +(type 'q' to exit less) + +Question: + +* why are some files named ERR*? +* why are some files named salivary*? +* why is there R1 and R2 in the name? + +Links: + +* `FASTQ Format `__ + +2. FastQC +--------- + +We're going to use `FastQC `__ to summarize the data. + +First, we need to load the FastQC software into our account:: + + module load FastQC/0.11.2 + +(You have to do this each time you log in and want to use FastQC.) + +Now, run FastQC on both of the salivary gland files:: + + fastqc salivary_repl1_R1.fq.gz + fastqc salivary_repl1_R2.fq.gz + +Now type 'ls':: + + ls + +and you will see :: + + salivary_repl1_R1_fastqc.html + salivary_repl1_R1_fastqc.zip + salivary_repl1_R2_fastqc.html + salivary_repl1_R2_fastqc.zip + +Copy these to your laptop and open them in a browser. If you're on a +Mac or Linux machine, you can type:: + + scp username@hpc.msu.edu:salivary*fastqc.* /tmp + +and then open the html files in your browser. For Windows, we'll +figure it out ;). + +You can also view my versions: `salivary_repl1_R1_fastqc.html +`__ +and `salivary_repl1_R2_fastqc.html +`__ + +Questions: + +* What should you pay attention to in the FastQC report? +* Which is "better", R1 or R2? + +Links: + +* `FastQC `__ +* `FastQC tutorial video `__ + +3. Trimmomatic +-------------- + +Now we're going to do some trimming! We'll be using +`Trimmomatic `__. +For +a discussion of optimal RNAseq trimming strategies, see `MacManes, 2014 `__. + +First, load the Trimmomatic software:: + + module load Trimmomatic/0.32 + +Next, run Trimmomatic:: + + java -jar $TRIM/trimmomatic PE salivary_repl1_R1.fq.gz salivary_repl1_R2.fq.gz\ + salivary_repl1_R1.qc.fq.gz s1_se salivary_repl1_R2.qc.fq.gz s2_se \ + ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 \ + LEADING:2 TRAILING:2 \ + SLIDINGWINDOW:4:2 \ + MINLEN:25 + +You should see output that looks like this:: + + ... + Input Read Pairs: 100000 Both Surviving: 95236 (95.24%) Forward Only Surviving: 4764 (4.76%) Reverse Only Surviving: 0 (0.00%) Dropped: 0 (0.00%) + TrimmomaticPE: Completed successfully + +Questions: + +* How do you figure out what the parameters mean? +* How do you figure out what parameters to use? +* What adapters do you use? +* What version of Trimmomatic are we using here? (And FastQC?) +* Are parameters different for RNAseq and genomic? +* What's with these annoyingly long and complicated filenames? +* What do we do with the single-ended files (s1_se and s2_se?) + +Links: + +* `Trimmomatic `__ + +4. FastQC again +--------------- + +Run FastQC again:: + + fastqc salivary_repl1_R1.qc.fq.gz + fastqc salivary_repl1_R2.qc.fq.gz + +(Note that you don't need to load the module again.) + +Copy them to your laptop and open them, OR you can view mine: `salivary_repl1_R1.qc_fastqc.html +`__ +and `salivary_repl1_R2.qc_fastqc.html +`__ + +Let's take a look at the output files:: + + less salivary_repl1_R1.qc.fq.gz + +(again, use 'q' to exit less). + +Questions: + +* Why are some of the reads shorter than others? +* is the quality trimmed data "better" than before? +* Does it matter that you still have adapters!? + +Next: :doc:`m-tophat` diff --git a/refTrans/m-tophat.txt b/refTrans/m-tophat.txt new file mode 100644 index 0000000..e3bb17b --- /dev/null +++ b/refTrans/m-tophat.txt @@ -0,0 +1,66 @@ +Mapping reads to the transcriptome with TopHat +============================================== + +Now that we have some quality-controlled reads, we're going to *map* the +reads to the reference gene set, for the purpose of counting how many +reads have come from each gene. We'll be using the `TopHat software +`__ + +For this purpose, we've already installed the human reference gene set +on the HPC (as part of the data you loaded at the beginning). In this +case we've loaded in the `Illumina iGenomes project +`__ into the RNAseq-model +data set. + +Mapping reads +------------- + +Load the TopHat software:: + + module load TopHat2/2.0.12 + +And now run TopHat:: + + cd ~/rnaseq + tophat -p 4 \ + -G ~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Annotation/Genes/genes.gtf \ + --transcriptome-index=$HOME/RNAseq-model/transcriptome \ + -o tophat_salivary_repl1 \ + ~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Sequence/Bowtie2Index/genome \ + salivary_repl1_R1.qc.fq.gz salivary_repl1_R2.qc.fq.gz + +This will take about 15 minutes. + +Questions: + +* What are all these parameters?! +* How do we pick the transcriptome/genome? +* Why is it so slow? +* What is different about mapping RNAseq reads vs mapping genomic reads. + +Links: + +* `TopHat manual `__ +* `Illumina iGenomes project `__ + +Counting mapped reads percentage +-------------------------------- + +Let's ask samtools for the total number of reads that mapped:: + + samtools view -c -F 4 tophat_salivary_repl1/accepted_hits.bam + +You should get around 179,312. + +If we look at the `FastQC report +`__, +we can see (at the top) that the total number of reads in the R1 file +is 95,236 (which is the same number of reads as in the R2 file, because +we only looked at the paired reads that came out of Trimmomatic). + +In total, that's 190,472 reads in your trimmed (QC) data, and 179,312 +mapped -- about 94%. That's pretty good! + +---- + +Next: :doc:`m-htseq` diff --git a/refTrans/model-orgs.txt b/refTrans/model-orgs.txt new file mode 100644 index 0000000..a760bf4 --- /dev/null +++ b/refTrans/model-orgs.txt @@ -0,0 +1,21 @@ +2014 / December / mRNAseq on model organisms +============================================ + +This workshop was given on December 4th and 5th, 2014, by +C. Titus Brown and Matt Scholz. + +.. figure:: _static/model-rnaseq-pipeline.png + +Tutorials: + +.. toctree:: + :maxdepth: 2 + + m-quality + m-tophat + m-htseq + m-more-tophat + m-data-analysis + m-hpc-qsub + m-advice + more-resources diff --git a/refTrans/more-resources.txt b/refTrans/more-resources.txt new file mode 100644 index 0000000..c85627d --- /dev/null +++ b/refTrans/more-resources.txt @@ -0,0 +1,18 @@ +More resources +============== + +Informational resources +---------------------- + +`UT (Austin) Sequencing Core prices `__ - costs and yields for sequencing. + +`ANGUS - summer NGS course `__ - lots of resources and materials and book reference + +`Data Carpentry `__ - intro to R, etc. + +`Software Carpentry `__ - more scripting, Python, etc. + +Places to share data, scripts, and results files +------------------------------------------------ + +`Figshare `__. diff --git a/refTrans/s-building-a-reference.txt b/refTrans/s-building-a-reference.txt new file mode 100644 index 0000000..ae8f9e1 --- /dev/null +++ b/refTrans/s-building-a-reference.txt @@ -0,0 +1,121 @@ +Building a new reference transcriptome +====================================== + +The distinguishing feature of (what I call) semi-model organisms is that +while they may have a decent genome reference, their transcriptome +annotation is poor. There can be several reasons for this, but generally +it boils down to lack of resources and/or attention -- it takes a *lot* +of effort to build a high quality transcriptome! + +For this purpose, we've already installed the chicken reference genome +set on the HPC (as part of the data you loaded at the beginning). In +this case we've loaded in the `Illumina iGenomes project +`__ into the RNAseq-semimodel +location. + +See paper: + + http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3334321/ + +Map all the reads to the genome with TopHat +------------------------------------------- + +.. @@ add links etc. + +We'll be using the `TopHat software +`__. + +Do:: + + module load TopHat2/2.0.12 + + tophat -p 4 \ + -G $HOME/RNAseq-semimodel/reference/Gallus_gallus/UCSC/galGal3/Annotation/Genes/genes.gtf \ + --transcriptome-index=$HOME/RNAseq-semimodel/tophat/transcriptome \ + -o tophat_all \ + $HOME/RNAseq-semimodel/reference/Gallus_gallus/UCSC/galGal3/Sequence/Bowtie2Index/genome \ + female_repl1_R1.qc.fq.gz,male_repl1_R1.qc.fq.gz,female_repl2_R1.qc.fq.gz,male_repl2_R1.qc.fq.gz \ + female_repl1_R2.qc.fq.gz,male_repl1_R2.qc.fq.gz,female_repl2_R2.qc.fq.gz,male_repl2_R2.qc.fq.gz + +Questions: + +* What are all these parameters?! +* How do we pick the transcriptome/genome? +* Why is it so slow? +* What is different about mapping RNAseq reads vs mapping genomic reads? + +Links: + +* `TopHat manual `__ +* `Illumina iGenomes project `__ + +Evaluating the mapping +---------------------- + +Check out the details:: + + less tophat_all/align_summary.txt + +Build a new transcriptome ("ab initio") from the combined reads using Cufflinks +------------------------------------------------------------------------------- + +Now that we've mapped the reads, let's put them all together into exons +and gene models:: + + module load cufflinks/2.2.1 + + cufflinks -o cuff_all tophat_all/accepted_hits.bam + +.. @@ cufflinks diagram + +Questions: + +* What exactly is Cufflinks doing? + +Merge the new transcriptome with the existing reference transcriptome +--------------------------------------------------------------------- + +We already have some decent gene models; let's merge our new and the old ones:: + + ls -1 cuff_all/transcripts.gtf > cuff_list.txt + + cuffmerge -g $HOME/RNAseq-semimodel/reference/Gallus_gallus/UCSC/galGal3/Annotation/Genes/genes.gtf \ + -o cuffmerge_all \ + -s $HOME/RNAseq-semimodel/reference/Gallus_gallus/UCSC/galGal3/Sequence/WholeGenomeFasta/genome.fa \ + cuff_list.txt + +Do some cleanup:: + + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/remove-nostrand.py + python remove-nostrand.py cuffmerge_all/merged.gtf > cuffmerge_all/nostrand.gtf + +Questions: + +* why do you want to merge? +* why would you have a list of more than one thing in list.txt? +* come to think of it, why aren't you (re)mapping all your reads every time? +* what's with the 'remove nostrand' script? + +Extracting your new transcriptome sequences +------------------------------------------- + +To get a look at the actual DNA sequences, do:: + + gffread -w cuffmerge_all.fa \ + -g $HOME/RNAseq-semimodel/reference/Gallus_gallus/UCSC/galGal3/Sequence/WholeGenomeFasta/genome.fa \ + cuffmerge_all/nostrand.gtf + +Questions: + +* What's the difference between a GTF file and the FA file? + +Checking out your new transcriptome +----------------------------------- + +Take a look at the top of your FASTA file:: + + head -30 cuffmerge_all.fa + +Head on over to `the chicken genome browser `__ and try BLATing the sequence! + +:doc:`s-tophat-round2` diff --git a/refTrans/s-data-analysis.txt b/refTrans/s-data-analysis.txt new file mode 100644 index 0000000..9928247 --- /dev/null +++ b/refTrans/s-data-analysis.txt @@ -0,0 +1,60 @@ +Data analysis & differential expression +======================================= + +We'll be using `edgeR +`__ +to do the basic differential expression analysis of our counts. + +To run edgeR, you need to write a data loading and manipulation script +in R. In this case, I've provided one -- `chick.R +`__. +This script will load in two samples with two replicates, execute an +MA plot, do an MDS analysis/plot, and provide a spreadsheet with +differential expression information in it. To download it, `click +here +`__. + +Links: + +* `False Discovery Rate `__ +* `Learn R with Swirl `__ +* `Data Carpentry `__ + +Running edgeR on a data subset +------------------------------ + +To run the script on the HPC, download the script and data files for +the 100k read subsets. At the command line, do:: + + cd ~/rnaseq + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/chick.R + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/chick-subset/female_repl1_counts.txt + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/chick-subset/female_repl2_counts.txt + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/chick-subset/male_repl1_counts.txt + curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/chick-subset/male_repl2_counts.txt + +These various .txt files are produced by :doc:`s-tophat-round2`, :doc:`s-more-tophat`, and :doc:`s-even-more-tophat`. + +Next, to run the R script, do:: + + module load R + Rscript chick.R + +This will produce three files, `chick-edgeR-MA-plot.pdf +`__ +`chick-edgeR-MDS.pdf +`__ +and `chick-edgeR.csv +`__; +they will be in your ``rnaseq`` folder in your home directory +on the HPC. The CSV file can be opened directly in Excel; you can +also look at it `here +`__. +It consists of five columns: gene name, log fold change, P-value, and +FDR-adjusted P-value. + +Links: + +* `edgeR tutorial from UT Austin `__ + +Next: :doc:`m-advice` diff --git a/refTrans/s-even-more-tophat.txt b/refTrans/s-even-more-tophat.txt new file mode 100644 index 0000000..7c8a591 --- /dev/null +++ b/refTrans/s-even-more-tophat.txt @@ -0,0 +1,65 @@ +Using the HPC to run jobs in parallel +===================================== + +Type:: + + cat > chick-tophat-3.sh < male_repl1_counts.txt + + EOF + +Next, do it for male_repl2:: + + cat > chick-tophat-4.sh < male_repl2_counts.txt + + EOF + +To submit these to the queue, do:: + + qsub chick-tophat-3.sh + qsub chick-tophat-4.sh + +Next: doc:`s-data-analysis` diff --git a/refTrans/s-more-tophat.txt b/refTrans/s-more-tophat.txt new file mode 100644 index 0000000..2aadb26 --- /dev/null +++ b/refTrans/s-more-tophat.txt @@ -0,0 +1,48 @@ +Processing another sample with TopHat and HTSeq +=============================================== + +Let's script this! + +Type:: + + cd ~/rnaseq/ + cat > chick-tophat-2.sh < female_repl2_counts.txt + + EOF + +(Be sure to press the Enter or Return key after pasting this in!) This is +`called a 'heredoc' `__ +and it gives a way to write a shell script via copy-paste ;). + +Next, type:: + + bash chick-tophat-2.sh + +This will run all of the commands in the file 'chick-tophat-2.sh'. + +You can use the 'nano' editor to modify this file -- type:: + + nano chick-tophat-2.sh + +:doc:`s-even-more-tophat` diff --git a/refTrans/s-quality.txt b/refTrans/s-quality.txt new file mode 100644 index 0000000..5e9c35d --- /dev/null +++ b/refTrans/s-quality.txt @@ -0,0 +1,252 @@ +Short read quality and trimming +=============================== + +.. note:: + + Reminder: if you're on Windows, you should install `mobaxterm `__. + +Log into the HPC with SSH; use your MSU NetID and log into the machine +'hpc.msu.edu'. There copy/paste:: + + cd + module load powertools + getexample RNAseq-semimodel + +This will put all of the example files for today in your home directory +under the directory 'RNAseq-semimodel'. + +0. Getting the data +------------------- + +http://genomebiology.com/content/14/3/R26 + +http://www.ebi.ac.uk/ena/data/view/SRA055442 + +http://www.ebi.ac.uk/ena/data/view/SAMN01096082 +http://www.ebi.ac.uk/ena/data/view/SAMN01096083 + +http://www.ebi.ac.uk/ena/data/view/SAMN01096084 +http://www.ebi.ac.uk/ena/data/view/SAMN01096085 + +Note that each +sample has two replicates, and each replicate has two files. + +**Don't download them**, but if you were downloading these yourself, +you would want the "Fastq files (ftp)", both File 1 and File 2. (They +take a few hours to download!) + +We've already loaded the data onto the MSU HPC, and you've loaded +it with 'module load powertools'. + +To log into a compute node, type:: + + ~/RNAseq-semimodel/login.sh + +If this doesn't work, do:: + + ssh dev-intel14-phi + +Now do:: + + ls -l ~/RNAseq-semimodel/data/ + +You'll see something like :: + + -r--r--r-- 1 mscholz common-data 6781517200 Dec 9 09:46 SRR534005_1.fastq.gz + -r--r--r-- 1 mscholz common-data 7023515467 Dec 9 09:50 SRR534005_2.fastq.gz + -r--r--r-- 1 mscholz common-data 7285848617 Dec 9 09:41 SRR534006_1.fastq.gz + -r--r--r-- 1 mscholz common-data 7542383700 Dec 9 09:43 SRR534006_2.fastq.gz + -r--r--r-- 1 mscholz common-data 7219923066 Dec 9 09:47 SRR536786_1.fastq.gz + -r--r--r-- 1 mscholz common-data 7467116873 Dec 9 09:49 SRR536786_2.fastq.gz + -r--r--r-- 1 mscholz common-data 7694614208 Dec 9 09:40 SRR536787_1.fastq.gz + -r--r--r-- 1 mscholz common-data 7944043814 Dec 9 09:44 SRR536787_2.fastq.gz + +These files are each approximately 7-8 GB in size! + +1. Copying in some data to work with. +------------------------------------- + +First, make a directory:: + + mkdir ~/rnaseq + cd ~/rnaseq + +Copy in a subset of the data (100,000 reads):: + + gunzip -c ~/RNAseq-semimodel/data/SRR534005_1.fastq.gz | head -400000 | gzip > female_repl1_R1.fq.gz + gunzip -c ~/RNAseq-semimodel/data/SRR534005_2.fastq.gz | head -400000 | gzip > female_repl1_R2.fq.gz + +These are FASTQ files -- let's take a look:: + + less female_repl1_R1.fq.gz + +(type 'q' to exit less) + +Question: + +* why are some files named SRR*? +* why are some files named female*? +* why are there R1 and R2 in the name? + +Links: + +* `FASTQ Format `__ + +2. FastQC +--------- + +We're going to use `FastQC `__ to summarize the data. + +First, we need to load the FastQC software into our account:: + + module load FastQC/0.11.2 + +(You have to do this each time you log in and want to use FastQC.) + +Now, run FastQC on both of the female files:: + + fastqc female_repl1_R1.fq.gz + fastqc female_repl1_R2.fq.gz + +Now type 'ls':: + + ls + +and you will see :: + + female_repl1_R1_fastqc.html + female_repl1_R1_fastqc.zip + female_repl1_R2_fastqc.html + female_repl1_R2_fastqc.zip + +Copy these to your laptop and open them in a browser. If you're on a +Mac or Linux machine, you can type:: + + scp username@hpc.msu.edu:rnaseq/female*fastqc.* /tmp + +and then open the html files in your browser. For Windows, if you're using +mobaxterm, most of you should have a file transfer window on the left. +Click 'refresh' (green circle icon fourth from the left) and then navigate +into the 'rnaseq' folder; you should see the 'female_repl...' files there. +Drag and drop those onto your Windows machine. + +You can also view my versions: `female_repl1_R1_fastqc.html +`__ +and `female_repl1_R2_fastqc.html +`__ + +Questions: + +* What should you pay attention to in the FastQC report? +* Which is "better", R1 or R2? + +Links: + +* `FastQC `__ +* `FastQC tutorial video `__ + +3. Trimmomatic +-------------- + +Now we're going to do some trimming! We'll be using +`Trimmomatic `__. For +a discussion of optimal RNAseq trimming strategies, see `MacManes, 2014 `__. + +First, load the Trimmomatic software:: + + module load Trimmomatic/0.32 + +Next, run Trimmomatic:: + + java -jar $TRIM/trimmomatic PE female_repl1_R1.fq.gz female_repl1_R2.fq.gz\ + female_repl1_R1.qc.fq.gz s1_se female_repl1_R2.qc.fq.gz s2_se \ + ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 \ + LEADING:2 TRAILING:2 \ + SLIDINGWINDOW:4:2 \ + MINLEN:25 + +You should see output that looks like this:: + + ... + Quality encoding detected as phred33 + Input Read Pairs: 100000 Both Surviving: 95583 (95.58%) Forward Only Surviving: 4262 (4.26%) Reverse Only Surviving: 86 (0.09%) Dropped: 69 (0.07%) + ... + +Questions: + +* How do you figure out what the parameters mean? +* How do you figure out what parameters to use? +* What adapters do you use? +* What version of Trimmomatic are we using here? (And FastQC?) +* Are parameters different for RNAseq and genomic? +* What's with these annoyingly long and complicated filenames? +* What do we do with the single-ended files (s1_se and s2_se?) + +Links: + +* `Trimmomatic `__ + +4. FastQC again +--------------- + +Run FastQC again:: + + fastqc female_repl1_R1.qc.fq.gz + fastqc female_repl1_R2.qc.fq.gz + +(Note that you don't need to load the module again.) + +Copy them to your laptop and open them, OR you can view mine: `female_repl1_R1.qc_fastqc.html +`__ +and `female_repl1_R2.qc_fastqc.html +`__ + +Let's take a look at the output files:: + + less female_repl1_R1.qc.fq.gz + +(again, use 'q' to exit less). + +Questions: + +* Why are some of the reads shorter than others? +* is the quality trimmed data "better" than before? +* Does it matter that you still have adapters!? + +5. Subset and trim the rest of the sequences +-------------------------------------------- + +Copy and paste all of the below at once:: + + gunzip -c ~/RNAseq-semimodel/data/SRR534006_1.fastq.gz | head -400000 | gzip > female_repl2_R1.fq.gz + gunzip -c ~/RNAseq-semimodel/data/SRR534006_2.fastq.gz | head -400000 | gzip > female_repl2_R2.fq.gz + + gunzip -c ~/RNAseq-semimodel/data/SRR536786_1.fastq.gz | head -400000 | gzip > male_repl1_R1.fq.gz + gunzip -c ~/RNAseq-semimodel/data/SRR536786_2.fastq.gz | head -400000 | gzip > male_repl1_R2.fq.gz + + gunzip -c ~/RNAseq-semimodel/data/SRR536787_1.fastq.gz | head -400000 | gzip > male_repl2_R1.fq.gz + gunzip -c ~/RNAseq-semimodel/data/SRR536787_2.fastq.gz | head -400000 | gzip > male_repl2_R2.fq.gz + + java -jar $TRIM/trimmomatic PE female_repl2_R1.fq.gz female_repl2_R2.fq.gz\ + female_repl2_R1.qc.fq.gz s1_se female_repl2_R2.qc.fq.gz s2_se \ + ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 \ + LEADING:2 TRAILING:2 \ + SLIDINGWINDOW:4:2 \ + MINLEN:25 + + java -jar $TRIM/trimmomatic PE male_repl1_R1.fq.gz male_repl1_R2.fq.gz\ + male_repl1_R1.qc.fq.gz s1_se male_repl1_R2.qc.fq.gz s2_se \ + ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 \ + LEADING:2 TRAILING:2 \ + SLIDINGWINDOW:4:2 \ + MINLEN:25 + + java -jar $TRIM/trimmomatic PE male_repl2_R1.fq.gz male_repl2_R2.fq.gz\ + male_repl2_R1.qc.fq.gz s1_se male_repl2_R2.qc.fq.gz s2_se \ + ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 \ + LEADING:2 TRAILING:2 \ + SLIDINGWINDOW:4:2 \ + MINLEN:25 + + +Next: :doc:`s-building-a-reference` diff --git a/refTrans/s-tophat-round2.txt b/refTrans/s-tophat-round2.txt new file mode 100644 index 0000000..4627dcc --- /dev/null +++ b/refTrans/s-tophat-round2.txt @@ -0,0 +1,80 @@ +Mapping reads to the transcriptome with TopHat +============================================== + +Now that we have some quality-controlled reads and a new reference +transcriptome, we're going to map the +reads to the reference genome, **using the new reference transcriptome**. +We'll again be using the `TopHat software +`__ + +Mapping reads +------------- + +Load the TopHat software, if you haven't already:: + + module load TopHat2/2.0.12 + +And now run TopHat:: + + cd ~/rnaseq + tophat -p 4 \ + -G cuffmerge_all/nostrand.gtf \ + --transcriptome-index=$HOME/RNAseq-semimodel/tophat/merged \ + -o tophat_female_repl1 \ + ~/RNAseq-semimodel/reference/Gallus_gallus/UCSC/galGal3/Sequence/Bowtie2Index/genome \ + female_repl1_R1.qc.fq.gz female_repl1_R2.qc.fq.gz + +This will take about 15 minutes. + +Questions: + +* How do we pick the transcriptome/genome? + +Viewing the mapped reads percentage +----------------------------------- + +Let's look at these numbers specifically:: + + less tophat_female_repl1/align_summary.txt + +Making gene counts +------------------ + +Now that we know which reads go with which gene, we'll use +`htseq-count `__. + +First, load the PySAM and HTSeq software packages:: + + module load HTSeq/0.6.1 + +And next, run HTSeq:: + + htseq-count --format=bam --stranded=no --order=pos \ + tophat_female_repl1/accepted_hits.bam \ + cuffmerge_all/nostrand.gtf > female_repl1_counts.txt + +When this is done, type:: + + less female_repl1_counts.txt + +(again, use 'q' to exit). These are your gene counts. + +Note, these are *raw* gene counts - the number of reads that map to +each feature (gene, in this case). They are not normalized by length +of gene. According to `this post on seqanswers +`__, both +DEseq and edgeR want exactly this kind of information! + +Questions: + +* what are the 'TCONS...' names? +* what do these parameters mean? +* what parameters does HTSeq take? +* why are we using so many programs? + +Links: + +* `HTSeq `__ +* `htseq-count documentation `__ + +Next: :doc:`s-more-tophat` diff --git a/refTrans/semi-orgs.txt b/refTrans/semi-orgs.txt new file mode 100644 index 0000000..b27c55f --- /dev/null +++ b/refTrans/semi-orgs.txt @@ -0,0 +1,23 @@ +2014 / December / mRNAseq on semi-model organisms +================================================= + +This workshop was given on December 10th and 11th, 2014, by C. Titus +Brown and Matt Scholz. + +.. .. figure:: _static/model-rnaseq-pipeline.png + +What are semi-model organisms? `One with bad transcriptomes. `__! + +Tutorials: + +.. toctree:: + :maxdepth: 2 + + s-quality + s-building-a-reference.rst + s-tophat-round2.rst + s-more-tophat + s-even-more-tophat + s-data-analysis + m-advice + more-resources From b2828429b1fee3609f4cf7fb8c7254d213aa77f5 Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Mon, 2 Feb 2015 10:30:21 -0600 Subject: [PATCH 04/14] 1. create new file to get the data and load software. This file can have another version for amazon. update the m-quality and model-orgs files accordingly. 2. add the image of RNAseq analysis model. 3. apply the literate-resting roles to these files --- refTrans/files/model-rnaseq-pipeline.png | Bin 0 -> 39404 bytes refTrans/m-dataNsoftware_HPC.txt | 146 ++++++++++++++ refTrans/m-quality.txt | 234 ++++++----------------- refTrans/model-orgs.txt | 12 +- 4 files changed, 213 insertions(+), 179 deletions(-) create mode 100644 refTrans/files/model-rnaseq-pipeline.png create mode 100644 refTrans/m-dataNsoftware_HPC.txt diff --git a/refTrans/files/model-rnaseq-pipeline.png b/refTrans/files/model-rnaseq-pipeline.png new file mode 100644 index 0000000000000000000000000000000000000000..8f476c7ea4fec92c5ed5ce833b5daccecffff7dd GIT binary patch literal 39404 zcmeFYWmFx_);|b@;K4lv2=1=IJy>uL&cPjmySuwfaEF7tYmngXlHl$z&2#Vny?18T ze4CFii)2Al-Ceb-YS;eB36=jQ@c{uJ0RjTzgXC9HMF7L_25uBsvK(`4zq!Jn?UyuUOCup^`rG zP(c9{PNMI{G0~`1jg-3}C4xx|#V`^5{S4;eFHq$Nx)nslc2^@T;7=>qFW1tZ2Q}6g z_Z#OIO-?=MTOj10{GfXr)ggV@a@8?U4!ggjV~qCQLgAu9BCA7oZ5kUxOWcV<%{;a9 zrKT2t9jLoGVem0(U8d`lr*N2i^^PUs6d{9wKoFbd;Ps-w*@a3_s_VyNd4(RiN~C59 zVrIe4h0iM|+3u^dyeZnFZ5>WYrr62F*BHlZW6O$Ym0wFE-^`DI&i$T3zG0H@0SjtZ zmP~`#y``F-zRQq-V8Tq*9owhKwG3qu!^ArYy?}y%OX%ZnvquiU=QhoPa49s|34;!m zG$=%yPchp)7&3l1uT#1=V5UvzmE8<+k(syxBSkUc_*B^>i9wQUSf!fA9#;??A8X9* zncjCUt`mAGJ@XB@J*u2j98F{iMPzBLTOO4&wqrxwkZGoWs_RHpfdE^<2<~o{%bgNK zT)c+>uVq{k23hLF^=%jySUAM4foPzWOw152Og*>93;R8gqpSn^7^XBFc={O zc+YNn<8n~K;(29SV4D8Wh42jImT!`T9kStr5Fvd>W-Z_b<$0)u%8W47(Pml^TjdbV z;`XxtK!o@`%z9Hsh(u}c1}WOeo;;5ZFMmD2vKtCL2G%hLfz6kJ35rhd1Az~#?YlK% zhyb6jsPB^AF?@wG6NIi6;$HxkAoBvGlfStw#Eh@oB=k2w4qKS-&}v;DY!Q)sBm@zH zx}}*Q*L=<3{LWCZ@`Orbuz#YY3eTZ_)P!>t5t8S^eAh1$8iN`qct1%|j7#HRnkypD zcSY)m-vX!Y|BwSaNz4ZG4CxtwmI5`c|6u{5!Dq3WhX78ZXLLi<1uFvXl^Wzx=U7&7HJe=SV+`B1pb-0h`qQWuiWRPf|En-Q2LY0W6#UTIu*c|ya)^6t= zOmMD#3NH2l-gdA?#HN9&{+xjn(`<@iimD9G9)~e1tM9#Efx&D|p*m%`fD=3iT6!dX zpQ>K8o>p~o&ASEtBOU_e!7jWl0$Z;p-0v*b?3K`sDDx4mK8L*@HYv_s+Gy4=PyJB3 z9nK6qK6fB^6L=G@M-vTnQpUh*LRuemUAQtLUrx? zZa`fVwII|cNkH)}9#b-i3`dcSIu-1zENM=nO~FUe{yBd5i;;vqu7;qFTqRjHu`W4E zoPc}+6$Pf86oHuV;x}%oWZK_ynVNF#qMlNoj1uvz$(PC7!xiym zu@;Gk$?oKG3G>O$$+`^Ll)9h2rCti3nxGWA1-q4e^Q%>J#W%q{e1Vc;X?ckTAqHs% zp$7h2?EMPlZE?PFjp%YEf`6)ilNWr(7oIL~{oMqX1%FurEirwP;QaYG9sj5Qis+#I z3LZob(%PlkCD_H^ts7B<>0#2lV^zc}h=`A9jflQv+iRZKo>=CnGShi&qfWf1v8M4z zBuO+*giMSsl~J*o#hlG7wJMdHVKS#M*Ee4@f1G6~i<~H$aL$^@)Mb0!c?$kp-4g27 z^k9D7`Wg<)}>x=H@Gx^WEcI0h|Fp~mlxO3uP^(v0eix>X9Nx%>** z71=*KB|5ddN4h47m-NDg(z}au!*V7}xC|TY!p`}(wJ54X zGbQ6BBZp;`7G~Cz^;#yt6G#?&rW9;t?5%@&)hikr<>GTW+`YU*9{H~6_q2$7;p*X; z*sIvObjNf`^oMjr8q^x^)vdsPs~>uo`qC`44Lt{g%(o5e#tT>X>MKg<*PS|#J$cf! zt5>mGzqL->Am5NYQam!l@&&sDwE6RkKA|av+7%2;9God zdZa~>Mv0XOA518U-^ZUo)CozhQj1lqYK*%zI5j^l;{U=Q${)#((c#ol=E-=c{|xo) z^sst4f42WH|0oB|1ziq*fbbJq99k8YACUzf{rwN<-fqb5RzD)f-@0VqXRx7!pG2Gj zdEwcSw6Pu>a2e>>$|;Uqn^au9T%4il4a7u6`dYt65B(mJ2#>^;#&#ATmgpCom5>#4 z7fTkCPFAPWsM`^b(2O81>C7 zjBdAv9Sz|}qEG9nub~({CA%_}SGvv4yU_FTt-+6R`vbK5k>}cXpr=^qN=YoqZ74aR`i&p8x);7XB?&PJ`P8f>QzgI50>F3RX<>5xuxeUuRzAvh#$L93 zCXtpAv)L0Bj6A`qSQ#n@4dd37s27>*uMuSF_8eCK90R5rr|dy{nzH9qZ9g z!3kfTuDp=Vlmx;G%)EUd8EPCu@NY}s##WqG`6TvuVfoJ^uKb!{-Z6}BZuz7lu1 zuHmV6zdCjkaNJ&SQ}F1#;RNTrvihgxQa#l`LuYZaazndfckCkXqOn!&s6$m>%T5o) zEZQg<%u-z2E>&eNFbp`t?^f%X=pd^ZvGUnt1DaiynVHkDZVRlqeYNiT*LmGc(Uw%$dN^6B##9vu!tM z^Sp5_;*L@FN1ipm`3{=vGxzz{;2qIv7838GYnW&0AIj(Kn}$ivFwNY#mATM{pB!fH zGVbete5Nx$jenvSNPqRXinvKls~a(qo`@aC%1GzVcJ{gu9tb=S55)=oqUOlo0^a$Z zSy^4VHtL&PrE}hFU$Wf(_@t6tW>U$m&EuKr*?t`K963kWnw6w`rK>!m_!FPi7>t<-dD->zMm+WpcMy;|COkan25NK zSaqwp>OUKeNS@vf?4*8axJ{fY%l2}9Y&S6&?mM%XI!qOiS}X40ewe$l+BZ7bnRYfi zEi4LKE9(%r@qNyHG`t#&{c|LDCzqHVDd6k@|8(I(VP(9xej9N4$nr(Q3#cjKZNF+b zKtLd3y#0fONXx(l%FS6bWpzh&S(z_}HdYM!Mm7c@1{W(^pfv;puge$UM=Ow{KCz3H zrM1Hs7e11IT6_V1e!I*_Li|q?M+-g@by<005gU6DF*^er0}}~90x>Z$uf37+7e!I= zf8P%L#z$i6=xFngw^>#pAOs;KMTL}IAdfQNXIwQ*1z!iHel#}zwsiK<*c^Y4aSZs%D9jDQH#Roj z3#!pKJ}WZ$ZRhEZ*uisJd3#c|ZZpx^KsOw{cQj?nMwrltnII|%{h#j~5qPAj++)=` zvH!l>1@|>KkvPN$?ti`=*&wOGpH@B1+u)Gspme>Tf4E%Fsy+S^dwsf`*bFA}I^dko zM?`}}L-ln+ixFOCxF}VF{)^gG?>U(xw8}5?fn7PpY00d&w5$?8%lqY?ft;>!#XiaV z@woA{cHP^u<9<(nK|lC9Te|^A+c7qi=e!Ric-{N8eYB?biEG2>T~Bd-6#ib8`$;WZ z$;9!3b?b4@QQgE^zCAuFMbsChD6h#)K~Zt1(F zVv4XXA9AMo?+$^H97|N~FY8{9^ldK>hnsDedEir&zdXWW2+Mr{FqUoFHuSa2X*bl%Mq$- zfD)N{;7GxR!m}h$B$JHrUvrO>S4+eHGC~)onU_{C8mIqg%>X59bI`RM}1~`1@ww`#dSYZ|Fr!U83=I`^(aWmG&SuN7;8u(_?!siX$c}v%sCB{Odk;P`@{!1 z-t)d{AigGb1?$;f4`#ISpKx?s^wq8EXMY<~HGi~j%1vZ{ILHs~_J>UmTx4BwTDH>P z48%bO2CK6GjJlU?T4p=A#p?B@u`?CxPWSa8nOc2XmQ#DdtSGbRcHP_NdK{Q_C+)JT zjuz_qMJW7w|F+lr3GayY2-b1jr0@2_wm^3R%@r^%na1X{4<8GXb<91PMyT=wN{X_U zDtDJ{I%VtT4KIf&vjZYI*A^a9U)px!q@h%wR-M=R_+As$bl_W$8?;gkqqut#l)jlg zJD>OC6b;@(DUG-OEiEf2qDV6(iLWWja-;d>40P@N3W)#tYT8kkRO8S}_p`hD-Jemr z0lYzD5@0=cfOR~9{c)fNbskho21O#z0yPz#csDt$X;Gn44+#tnER}{c3$ULMo&ae)X#goyNSNk1DcYsw-?)S6T$=7Gqbgdp6 z7YxD!coW!{O!LVVm}M_lUA7}ty{^t>qxc^VjIup%B=cAj>6$kDB$dtS*Ku?`Vnqd3 zje6iI!Mu;Bot@8wHf`tV<7G)}UXRCQno$)myUDuzp_Q;n$D@p$GI@RpG8~IW!$7n$ zL^&e@4uh7l)V}Mxz5pR#IxX%ZBaLKJtUYNtR4wbdiR3wH8t}7HUin>XoMu*N6u~_3 zwPDTekc{Sh*tYjW@+7jH>(cr8Akp%VD88Eo`1!V3RV|HI(Qvoq!|~q9sz1PdoKZL~ z6S^POIIT`8W{ZN@xOZYi$~*HqtuKJ?sCzLM#r@jP_O^b0w_IFJ@p|k1nk>407{DJk zQkMA$#Dm<{?^y`;Q-Wy0tP$LeLhQF?k*wk)3I$i07>x%sn6sxpMcH)gz<<_Op8J{3 zN6v&AhqAh#59@w$XrU)Q?#w93$!>KqhWb~MA@9W{>w0lp)pY>*(2)X?4qA$Y7>b;9 z1HlTsNB(%yYV)uN?CH1z3G*;Ra8AJp{I3ZbFcs{SW;$KNg2A6e6(v)f&gUBj1U>Ku zIcXG^0s1Q5@5?(s5?mak0x7shdAIxm=1{Np1_90Qvidp4S=PjRU9h|1r6+znAt^vFq+l|+IUnhl9Gt7Gl%%g|`s8;Zvu*XJaV4Um-69o#tzf_I}B23_0C+>_><=SETi(80jVs{8;!?b z*`~9NkVs*jlqoF5_Z$8w9%|_H7jscu;d+lWSWDuy1Cs53$`0JEuLO%Tl~htB(LINdIDf1!_m-wjIyHhU@trb11CdZ zadl9!i|VMsV+!)XX%e-rqO!3FcrK{*q-C`5kCn+2Yd@w^Krw1i@t8(xKw4Osgq}HW z_G@C4xhCuJ z_ocni2ByqtT9_kVvGC&hUhny(+5;rQEZa=tT44g*JBDzUZg#KBUz`n$ay+)%IdXY{ zO!UDa=+cbsjt@MHy;lMFx+%g%R-<%vGRXu@I--ZQ^kid%#YVjbBG{N)_M7lHI35jr zVrA}K^q*;-=H?X?IOm~uNND!5el;{D5|?Fh1uN4&X_9JO+Pf{OHQ|cVwl@=qhIi0c z61W}Y%5uaWYEqB;1t@*vCVv*Su&yXyVg$|_aVSO8WjI}xa6wWiJ9D0mvI8P)+fCR| zM$?SH>|(X5cv{xY1@{PBuxZzFS5J1OD^DN=<#Pa;aXwa0(rIeE+sx%X27=OHdK4=` z$o#7Kq@qA@+O&;4x-~Kk?8De35mzd-z9=zBvfQh@A(cY6j*3LrXU-`?`EnYJmyDl4 zN*3S?SU1l~P-z>A_KP6F*Jf|YB$UN%s8ZkwArR_R(K3Z~mUG2(d(qTd5-Bwj!9>Ko zqQkjRsm+Z(Yp#&d;OKmLun0-g`B*(U;e9`<3-wYr%LRgkO9(}mBr%YvpZPlXlws+< zFkf7STh82-?l8{%F84SF14i`@(!?{o-D925btm?-X67`?>$|rtx=Wfv1-pBQQU4YC zFM8KV5XoAmjj@cDYvFs9=UL%Q1tIJUgGc6Hj`3;yAaZJQr}gDh{JqnbG7joHeh41Pug$StQ9areuaaY(#V^HV%xuD z6n&DZo36A}hI?SJI369XRvYCgQlioB)cVN7$ST>ZXk9z69OUv z@|k2h zzbBYp_#jHtJNva@z254PYK$t0*Ymc;9tYy+R6`Olz<^_&NT-`*s6v6Wh|4$X`$(a9 zxS#;#h61vA;IV-i_ht~Bx0n%>BFBgjd6p@>OkVwrYN-PqjjtrQn4(^FBtG~|JRs5F znlCqFOJrR&6r+J=?Z%RvZ!?WoOjkooV}kh=b&p!piZ9#>K6wWe6RmP%2dUJu)G599 z`cjO&>mzbrUxH2TW8Xwr8C$JRXx6QEG2i;FDmw8pfPA#@fIn7r)XYf@eU;(pp;iF{ zD9fPW=&@lTs6PsOqdgPT*qYtD;_fvU9DaRz5mlv`y;^Qq^p27`-suh_k_ocg%hd~w98 z<2tjXMG1oCxuHn@T+6u4)jSU;kQ#A(bRrzm%PX<6c4q6BuX7m&S2*Cotz;jY}FJQYkI`7*F| ziBIqY@U|mC{9J;Nea`sFQQ{1U^wyQ77(kg1`Ig$yq2I=Z>YD&jmDrCEzKg!D;o)5c zWKka*hVb`^Z!j+>0@^5--mt1zmV~Jb?z_uoV1qhM2fUS!wktQ|?Jv%C@R=%IHE9G6 zqo%q|`?at~od);Ea#dh)A-Mp3f^CC?`)?S?0U3Dw_vl6G@IQD_uL~nam}f&w5aPdJ zW!)R{`TuOUfdaK1$*HqonjeWUXxIH7NBv0&8$WhS8t#bmTq-@$x@i+eYb?RDRi*B) zVzVdZBtq(eJFSWz0OwNSsVvn}-6gCre6i8hfjX(jvhLU16+XOi=Lk(%0Gc>%uGtSX ziNr%Z;}rcg)nEjMC_!4Iyr#KQG9XYa=<*5vEH`0GWKho`5ZFsU=C$Y7G}#x8;@O`^ zfU>udX`E|20DwrWPuEr6R&Ji9|ANIi&(J&elP~=8h@v?q(BlHH9`5x=Kp9Ffr|a$M zzDr{#wwNtTl}VXl>vufjj-2Y6KO;q3fk~0+_?zjp?K0>cG8geo`rpn!h3R8%4?xkW zsSSN0@OV_$tD>fj<Nz=UK3@AGn?Ore6O58>7m&@!5hBbzpNfDI$O-J_$ zzS{sIuedTTv)hD3&>(BkON8lMBIM-bA>o}8KOXR#_7k{Uz^q6^@Q4B$0rzwyMkPh+r(fnTk81ETZd`FJR z(4>c+s7ZU6{rXfidp+>cX2Irf(TVe#>!6pCjWSexaSl0v#Z1w$VvLSz2EPs98b<6C zWqZ}~0+0$@)AVy|RIv+XQHH$~Iu_2OCqRN|5x&;#CTUSju~-P0z55|1#u|EpIczuh z$;9gCQz&H?S7=ztZn{AjbsSJQ^Wd2m=f(k@;oEHaI!NeVp4mA}kxrKU-Coju=yNbx zek4cQZYV*`ZrkD^kWZ;;bBJlEaVEeH}k;HmLk>NUKxFKjI^xRLC0>Hzv{%^=h* z)jPhQ?|cD0UVit#{XxL>4}KHY?anw1dl$*DJDm*x+RWp!gizjR&q*Ozjb>owC(Y8DYp>H?!8F zXn_|lJfH}~P^0pEt!P{`ca`HgW?ROhk1bT@2gjRHwS0wDa59nc>{%|njNW5Q%2%g*@!G*n8RtfcW`SY5ECCVfaHd|AT1@qo#2M?|QfL3y?%Gk~EC* z)KSE9CtqP?wP~Xn75zt*oYgGF?S+t?E8p zH7~1>fw}usu{jWPv2sjIll@sU=~v^@#IW_{6(=ID-LDXT44y>#=N5q02sG7h@U-3e z(H_pwV%p&1i37wi82abrP4Bl^@=!r?$b>kq{x#68vXU&Q^L(^C2|sgtWDh_sXj=(C z9)!*x9uBKyhbvUA>pe9<@!M49WC?WWQnx^l4WvvPDm`yk=cDrF6sb@B?Py&j1? zbZIf>v)gaVg{5>b#T6a1t1DhZ=fp2JJ%J==r3QuyN;Sz^L@Cnp&^s^rHELGXeuLnC z+V1QEY}FoM5()FJd2>S6fh?f97eX3^D?hoH;V@qPe6yTh#Snf3kX?<|{OxNYZSy|% zL+njN^ip4B@)5(Z5^wR_)})d?T2>2%R6=gPJf1ev0Zffuh&5y3&qwY>tF@OxB#%N# zs+v2@-Vm6!2Se}{S7_5xwJrh zUo`n0PoyanZDo@|BZ(V23nHRX`wn8QIlWJq(zeqV1o|+2&;#&FTIbXsZop#r)3EL= z>iq6;1d0+nzxwG7vt|3-iK$8HN6CFbq{jydAGx_pr0eYE!zDP8;3 zZ{t95Y7-%X@W0i<`CaRrs2wQDdv%Q+PfXGW`i#orf9C*RxN)D_r5_0lvh*j0pKZWHn2zh zlRqp9hcx@hCct@{&rfs!*`98>@5=ox1S`@9a!Vcn+%BBch}vz3xlol>Uzyjmpm#8A zI_p7{Evaza9m+~B>Ku-Chh&8pozpY!s@SrNH3$ixk~x?7^ETVK5YP`*?$wnt(4r7f z9aIss=L*zPd7wFxPom_cp<#juB};kLz@&qIV}BESPjwG)kmwENwiA_Sw}Haf*nu^* z0_k}Z#>wvm6jJL$V9o=PBs2at!og4VB=?oR2}K{VEA;pO7SeW)Dr- z9UxghC19zle-6bqF_gUuvK(>vkm#Aow~k|{;01;@K}qH9nIzTS`D4TABW6d(6t2JG zp;r5gV_7|0w4ZCW2#~g*p3{Bu2(C_FW@kfvsKEhq!jPie!ghk2qM@SH7s@4*d)$-{ zBWb_IQeRR5ink>|wp7`kMh&cFLmHUGF;L>9)qQfEGX_wam1fEi!h)dQ7`rH1@Ofef z{77@QbO=L&vRpl;CZ?irjGYcRV=2i^-m)&#LR8`MBLAZVx?ydh#)Yl0Z$$W7IZ}f1 z2+4?W#8i3`hBm59sD$q&Z$RgWxb0Wg*IVf|3u1T2#cH>@k=*3ueSDh7rHN~+tLgd+ z+Iy7|==+b?jmeUCEb3y($m@xHO}z~^K>a0m`zIVIhhHgjI0O56a^v8xOrT`62XSjS z$O$Wo8w@hBo@Hj|r1D@jV5?IvhuSWp6lfmxbrfp%;Rk`qJ1!zUgnEh zNYD<1iOVvaxaIdHY{pBJzAqyKN2ZlU2o3fdu0b=Ufuu*UpwYo%fC6s&_xGeF)arg; zns)vOkorxjec!~F52o3TTWcBdjYq0`q%Phb@yYEmxrx`M4F3W}bLfB$V&PiG@z1aS zJQ~&44*rg__1|{`j}%_kTq~9R@4NK? zn1l#Zy7eFQBM7Z2Er?+C|9XFIDLtU8RqR$Q{rBAn6vTlI|J%ubpx6IK z@;?jkf7j%HV)4Hl=RX-E=l{191Vii#paxzK08qR;8~HI@5BGGlY~vdp2OKrb&`YJ0 zZv?x3ihfYd#G7n{RvR$>Vp6-;ZgfW?d8p1&I~q) zv&#vAWfK7UApnp||G^zl2@eO30R0PTT834&@foyk_aFMr;d}fD zfH@dcYHf2?_!WTt+5|-K`iXR#4m{itE1;nJRjvryC+>ondwef42nhX_w@?DtdEpE! zFxIWAdC$WlKC0hkCyGS2Zn1E9Rk7%-15ALrE#BdY<8 ze5r#?X}ANdl55rHhc^61b^&}3L){qQ-)PtbWcClh>abJ?B+6d>;7{0v0IlLK$a^`) zs-GK(jqt{OaP6fSmLIg;4|yC`GzJ0g4rV5Fz-6J*1K16WKulr$e~a^b04QmUjV&$# zET9_dxw&J46 zK*Lw?Bj6}S2BuXqUm1Yn?lOcwyX=JgLn226kilm&?9{L^W^}=|1bFTq5<84 z0A$N*sQjiS2fazeN5F#hh1A?V;I(NUQS*WVsw~ik(7+ga=wi{E?}PxbQ)~_rd4oMy z@|$dA*``DbkMbg1S~i0|9%c1CM+=<;%G7o+n4Z~J(=+U%8o&$-o<@OKs?{c`x-hrI zC-5tVs=A)nwp*c;FmzTwPdgqN^q57&v2|QF&wx$zV7A(zZsDr>W*u0KGLv|OW1{^^ z#D6hFx(P`9jX`JMg&J|47W8E$h8nB^-WBPMhJCYz)V3l%M}atOdY0ONdKU)MGD{yM zE0S2R4?YAoRHOF=ed7wto8hC7mjvm_pS2X22|K9#jMvA`S46lj?)?nIo|>nN(O_B5 zRW@gSLf#co1^@^$0edIa-oGAK4HEOnq08lve>q>Hb zAf-?>9oGU`Ql?)`fDs}Hh$YW=@8@eLQca^+F**Tb5A+}7%!33s)nD#+;*a%{RbuGl zU)VCc06Pahzyv2<8Q1-mqD@8jV?WLfCPa8M9W?=cdgnPN3x6-n{h3*P_@VpxtQW=1 z-PMP=zy%0WFp#vC=`L4ts(-U`8CEy!XSo;NUl(LRr|k6_Y=s6ZX0ZVaFs*h3lj+^r zo5_Z*x(iA_UY4U~T0y5-E@y}7(=c%gbthm#Svt{Pv2oe3>l8v8=qXdShPr-8*oT{d zg<-oQPjb&vB0FK^wP9#G6IN2M>9GP7At!*j=3sU0I`Eq^EGbAr062s+JmH}(zF-rQ zKY&@Kb{`9&PSF((aQVE}2qE~-1LSr$lH=X1gFu&0P%KiikEMaUe>nOrWCT5*B$8S8b2j6-~b`eQU~GQ5%4DEX>In zVlw}%zbMu`U=WdG@W)USkM`t^<`^J!x3mZcMR}j)ifXg~J?RBtH4f8j{y`q}wb9Pb zuGz9XlzOL&iyq|X&6l1JFo8({ypvTPDRIrWYwpwOC@6&h*Sd6Z4vT{?$7?Dd)MQqZ zeQC!hgeM-=RKtv*M#BmDwEl-y4j`AaKG;7wx6?}@yaSzl7=lVl%XD5_g!y$eXw|J9 zBX8R5e2xlCuKSx!EQ^L${WM0pz#gM`zKGXE-eC22*0qEYMX&~M`MJkRCn2JM6f6Q= zSRZh1;GAuWSEihQK7SgeIH69vA4WHuaDJeBq%%p1dW3*Kr2F!I&`!GNuK&6F1Nly> z*TbJ*AL1Vavj<;zL*2jLH1}}>oTzNF(v%!O!0$6W#xSaxrfToS%Q7S{dV(Hw+yb!? zjqCda=tt<_&2#X!24eqBQR12Z@g5i$4Cqj9D$#BPwp}nZlPP# zXJ)9uv~8uZq?(cJH*46N8>ev0KIO7ry&Rz!`zVi59lQ-s-usy4sU*vuNK_j<@(zVJ zh6b#@(Y)!h1sQ-8bD%_bqQ0ZP8CAQ*i8xVZbq{Z=w6DJd9|&ckl@&lT))j(MCp~ub z=Emd5?|X+r3-7NjC&PCcqn{HFn#1t_TyXIg!LdCirt_=h&L?ses_KptKs)fFvW+mo zMMQ}WJ-UPky?;M6{d~Kgz1E3q$?RjSTKrph|GwVp17y}|2V8_N-66-jQwI7(cPT%Y1wr z>unkUveLFxa|(#Jq8KsCN}wbfVd8|6yKBgs%v*^0Zg;ai+n}x{Z(NPkc`-t(53tFV zN%X71RmVW?&G@F~aA`18IO7JD_$>FTUvs;YL-Qh| zhj58w8X#Y^d}xXhaTYNZ`{<@A@ds-uOGn=3 zVtf0IMoe9X+BZ}+VVlO-^$+B~8~W_y_(hd>1@1-_X>}H;;MiW69YYrf%K%j_RzJ|Di;3xH}zwHEA0&Oem>MKk+ID$P6#yz znd{Nu+&M=H%l_Ei8RH{_Gf7ndXZr?hWC0R) zS9B`<6eX{|8M!Jw3dTcC?DgsdX&*cqr)~PCc)d`@R(B2IIyap&6(^l3FR+cDF|8ia zAW14S06Sj;Q^Hq@yAeas=e!^->+>Id1%W~+EL(MKce6EjNd3v((G%F$q`c1;indtA zvqeBKJLsWyUV6H?c}ywXkS8Fc^oIwnAZ|Ib$XAG(l5J1V+)=}mEy(Ho#>kirLqM`6 zywFT(dVM3v#-^X_{eNrzlECW4x0cgrj3=ptOHPn##LRNPb5oNBN4uk0e2p1oDHJ*< zIz(YLz|)%F;$Su(NPgd+^ifg{1Dn%>adkb4I`60$CVOPW{=6IXUSCd)`37t?qMHqY z_*9*7E|-3c_*7>nl$)LMVz@CieFe}41?u(JjE{roel5l;8bts;?^?84II6QvnGaHfUdXaV?- z#6B%89A=ZyoSvSKF!Wb3JF%Yqk42xqeLI#C|MP>^hJ!s5%~d`sybT1;SrtEEcqe>}}0f{1VlG$Y-Dh_e2= znCFsjdile*9Hd(m11a?IgoOc*Q^_9hb&T7k)pwQC7~v}CQ`VTSu&Mpooh&NLYNRsyR`NKj4FElrWK+x+^RoS3wf<}jR!1LHLh|2W6PvF6sXBN!+wOiq1t3O zQLO=MSd?K1!#p_frzoG^QmQBzb4&RW`^#3?+bMISDs;mxMJN2L?Ci(Zb5Z8=7{U+L zjBP(rHcWm|rr{uE!^AJHGy;q7_4=^Ludw;2+WLY+uw5U^6<*bN#=rQC+?=JK1>7B8 zWKt1`O6d90gk#iwCWC*OfMV~N{5;(jm6!~1#qlF=Wwu3CM|<3_!!1-zi*p-me}(g) z)e@IpXEW-%TS}vwr(Alh*_!eYKH3f{Cyp=-x-VsV$XJ)3=#CzMVm`zVJSuX0K|%ZC z5-Pupp@6)Undd`!Ucn_fb8*N)-cd>)sxSL@kNQE+qre(nLZ3Z6(&p^zmVis+w(GJR z042%KKNzl#&m#kBum|*r?Q2N5TU?VpUnmNvg(J9v0=EDz0G8$M(r1!~t zn;E#g>t2Jf66jb^Bkd5n!juDJ^%cClxuzD)~pJFOAtkD(~k@aAwyP= zyWjmQ%kAV82Lv(WPU?uP zqo!hR)6pC0Iz_aA<5_zq8CPuhfngx3o+P&7AULy+Q>9R4=Fpv`-!@+oNh7VVn7|2{ z<-0e|>cq+Yr;(rRK~F<-;3bHqw)UrZ7FZM##UqPi(O@i<(! z$anqJH-tr>!*(u2YuLbyKum{_E9+8EpKoNq5tmYE_!%W19!e48kR7VQ@M)88X;xgh zPgM(W4gBg6F2>A$Sq9FTX?gCpb9h>Kn+tvK0@p0T{ zZVPyi2(NYW{%E_S$1~tDhqA}(7&xwOF5l2Q0Jci|{dNc?n$qsUHboQlr zjg9F;T&3u=6)A@^>xx%??;lWSRcl@vUFN730u?IJ7LCy$wJh`^G_3eo1EZ1%Kg4*6~E#=Q`1Atb;I1!_2^V+g^wyiNikA zM116%>F;bsDlp&(i|F^dw%Adzr7D$KzVw6A(Tq`%8g>RO56d$suQJd@>p-!3>}Y+J zvZ`R4O4AI;My+ktK6lx(Iv{DLBcK@O<^}xqg_x|KU$?sb^YH2GZ#KH<{Su4Ym}$yV zw)|@Wzs$miuw-&04BxBmkeQcMSFIVqTW@?!f>M^eUCedXrSv8$a{D^13p5Gj6r&Qf zvuyHrVKPknP36TC?s2}ezFLx=HMn^R%VuFU9L<^+Gq zn#%^_JvCfK1>XGI5wUXZhQ*~a^)D!R2^SZH9;eJhxzbG^rAS@-HX3^nXw|dz6Ut=j zAF}32w=Wd$e`7zOl4Iu+!+%0_t&jny;wZWrFhyI>%*}kete$VbS+bx!hdB|52avZw zcSO|v1AtRm>fD|NC2Obe&V^z{U#Fq#c!rCVx6*&a9W2KtVXarkiqhIVv@YR{9Zfm+ zEp^?xh)5OPGPSer;?TECJ{8HzMI%_D6Ol*~NRm(OlZC3T8r|;7)UAy zOgrQYr;9?$cn>48hh@BAN5SbQg52#I(qLe&+26*22MUTAO&TI-|~Q z5rNR8%q#V)MwvYCjAy50f82PsP!y!)MgYK3D&&5>yyT&2y7`#wU!0I~2su6c)9scs z%6a6uX!ZH~AA1NMMGreE^hD%zzrkIIowagXpO_5j@HH|_}Fs=}~4D{oDb?`qk z(NIo{2YkMROF?Cj$6fTk0gpl+{)wLZDE1#KJR}LK$*25QcAE(+k|dG%lJ|}=o25Ty z+(V_C5=g(w>G!em!Z$_$qhJiR0Y1@kVlRLA8uRDD)QGjc_gjcJ_IKNEGUTt%CbJc- z1ftm*hgHQHg_U$(!tEUr`$KE-lO_i|ucEXberx^hA8rxkyfUo7GY=ARVdA$I~-t{n0hdJ--d9ygE<16KE=YH0}Jnd5^ zbcNRN$E18d%ca7G<44#EX!&Se<9OmnUVe=FRWJxUS+tV2wRn7@6A^^eWq#7l>OleHo`KAZNIQ9Ol!aLrLyF$A?mZ} z;Or#;1?!=_lQ~%oV3C$^BdL1r-4^h>C-kl0_|-LoWQ_?&MKsb)G>yH{XWTwxgcs8 z0*#B;E~1_m*U=@A@fegDF0i~nvRuykd6Q)OYl6WU27FfRmq+%Nx#49Z`JN7k0z(JQ-S0d{(&`g7*^hS; zg9=%(&Q8p*AiCNi5(-mrD9Z8fY*1G93(3rzC!ddU^5hYgen23`J$(Hxbu+XQq7X^X zs@tu*-T(Cm;rB&(Ilj)V2a_C`YM+Ql)VXx=+s1Frui(;bg`qQ0MPjwrHfN-85D{HR z`53O-)_vDy{zdO4)i-ijI%AD(I8}w6cVcWklz55xiOFHu|0PI1c#~*O%xgd6YxeBfAC192(xya#yhB9E76qJw#t$aqC0>IHf2uzvNxV zmV>x@K!)o6?AL{lG3ex;k^l+Tj_sVsPeu8*Tk*$Nu zf3qU6b!E`Q7sbT+bM0<&z}G?3+GFEqN9f(8IA}NuO)AWmM&#B{<0hj>rS#+JV|Kiop-|YO^2Z%YO98=-}kuydS88epU7ODG40*0ft+P z0V^BcBXz@`riMsUWSTqec-jqj3aSl(>vyK&ppwxVu2EC3M7yrUBs81|jf2YK^2NG0 zck4+l1@wJBa#kyE<7zS+LxN16qUG*p@<*&y{F?|=9#cfkMMh6f`*q=AxP=zh?NJQ}8^ ztH;{`1HuC&9I{0-Wqp1&L?jU_i6pJyDz&=K>rp#CG0Wo68%mmG9{l$tTs4&}s^uvD z7jBVb9BvO;tMocxwssTDp}R#%8}t9tLR1hgM#D(@t#yf{V&u8vPW$*}PVzRJ<=R2_ zkyg_g0~=eSOuI49oR4;Bi1{Iw9^oe~_$a*2YjuRi99BH{z%K?hxm|P1EpOdwHs!ca zdGFOP(6)oBJbtHg>AE(0p{);fywJ3!{;vH<#9OGlT;@uJo$uX!ke#Yf1#hE9=Tt!+ zqvNVOM)+z|V{9YKl<^u`3_GFfo>~8?kTSZ2g_T@DYB$+w4fgU{4mP__H4eb?XK_0= z7E7s+LOn-y3ATG*)sylZzvO$p=ti@)kAfLDw``q5%l=y7KY4SHT+>%Rab{}@`~OH= znI(aD5I$02a4Add>^vE~rs(kL-9l`<#1CE(D0JfJI^5fZc9- zS@tjd{O zhc0ELkrF{tLK;z;kG?l5~DDje5JUejX` zJ^e=leMEp-65}~0t+mDR-_)xN`kX94J;ci@HU68MCJgHa0%~sD5^NAZBLV$?aV;1C z_1G<_7W*$cFN2QF2$nh?rE@6%lPIA@@+c7>w@Clzqu>kJ5<0b_GW1FdQsHfd8aU!> zgdSJ-%;I)?C}b8K{nTsFTkj^ixR;kjyYi?38+e#WT9e*b2^W$LTr zX1nT#gQa%Y=kxdEeq?yAbOf0Q-}~49%_k=m`2guRkaZ zs)g+RSgm=A(2J}*L-A^uj^oDT@E#~K z3q7m&4aTZB;tLKOPPH!U;)%A1N#{ z@up_J92zb#gy^UZ9=erE7oD5kDG}FRjK*h41gFy=r1RO7$-3vng?I$Sm<;2?r;^S0#F3Ve+S9x9-!J_cdFA)$bdBYH z;;6>8Q7$`$`ka9D6TZNgM-i}hqqb8ssl*gO2HPXP4<|7XKHqo*kn~ci#s`N5Y*EnR zOnAB;rvDh%pc;WCyEAAG%>Oh8IP5L?c+b_!Nwe}k)c^nhMuXK8vv=Y#om9?ML?CC& z(@@{G9dpV%&ybp#dE5krlgc(9+bNqHkHXWO)8XEp`%Ue;shQIP`|DHz9K1hBjSV~{ z>QPjaO&p}pmLrhtgg1Xa%FN4)Z2(CK;)tuVNPP9JxcHXhjV+n6fx_k7i$OIDEa+z7pc>?&h{Tuy@;lO;Lz9}BcjI@w zhJ(uAZyAm%+9?<-y7NJ29_!BJ=Ph(Chu>-j4)HqnpM(eF(|*yJ465D!>en{|W6-mA zmr9bUnK$2(U?1b)i1DX|K8pRY7cSc2wR`M7JYiM*+d*m;r$Q>o#>hc7#w@R$YsPv> zm89<1m!L6_FQ1u1e!m~PxNAp7X7OpSzj0Rrm4}e6VX;M#tx=*m{ZA zWi>5#CcH7ie{-50@_V7(S!wvsklM2-@{L}YU%R4;_LeE?QT{2jiGh94d?sfa{3yy? zboB3CJ~YQMsAd%KK4&)-dD9QgRlpHCx`i9v^>M01ahgU}c#$u^g)HVTX14Ec9|wq< z2kTqyInRo9NTUBz+aGngCI9x;1tK+o2vSV1;WyR{y;)1UUgE9L{6xbmZKAYYqbEcz z8`^$!^X0@%+fMg&CNG3N+$msqfRpzb+1$|kHcP#|nfE@v>&;qn8MYU$c>tZoCr|% z)t1O-v3Vz26rr&CIr+}Z@fOAA?qTpIyUMDQf4v=cz#n8{JIhfii4NxtGMSD0i0x`; zf4ygl+Js)V{zaJCB4cLMjJ$q+NKe~%% z$|nEWcGoAg91TMM`n&NbmnGq1{ES&D)WvA<7npgoTef!ekP9BlK9FNpHWkX5>@p#t z_WwiDb0E`uFMz({9vTi`T9r!((@7%46wyQmq;E1*b2ac zZh^?G+th1>O!&KHMDq3G&NideCUTksqv8EZC*Veb3uIoVKb%g%+HWnGK&D5EBbu=2 zNh0pyZ6PSFV^*gs3@n0#Oo1=PbaGYI)C(i$BqZ~2BK>X7xNSyX;(X7?A?e>?o3~&X zQhFB%lzm8(8cZ`bQ@844;U~(o7{{h6oL*^Knr8jU$Fn1S3br> zX6IM-Lv90F!PwjlY*mPFuJ=Qt&3u0d!NG)X)lci&xC^)?Xz&tHK~}W_eUC6w?rK#O z%fgSdH)9dvdZ{UOzTWv9OiuVf{8hUv4uFgv5S+N>lW0otxqcm2<8>WEAfYk`jLj5% zbBNtM51>#!2eQixDlA#O={ktN zdAof;nPCD;0n6{6!`x3D`w*||V*_G_sm?Vek0;%SxmzMZ2*!jo_(7p8yj5vh@5Q&j zkvN{z4MQ0PAiiV$ldNjX8ICm&?h!awf3OEcolrMWZRSpAu{U#Dl#DgsZQ}qrRB-eE z@M(%v%ffFM`!ZnRgOccFuZ z(3IZ$?=DSK`c?oH!W9@1Ok;=T>|=5ATn-}X$hTAY+niXH7Sm2K@#qs%;DAjj0Q~n? zE(+jCK3>&%DdJY3SVJTCzAYFKRs{-CCy|ryHa$Rp>WYR}b@dJe0`pwg1D!f*++-Tk z*WUBeNO>(1V_F}4d38F96!@v~MsO~uC1cF3mA}n7<69bs0*%+eJu9S^?aO$4d_<{} zU4aB+4Q|;RkIsBJu45U$pq%T!N4zi%hR}^P_G=^iEM_FMikJO*13&)$#9-bV{u>Z^ z?Xc>O<_*M;18P-ij)&a;yhz~XqrRQkC?_^p6LI|a+AlBa2U(^yzAAcXbG*7ku3c4{ zNW6!hr==B5Y`Iat=rXY-M~9(b1JgqYJ>Op>C52XE#$^8?(C_JYjjq$U7@#rvvm{RaH})^ zu8t}}XdZaj_fb_&VB26eSN2W|3VG~?nO_vMODr=KlP=`S;{tRL-!J*ih`7Ocx2jDf zzSe^@cCp_cRP~>Wy+qs8#pl5`!g*m_)=sL?=|i}H-VQEGfq%ZX)wGzCx`~EEz%YuoF&?VTfTw4 z3p$KTyL^J`*+Ejg4WE>}BFhtkvyAnpkya3piuz-}uiY-OqU*1>d}lDhjqp~ThtnmC z$L6pPXef;-?!FH`;Nw)|urY+z*E-sZzp-CK>AH5vc&3kKlfd`HvQD8vXIB1DaQmZp zlPSaf+V9vFhFp4Q1gpba(B>yIP2hf=hp^Q#y@wc;J$w4W6kvc~@`;=&>BxzYY9gP? z0-C+dM!sxAy-2dG=QO8J13k^Y4gH)rwZ>YeMh3|yg-1Sl3VtGpRp)r?Heg4UJ>^|i z_^Lvz^K;@lQUN@O88jpZLxO#p&m14DM|c6Cmei#8)$I$$^MYI?yW)3!JfoYXisQ1D z3FT}w<>Wl3m_=CGb|zX)Ju_%pa<%=|SV1OLa_V68;dgyE^jY5-$4caqMOt@#WBYT6 zUx8`uukbfogQIs!Dsz$35IYaH~w7 zNQ8DmcYCHZHDNUMniDm+mF-F3^W%U|k^J$a+!IS&5{TZ0jB-a)Dtkbq$dICwqmsEQ zATAoOBktgQOeVNTe3mnwlt_%02_7vMuL@G`+uDNOFC6l0PX+0cG*c3TRK3gfSpDjH zhH{ExTfpiF(As82N5hxwNl|ggj$d}lPZ${T!9s`C$Hs+2G48EjCaIFU3V0q}mDVCw zN>w>JS3%kUt%GS_(>kJUeRm{@qMa#LS{L5z|fD{Hg!366+0Ls ze@=<;c;-wTR75op`-4rAbzZcH?ZT3HkN-A}RQt z^cvdy#pBoOBo}S0>=)KX7Jb>3X1W3!%++^v>WPAeU3~RCNo=u=Vs9-~8iz_} zosXX#nW(*73i7p~ReGxuIKcF+TvBqhW2|(Ngbou6K~%b%D&nuPD(yvFYM1(J?Wa+i zv1e>mFHz$lXT6HFUxiPCR2-)=w~6Ejxu1?9Iqy|hdxT?gP*X*Y*ZNNV=gSe}+c>XT z8Ag#tLM;qM%sjE}cFoe&sv1zBpJX!TuaXvhvKkV8Yos}F_w=3)>?WRowGBf{Kd|8Y zt)2Mg!u(?Pu1tF7qC3LN&)$70Nnd2?na^-g-0f!GjiwT#_7uYJo}z(-_-BEP%|P*S znbiAvntTj0Nqoft&!n=~aAoBn8L@#2ntl(CgbQ&B?)!msa?+La4{kB)3ljld&Giy) za^ou#N?Ak~Rd*vEqmq8`G>)qdFC5+D<`wk|5uE5Q2Ph%O6!Amt1KqB{w<-1g7)S8q zt)yIifh(*VCPUv^;fA3DFs0!kZ2`#Q<)tQMj=Oyn6j^ImGis@|tr}7&e_quySE;(Q z)kmGeg&~6onP#&*(uGL^x z`EXG~OY%&u9$HV1l1@>(!L1&ezMUeNI$NKeae@p^ymx}s#@XIC@Ow~2_+Tk50E2V5 z?j_MBQA?zsKno>$;0sB}|o%wS($PxlG5pKll*Jk_0mQF>$+a)cc ziZk&IQ&8g^;I%c94u+p%DHm+&?6y0f<1vS~zz;^@KyHy;PP5Uk`|09C7vW7LJ9Rqo z+##4{F3ZDY<|NjYIbn@D7g}lpy}YJqQrVafAMzIUP2-bbyv}$P*}c3bQNIb@5T`L^ zx0&LZsRMyHMqqDi@PlbZ=2ZcGjpvV;yU|eY-k(9!NB@(g*g=wz6NaMd0BdJBq1#O) z7j?G(RQ|X+S{6GtWP!W+dpD|iLrlVNN*(SG?-H4;wJsul10RIPM$7{*VKy0$VNRZ8 zN*gT0Yz7kVFi+P!IAF`f*N!Fdqmz205_2`80(El9Ws_G)R6{a3v2_6yIIdLeG*cG@ zb}ZC!(ZOHdJ8gRFR*zw~au(yR` z&?IVWXrnstm zs!EBM6xkH^Ojs1(He=K+E_A1em|(d`XJa6x>t3v)DOLuLU+}gB)CIL$fikX|Nv?Sa3~3q=ACqNlQ6@y6=NC z^;42xlzklo>-F6-=<)0m?YREGIA4F8d9e73ki}4Y`xWri~a&;rrVG*}!{wRl> z^FLIA%x6#XuJ+I&zXy*SId}zdDO~Rvb40+GfBp+_Ny5txKtQ~2=wG1JfR2-bGXm9 zaxmfA=s-Werys<6gadUoUmB33B3~lnBzT;p&%7GJz-d_J_l!PhpeK)@=EkPb_)~k) zr!(`Ow}eKx))%Q;gUlP&+@C8F>PNflU6!}}`?s2~7JJidJsyyIkeZNQ)eovuu9rMe zl71DB8!tjhYIFX4EV2z7jn3C@^+AK|5k#py<{34k356S}+3Ks+o{`Wv%1}V*>osXf zG=1Gom#Ue%k+>9BG)I`KRsO;4bq}#ixo^+XuyHFwpLttWmcIYq>#Lx3nMJ($G9+IheJjz$Ny%sdpF%?hljVgLs zrLQct_Xj9)^oea8S=S#<;gC??CuH?VPVbe~C7#JO48GBur~8pL?jGVXBp$nso_hG1 zdbC^lf%ll~Z!f-&ub4f4K8zhI$e_j(cl-c5)3Jxd-h0LLgt%_3C^9_Z^Xv7I+gBt} zPn!-Bj_gde-uA`_s3df!^$N+6$MN#uY8C}#9b=7_wbF~kZb#wrSxqd{;R&3GeGyDKF}$KX6gE@?g0=UT{G zuFsa|WLGv?sqJ_pW%~Wiv5jz~5m8Yop*FgjWCJ6Y(`y2k>?XvcWc8LOcUw$ar517) zGNPIokGmCChk9{I<6UE4%Q2CVvrW$>W#>j$A{PpFB^IG7S|M}WDIo`IoBm%#H)+~J z2(e$kWWO7q*ZYN=NOP~_Hh6;dP1eIH>`@9Ejrs0%;P|y*?=*npOy0h-I3VX~A7a?J z&0~LrQz#$Xs$g9wfRkMO%0NJnZBxKzJ$8S}RNc?clV;!z7a>s$rNWS1BJpPJVocw> zHcHDuM7qKD9Gs`QX={`*30`0ip6oG!p0DZ?4FgBZY!S|-xp7o;E$2>h_?sbxgP4~x zA{?Z1JFA8GA8>^l!tXqmze%Yq@$awLs9ETLO0p6E9j+cmLT6AOLZX%5(jLP1XnClY z`i(oS9Ex{eWDB`L`jy|QHKf*nzRhp=U>M=WUO}%255lG*`XiW9oSK{;yO(J8Y zS%xs;FyhU$p&qjGDnGNdoI7TQq3rI9dQshp404x2iCFy}pVo}OYi1pje@+ddXZCY< zd$M8^ggLq{Vv_dic`so&b!s$T-`&??JkDVY$FHT7zS+6r#+PVcxF=bTby=G)K2Q}= z#x-Ix9d#G#EzLtA*_P0%@QAnhaTn_@(rH6@*78ob%w0x9JBA^Ghcb=?@n+@OLaa;) z!I5}<2hlaW+hT#z&e|L`BltI`ny`yM1B^3YsL*suxaPiw-w661YselL$ z>D)flq9o@DuseQj_h;(jYhtEs@_6UN0x7p33z8wCl5ih#xC8 z1>WM07upO#FYhU4#4xf&DcH(s+kQeFVnItYW~bgu;n$LZd~Pt02*bj54YVo4OQmSE zELVD8`)(Q=Z)iC(>MKQ?16_^Clk#eHu2e>RDSZRdITvg4*q6!m*NbdNOy|y1n$#78 z?3;daFN{p%H{Nv7*1UcK;Sn2@GZECz06uuY?_tah({jq1YDON7m@__wJoW8(%i~G>D`S0=}r%uDDVU2WJz4aUo3oFjgw5 zCQOL)W${6R(Ti$9`_seRv$a-!OCqI*CHd|IMqm+F5c<1d*6|0xO&#e+r5{s1~#Y+o7bxEI9?c?fSDT+^Ob1pC)KiW-F!=`IDU|ArG<}D6|0!a zS9*|TRxr{kdvM!B1@$N2bIFxF=ZR`&ZA>@#3msU>5G}*P*|5$k*9-=I$-S*Sx^?0? z!!MbLw-SNPZeDU=A{c{wk z9ixxvGo5_GZC>Ia72U%?DR=e$oG)dBe}v{BiYJScp(*q?W@lD>-57lvp?QuK~imNiM z2c8TkAqFL9Nj6DlLjxY)B;`5np02rkcr4H~axJwjkUcW-dfyIu|2+a2kNL$*g<22N zzdG^n#*XKSv{ve7et2okmRBAG$@7;Spo-&!s+$@lW0~73bQ1*R`b=XXt;?V9Hv?H_ zR`Grcn)%1I;Q1f?x4ufBoQLytkv%TvnBqO#+EvK&mb>3;>i#fWZR6oe;e6vjP&|Tvs8F6M&4D{ z{|j2VM-$fmYMt7>Q`NI8Qlw|8Apt0#gZVoq0`@+X_8(Tqx8}W^$={^1EgX?of&4pJ zMbC=-(6CWrJ(>EQD5xeS$sSwaT*jPPNZvLyN|=Co16eoVfBvx{9Z^blsf2*9(Z9?W z!B2Vuuo1$>{y%vgWZp@}5qL~;yUvtZ{o5nvuUQNPkT?I>_up?&MY^##wmrf7XHsoI z9&CJ|I>h)t7Bk@ETE@W3MQvWl>CJz!pBR}3VB;Sxod^FtM#8Xh-~?kbSFPRi-(!?P zMSQ(RWRrM%Q0;LP z^vmpeM^Lup`~6(S&daeD#m4q}E`a?oY-v&MK3op^HP2wIZ3G4M4v|k;0!Xz=9k~+0 z%D9FB4iqD%_axwN{xU{i54r)M9|c;DVvv3%2r_|@0d`-mR}->p2(&JoPKu6;`$8EWQY)x&dh%X!FTyl~}CD{*8! zKk(j<1SaUMu={000Z*2J9|AVgjgSJED^7H&=$G#rIZk=0nOAuI?geL|vD0+YWP})( zI>|L`ST7+7fCahkmBT0X96q2M=xK#&#ry>740n*;dId&})*yqjw_YR#`GhgRTDqUB z;mGRS5LcqCM&LFY{b2mV-(}LNnXk91n;+{V#!q;DmE8jjN=4h;;`hY?aFvBCdyB4G z`QyX+k*O7A9%QF=?X_!%p-=%|7&4^^s0nJ|ISESOX__5++ckQCLskC&Jk&(Ex30sd z^Xh{F%tJ~*F1xlNhB%`^CW{^k@Z>$4CnhckL6*ai3Bdzo%e3ziBU)NjNpsnIHS;)z z**s4s-rbeO@+MVtv8K(r2>@yYwE8;1K4oN9{3Wlh3ceJPf)s|k)w`4@*wn6#+9Jrz zCZM4-Fat-Xd89>wux;f8jhFkl&Me z&YEq29YOB#vC28RY z5-Z(o#R=B>oxt&2^til=%b+8kRmOrybr%e4_#k_2_#M=f-_sq6qmRgK;FiS0%)LiPQhaA3qA*|c}3h6Z6#Wv1WEM3yCwBGG!<(^5H{;aC*rb{t0c=#htatG#cRf`H@ zm1(FBkSWDS>1sz#z<1v<3X{h{BXW(-pG>oT^SPolxhl9Ul?bom){Q)GPldP8)bf+k%sMWvNG{t%eSOr@0U;wo4s6K;Xo+jm} z-73CmiMiB~^S8nIb=MR)>C@KJ4agE;QAXYAbJ-!3%AbCm zI)n4fU6J7^%tDN76)LhTVj&c4&Me8i<1MGfws|m`t{gtmxCUs_a}u&8S1XLuJt!I` zGdhKpbeNti7M3J1q=D<8!!%_*glGpX#s11Mga@=#en86^=#HeLjnbLl)gkBk?~Nl0 zO#V2a57Nc3DY2;--fHDZRArhtG)x<{bCM~AwwIA&U!ZMFwNx*REA5;H|M7} z9S}CeJa@gzxn-8mHrJ+VR+QW!J11#%)yE>aj)Jm-mA2FSCi}yyCP&fP*c<+|fZ5U7 zGXwG2q6c-^a~1P@PNS!LADMYStqE(5I)Uq)C~tx_o1tT%SJ!PXoe@{HIc)66IXt<) zNaAjF0OPD+YYg0&GyS#`G;;nm9lpLcSz*e<+CDObB)mS2i5|D+oPlvv6ggV} z87BDkwfV0cF#w~Qyo0anIM8Nrmbv^O&VFVdl)w4sPa}68c#N}-mpK!@lx(XxBNycP zwgCvmH_sZtAo})cYe)t!`rA0Trk`HmZJj<{ESi4DMbazrdF~k%;(LBHrmO2 zIV{_|-5teWQccT`ucP!A;=7r6KZWWau=C)!T(_;4SPop5iDVKDs!4B6vY5xMm&CeX zT74Ty%`|6kG-|J|5}eDFtlsY23&y+L zcsB5-Ugvd0PH&{c8(;gt-qwTH=A0qor)uVNDz`}=k$+#zh8`Og?mjIM`)33qgC57A z6&QJ4Wc=RqpJ9yxJ(64&rDq3*ga6^7Wgy7t%WbL5?w`?@C^Gs&Q8D`G$A^&77v%B@ z_y3B%kkmB3@(`$!^fzon{?sxDf!6T+W2dZt?gt-i^pepb|7TQ$Kt^AS?b^ft+>a+R z`T|tue~%GFUzs4#LYyQT|F7sP0&Jv7P6q`@{u5vpz>(2ci`yjsKleieHvSBK$MxT1 z1T`g?0kI7aDu(`(4PU3Hj}pdUvio0wEH7xaja}-d{}sso-=_TU75_h@DSxHOT>^#< zgutG9#as+pxt#z4H@cRQo~yh~3!cp4ZcRrPBl}yu>z66&!l?yhBLn2$q`mhPy$y#e zIT~5>{H&3dX|2HSY0La8fJ^Qpeb_!iwrG*i%U_o?UbP+@wn^4pZyg!GiV`f=oJSIX zEuO+Zc3c9^SETQO4VF~KH#*IqbQH7YQ*^;2z=uHP4d9Bp?@_xfxA`~y0JL5=FcX-l z^fK|$f!(Hox=f5aAju!~nDH`X%3}Qp?6%v$B5@B0mmPNh*6{*%9OI?iHGCtyC{j`8V3TOqf)}`T4r(|yXEi(Q> z#>K3dYz=H+Ug5me`5gI?KU`hKR*y+B5SM7PTSk> z3(W>-nNUB|+W5iMX!v~*;C6qpH1J`rD-fPKK11STh5&&?vQxB3TNj=Q>;G4HPkbMo zEdikEgZ|b*L)Hx-H?mDR0NHP&=zUTJ1DGSCS0v)#vWr=V~Tojx!^a|b|6)^Soa-4K$jqh`yFUi=?ZHfm8x5TffZ zF{P;!gC#`@PFA;eq_3RSjG1qf7qWzxJCI+!0{(77>W{k6tw66NNElasByW}f;^Xgb zaL9DX@?QsF@3VC`HyBV{L8M#NBjDqO!R7QQt=wzT)$1E8Qe=*8No=q z31H*T*2Xn;DSarBF%QDDG19llKv`=Mppm)uAAsjt)q|UmBbM`1O*^=8}nntPViEzW7(j!f-UfXDg1wUnTyi>tbu;W;r4h@(*rD zZ#R(@`2aSDFP;7(KD`92C`p6f!mWRTWcfTXLdL;WeI8*|LjY+*5@G-jq(FKB4)@^R zJ8cG?hqnsIFIJGy>2*fP)}RkJ`SHcaJh;I$Yzg!q&;u3o#lwJmT!DJkuln5eA+e?khm1V*0C5=xES<# zUO&z?u)s$}+eWhp^=)sq8!=sh3OxDf3{ZL*3MBr-KX6maVwNScd)c23i-SjBCUB9> zh{d^LRb`3ka$cIlWm)BM`4ou71!Ln}>H&N)d*~OTQ8b?Ox3*klK5|S_5&XugU8i1D zl;#2eI~6?*Si;$u>JxCpOj2LcL~`FYNAr4++U~~0r8Y5gXsGBtNc--gq8F2Pc1Nu! znMvD2V|9QnmVm@rdY2YxQBuornE?_TL<|drcCelU{R)*c4!npV)fC(_gMdj*LweJx zG)m0#`6)6B$LeOQG_n{;2JYAVzMU7(G#GYbPN>0IMI%$Fr9G!R>r5FJAtWQsQO4Qo z3=>Iy5{;j8-kBBmOUBrBj$H9D%u74Hl0F7oOzXL)EKN)SzIWc`^_KC6Sy7c~ekR?w>=NL!uqhcC%?XXn|20^^sB z@D#fvw5PSZr^teUKHKNW0a;e{<(TcfE#VCBrz0KCHe6I$9wPmp83Z#)`fV=`i}Yu> zDRVY|g!tln61=U@QfxWJ&vNjP>~Xo@gkXCov1YxHWO!YMAEBY}5cd z6u9ioz3P=~kgIpe>h4!emo}>3kwS~Rn&}3Gf$Tf9$*m-SJ9`B5c1H}Yg3Dk9xRxuQR}yY_fng2sAz-{P?C-fYRYvbr{A8DTKZahNvU#` zMiuPTfw$liSVahIsx#Q$-SF1`vFX~CfjVATEq(EeQ6tg8mo2N_ea^PKWk)>apMnx=U*fBrDG|ty}DcI{ocI$#-qiS^X8YMBile6 z73<3d`7rqv*W(Y37A4rb0O-393PlQycC!R6);PiA3HHz4y&uydP}gf>AZ~njAC&5w zH5}J1)##IYZ5@*!6ALCpx2N{>J;vC(EUsAcp;te_;)g|0)6-R?WhHgz8D5EN=HP?f z{u-n>=;c>w+l;Ec!-p zc|X4d4{HBw+W1VPqGDXa;VIqn6Ru(3bv!fouDcOcjeIKX{CvrdTV98DxwQ$~y^mLK zI{Cx7hS415bJVa-`d{}R4sq5bNXeuBQRbWd^TTVqUB81VXF6h&wbBLmYK78BRFqU@ z+UHNes(qztBLtdas4Qg$kd%S-Pb_CMt}gSIPa`*!s5xFLwiv?%Z%`tVH|# zKbsSb)l&6w{c7u70j0RP zV&F2&(G)BZUq!VsA!vN%7Y6A~^Tr$xU>nLgsU$p@z$CY{4`cH4mpZB3*YVo};zK&A zINYLZ8x}(zebMA0*ntLz4V||;pk4|;b3&QOJ?Sfc8jV+sXng2YsKKdc7ut}sm=CDK zK}JtRDEO4ul;qYeRL}5N=RXJCPK=kafBWsF`*i|7u~{l+qZ!2Wbe?SEi0b#UuKb!0 z;`46D=&?;K@A+f5hCQ`&ZGuMT(sG|Mj{8y#>L&-l|e%S!ICZO6xHD zNto}_be$W=p^cLXYGN45y(NDw9xZpfpgMm>5!;0K0|oj%A~8&e7y|>N+dbk7@fi)9)4D8!IlMx6e8ybSf%hvD~4wM&W`a~vCR?C*hqz4}>2NZA!G z#tnVelxs3Hso+_K$&Pyh%cE{G0`6>l<5Y8yf&Fw8gHUC*Br>CdtbLBdAWBGIt?pK% zI?^un*$(Dda_Bc!ty`866IUH2b;^hrPO|RiqUnmfBcI%Kp~9`njY44fO0pHNk}+=k zaIK$PE8z}Hkqx_enj%XiOe(teC^+M%3aKP;;8z^tXd$?Kb(B?-|Htgehjao6Y+hIr zS2GK63Zs+`JZMjnMkU5Nv0d#_C0I(T3dPtS&uNZiCXa&3DCk#V;-gSTbC$3~ES}sPTm?v1a%~tC|%1Cx2+g~w!25U_b%{`}$7p8MH;}Yzp zflo&6Vh;1ZAUVqIQNaF1ZFuvZ?%G6|(9q#WB8E2xsli?$g$F(rj$iUule73#)TGXU zTs5ybVr=dRDs1k^4CnSkiEg<*k5OtkS1BJQ!o%s9B!oqJoIqxq-jAJ? z{dt(67ss(^Ug4*b>!GCxjJFDRKXjj!cx@@uF&ZXy z4hdy^9vWX&(Bysqs(T9sBEr7CvqES5{ed%2g4#YWX93gYh0oC4tLUKa_%PU(sdckWtK?Wl{?mPFr8Cd;b8F4c01rDKOXeF|#jz&#PFT!RE@Km~D zV|WHqK(s%K+EK!G7basU7Jl03Pl8@}w>I1`vMSuk2ebPdC_-zitcuvcR9{6GTj zQ3aj}mE2J1R^u&Afz2>k=RIO8(riME?n>RwPE&lYa&GQ&F0XS#|HTZhbXBE66*A-H zpE3mxiK&_IU@BWZH`Y;rl#acoS3am(_=p;0{sKo}vGE6)zkN^s zL~L|gr3zeDj}6W6A}hXQn4sYdJq;GH7nM01%I7hDeqXPbQJZwRnnJmgp-*nC+{@DJ z-zT|^m{{@xrOC(*C94w4<>HIBvAOm`*%OYTQ=K%;m>Zbe@6w8LBAA>T%YZU9S8IUHYPZn9F`hAWe0}Exrpy!}o2P*Ew1kz5DMInXfRpnr~Wg z4URhdB(5TA5+;Vah#6vQ=_O;Ja1f>wah362>EJ#J%jubziO>$f6)^CgItYe~-T5y+!l=9)=8AGToeT>XxH7 zn2WzG-+7aRKJm^(b+gt>hN0+afZs?xD|CggqHBL@Q9dkfw9GYyn{M&!`-1lwdC{WX zMuH>ZhDI{GV9Cy5-+J#u&VJBKqV7nr5@LF?J?)5I#iJrX5HfUPc@qVXM)gov_V{$n z4}Tmk-YF%TTA?47mGt{TW81Fw#Tnke?^P}EC92FpUHcuoncohe&wS8cn`@q^-0=Ed<7`kO z+g8~(dS@3c*_VEY{(?=JxVpCw<9W{H@~M-`qr8I2i=VVNrkuIN>N+?TqZ3@TVNL0p z?!Tvpgt#hf#8h$%}N7MedPeEL*GVhepku%&iSiPpG~{D6Q7 zhma3HT!#|XINnqgpJBdlv(76@cq2J5`io?kku&yrij=N6#_Je6rd?kbkrq$`el|}@ z_bgrVG4*}8X`<2v4DDyQf|`DHV%$%DWu8aBhV#9jEyN`UrTS5b9KoX&{(@80N5MjGnDYAXq{$y| zoWJkTeoLB^V$PAL8TQKBi`??yyQG%XfQQ!`o$$Hvd3y13w1lT@C)|E)>9fQ+KZ7T& z%Hf+AO+4}(LJ~70Z(<%ZXoZnz7vpy?HU3g`4QeCqN+RLYhH5Y?lkBa0(G*LT4|tbM zXFP|As@!sUh32FL{hkrDA6v%sQMa|`71JHtq#l0Ag=iqQ<0p{7K6$W@s=b=tlkfMz z*qwv_9K+}>9_n2tv~k_yKTMY%!!*i;-tW6Vt&d5OmQUU*dyKMTLX9!^{)A@ckD!T` zD3F3hE1Hr#=$~cXQqwbhHo}S3uwg}-n5Yv<~ z;!WN=MJKi7d5pLqZ?;7iw6LL!u|K8%PCTidXTE(c!OPFyH1p=1NYo1V0vK>Sm(6K1 zl((qt+HYDM-ftQ&?_zq$`)loiL`G@WEejz|V_3Of#483i+;^CkEo?ocsCNs|PnCiQ z-(V2oLc3R{mIYhPqIHKlmPQG`V7yxZDq*5fNBOS;8bk%Gbd&TgiH7K?grhZiW0Noa8BvVr7OhjF5=c$MGcxOxzkjH6Y2u4M zK>+!#(F2(oM^CaPqoF~UbTkIx4MI0VD5xb^rYIfW#F6}mKEL$L%yLV+Pz!H7GN9Vs zb1x;CeH|nGNyG9ot^6=8N^4qrCxr2%uE0ov*<)Nrf+%Xlk=Gfi(Xx6uy)l^x(K%Kh z+jp9q-eJWESXa3NTGXANi{R0eEOc!gSBhiH(dhF0rHG7~Ncva3*86jE*t@#UXqfbK zqTR1v8U(!N9h`8H-57h_7}C=8crK=j;?=jZf_Qq1&H`pV>ax4;ZWQUb#B2s^JJ|1Q zF)djPa!;bsmeazEW3N7dxiz9n=V2x_BnS;fQQTwvBtzf#ibs}q<9E~I=pUyh`Oyx# zQCECR5o`F5bko8gg1I+0x1ybI&?RVwn8ohKF@NS$hF-d#>==nFRdXJEW5GzcFOyLy zje-^Jfmfzmc5fRxO!!$k0YY&rlY7VDJ04%tAtsmSkr2xdq^h1$U8cbyQFK{6?S7o- z+I@tlqae=4ouxt$k9ZTV54i7q=QBZAIrsB*%jUW+8oo%bw``Ye#jG`Z5qmHj=6=!v z1j;TL{Oox@6OxwFt7UMheG!-mX_bae1K8zAH&oPIFdNxF4-pS#PoAy<;R?`lQ zo4m7Q_R`DBpfTr+e8h!IhkS_)nn3^3*rrYV%I~~5Cl~)U+KvJ%b1(m-y3k&BZf5Am zhv8J#zT`KmuA#!+FH5?g)3|cxy-#0S0aObm`i6kRv~x9Y{UY9(_z$qq4P{mH<@$4| zM2}^%P=A+Ot^9!%*YJXjE#Gd6)Ur6vp@eqlPK=VtUZW&F;V-#moSVp@{-Se82{|q3{}zqHAHELNW6`qa-dY`Ct)OenS>2ReGIHV^GQRX!y&EQ!*zn9 zad#)F!`8! z^Zx`>1+4lyMKVw~@a+*Sr5REBwG>@1iirE6zC@df{RonHc{}|kv z)zBcYY!SV=DPd`q+1!*^6j>ltnEg1_1|Dy${+u7>iZJlaPNdiA;5bmLmsSbJ+JL`1 z=Q-K|!iYA-ADER|hW4bQu!|RTf-;-77y1bS8|159hGJTzM;p89s;ko1`}r0)z-0OO z@#7|=7*LpKuy^_xAt!RV$2&EsqddLZ`}uSHp4tve@Ahk*1Dyl=&jEt*O2n!&kY(sF z=rZXlv;3#yD;O{^nvKD_#TWiLDdmg2sEtQ*Eg0k~WfofWAh{?Ea&>5qIpjgv<=bmU zMa54o4%Tv(LJO1SY&^(M!`MjW;6eb&3EJRPuuLG}S$cJ9L7Mn8=@B81g2F5VqA@Uw z0p2T}*`)7GIj2J0a>oosey{dDe3Jst%!)iEozgS%moscWW?)mHRZ8jfYX`@Hscm5D zbIh)pY~fx5$)DQ<6O*_W0cG>aC!bjV!9s_15UVarI`T!26D87RW?Y+Wweslq+8kh_ zPP93Ko(cE8_ud<2TPip5RFF+Xo+Cq$ee$ETN4QRV$f1t7nWk5JKhKBXQ`=$b-F~fe zpmSii96)zrlFVc`HOgs6JcEK$0N3!r2OkVK-E@=P$8C+@5Mv334kq1DDnsl-;a_V! zdjI|R&E1ZbDFvjN3k14cYs|ql0{1yq4AGZRmXUi)qFFw)_*lrX*?dD_1WKn@49rL1 zVr}AF+Ze{BQJ_&yo;+zgt56YW_AK|;Tn9+kAwQ@@Oz1QQ+yVf#UYCGR!AAuEX&M$j zi&Q!aIR3nfAYM0V^Cu=v{DfJ`!7qi~X*w;J15?|;)aPj1EZxvnS-M6D(cVy$S_BLc zN-ctpO*?pGy*^!Re39O<*ii$c#DnnKhOMbM*BVsdz$YL=X$SrggbSvakZOj2e zkFgD139GUY4-C33$#UnX{~}6lZYg*LaBHlxSZx`^-Decb5uAZ=A8T-fP`s7KBePmCbi(4{m!LMRXIw3QZ6ZNdKnA!%WKF610rN=6Za)hEZm6##ZE?Fu{3FQihd?t&1l9Y00`{!6$kyWWrGr6R4 zBhQryF|JpsAGbe)D+BeH-%HbLz3=JxaZWq*clxn&pmSjN9ANbaMG}`Z^aO}tVODoB z@^L?x6BsejQ9M9C{fbnEq}TKF+IWO)B6tpj8m*$-E(fLLfp177>`KpUdW{*ul9gG$?*YMa zFzW*C03{ax;;zhYTJ!z=&VkN>&VkN>&Vj?t0amVYqs&)qPn@hpc.msu.edu + +use your MSU NetID and log into the machine 'gateway.hpcc.msu.edu'. + +Then log into a compute node e.g. dev-intel14-phi:: + + ssh dev-intel14-phi + +Assing the path of your working directory to a ageneric variable:: + + workingPath=$"" + +.. :: + workingPath=$"/mnt/ls12/Tamer/refTransProject/" + +2. Getting the data +------------------- + +We'll be using a few RNAseq data sets from Fagerberg et al., `Analysis +of the Human Tissue-specific Expression by Genome-wide Integration of +Transcriptomics and Antibody-based Proteomics +`__. + +You can get this data from the European Nucleotide Archive under +`ERP003613 `. +All samples in this project are `paired end +`__ . +So each sample is represented by 2 files. These files are in +`FASTQ Format `__ . + + +In this tutorial we will work with two tissues: `salivary gland +`__ and `lung +`__. Note that each +tissue has two replicates, and each replicate has two files for the paired end reads + +.. check if the data directory exists +.. :: + + if [ ! -d "$workingPath/refTransData" ]; then + +make a directory to host the data on the HPC +:: + + mkdir $workingPath/refTransData + cd $workingPath/refTransData + +Download the sequencing files of the salivary gland +:: + + wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR315/ERR315325/* + wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR315/ERR315382/* + +Download the sequencing files of the lung tissue +:: + + wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR315/ERR315326/* + wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR315/ERR315424/* + + +.. :: + + fi + +Now do:: + + ls -la $workingPath/refTransData/ + +You'll see something like :: + + -r--r--r-- 1 mscholz common-data 3714262571 Dec 4 08:44 ERR315325_1.fastq.gz + -r--r--r-- 1 mscholz common-data 3714262571 Dec 4 08:44 ERR315325_2.fastq.gz + ... + +which tells you that this file is about 900 MB. You can go a head and use these files +for the reset of the protocl. But for the sake of time (& memory), we will run our demo +on a subset of this dat + +3. Copying in some data to work with. +------------------------------------- + +.. check if the subset data directory exists +.. :: + + if [ ! -d "./dataSubset" ]; then + +First, make a directory +:: + + mkdir ./dataSubset + +Copy in a subset of the data (100,000 reads) +:: + + gunzip -c ERR315325_1.fastq.gz | head -400000 | gzip > dataSubset/salivary_repl1_R1.fq.gz + gunzip -c ERR315325_2.fastq.gz | head -400000 | gzip > dataSubset/salivary_repl1_R2.fq.gz + gunzip -c ERR315382_1.fastq.gz | head -400000 | gzip > dataSubset/salivary_repl2_R1.fq.gz + gunzip -c ERR315382_2.fastq.gz | head -400000 | gzip > dataSubset/salivary_repl2_R2.fq.gz + +and do the same for the lung samples +:: + + gunzip -c ERR315326_1.fastq.gz | head -400000 | gzip > dataSubset/lung_repl1_R1.fq.gz + gunzip -c ERR315326_2.fastq.gz | head -400000 | gzip > dataSubset/lung_repl1_R2.fq.gz + gunzip -c ERR315424_1.fastq.gz | head -400000 | gzip > dataSubset/lung_repl2_R1.fq.gz + gunzip -c ERR315424_2.fastq.gz | head -400000 | gzip > dataSubset/lung_repl2_R2.fq.gz + + +.. :: + + fi + +4. Link your data into a working directory +------------------------------------------ +Rather than copying the files into the working directory, let’s just link them in. +this creates a reference so that UNIX knows where to find them but doesn’t need to actually move them around. +:: + ln -fs $workingPath/refTransData/dataSubset/*.fq.gz $workingPath + +4. install/load software +------------------------ +Most of softwares used through this protocol require to be installed on the root dircteory of HPC. +This action requires adminstrative privilages. The IT supervisors of the HPC install software packages +in the form modules that you can simply load by the command: module load + +Now let us load all the modules reuired to go through this protocol +:: + + module load FastQC/0.11.2 + module load Trimmomatic/0.32 + +We will elaborate more on each software when we use it + +.. shell stop + +Next: :doc:`m-quality` diff --git a/refTrans/m-quality.txt b/refTrans/m-quality.txt index d9b7aa7..8337ce1 100644 --- a/refTrans/m-quality.txt +++ b/refTrans/m-quality.txt @@ -1,205 +1,95 @@ Short read quality and trimming =============================== -Log into the HPC with SSH; use your MSU NetID and log into the machine -'gateway.hpcc.msu.edu'. There copy/paste:: +.. shell start - module load powertools - getexample RNAseq-model - -and you should see something about linking. - -0. Getting the data -------------------- - -We'll be using a few RNAseq data sets from Fagerberg et al., `Analysis -of the Human Tissue-specific Expression by Genome-wide Integration of -Transcriptomics and Antibody-based Proteomics -`__. - -You can get this data from the European Nucleotide Archive under -ERP003613 -- go to http://www.ebi.ac.uk/ena/data/view/ERP003613 -to see all the samples. - -I picked two samples: `salivary gland -`__ and `lung -`__. Note that each -sample has two replicates, and each replicate has two files. - -**Don't download them**, but if you were downloading these yourself, -you would want the "Fastq files (ftp)", both File 1 and File 2. (They -take a few hours to download!) - -We've already loaded the data onto the MSU HPC, and you've loaded -it with 'module load powertools'. - -To log into a compute node, type:: - - ~/RNAseq-model/login.sh - -If this doesn't work, do:: - - ssh dev-intel14-phi - -Now do:: - - ls -la ~/RNAseq-model/data/ - -You'll see something like :: - - -r--r--r-- 1 mscholz common-data 3714262571 Dec 4 08:44 ERR315325_1.fastq - -r--r--r-- 1 mscholz common-data 3714262571 Dec 4 08:44 ERR315325_2.fastq - -r--r--r-- 1 mscholz common-data 2365633645 Dec 4 08:44 ERR315326_1.fastq - -r--r--r-- 1 mscholz common-data 2365633645 Dec 4 08:44 ERR315326_2.fastq - -which tells you that this file is 900,000,000 bytes or about 900 MB. -Quite large! - -1. Copying in some data to work with. -------------------------------------- - -First, make a directory:: - - mkdir ~/rnaseq - cd ~/rnaseq - -Copy in a subset of the data (100,000 reads):: - - head -400000 ~/RNAseq-model/data/ERR315325_1.fastq | gzip > salivary_repl1_R1.fq.gz - head -400000 ~/RNAseq-model/data/ERR315325_2.fastq | gzip > salivary_repl1_R2.fq.gz - -These are FASTQ files -- let's take a look:: - - less salivary_repl1_R1.fq.gz - -(type 'q' to exit less) - -Question: - -* why are some files named ERR*? -* why are some files named salivary*? -* why is there R1 and R2 in the name? - -Links: - -* `FASTQ Format `__ - -2. FastQC +1. FastQC --------- -We're going to use `FastQC `__ to summarize the data. - -First, we need to load the FastQC software into our account:: +We're going to use `FastQC `__ to summarize the data. +FastQC cans on Fastq files as well as GZip compressed Fastq (.gz) files. The program operates on one file at a time - module load FastQC/0.11.2 +Now run FastQC on your sequencing file +:: + + mkdir rawData_fastqc + fastqc *.fq.gz --outdir=./rawData_fastqc -(You have to do this each time you log in and want to use FastQC.) +.. note:: When you pass multiple files to FastQC, it operates a queueing system +where only one file is opened at a time, and rest of the files will wait until existing +files have been processed. If your files are big in size, it might be more wise +to run then in parallel. -Now, run FastQC on both of the salivary gland files:: - fastqc salivary_repl1_R1.fq.gz - fastqc salivary_repl1_R2.fq.gz +Now type 'ls' in the output directory to explore the result files +:: -Now type 'ls':: - - ls + ls ./rawData_fastqc and you will see :: salivary_repl1_R1_fastqc.html salivary_repl1_R1_fastqc.zip - salivary_repl1_R2_fastqc.html - salivary_repl1_R2_fastqc.zip - -Copy these to your laptop and open them in a browser. If you're on a -Mac or Linux machine, you can type:: - - scp username@hpc.msu.edu:salivary*fastqc.* /tmp - -and then open the html files in your browser. For Windows, we'll -figure it out ;). + ... -You can also view my versions: `salivary_repl1_R1_fastqc.html -`__ -and `salivary_repl1_R2_fastqc.html -`__ +Copy these to your laptop and open them in a browser:: -Questions: + scp :/rawData_fastqc/salivary*fastqc . -* What should you pay attention to in the FastQC report? -* Which is "better", R1 or R2? +and then open the html files in your browser. Links: -* `FastQC `__ * `FastQC tutorial video `__ -3. Trimmomatic +2. Trimmomatic -------------- -Now we're going to do some trimming! We'll be using -`Trimmomatic `__. -For -a discussion of optimal RNAseq trimming strategies, see `MacManes, 2014 `__. - -First, load the Trimmomatic software:: - - module load Trimmomatic/0.32 - -Next, run Trimmomatic:: - - java -jar $TRIM/trimmomatic PE salivary_repl1_R1.fq.gz salivary_repl1_R2.fq.gz\ - salivary_repl1_R1.qc.fq.gz s1_se salivary_repl1_R2.qc.fq.gz s2_se \ - ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 \ - LEADING:2 TRAILING:2 \ - SLIDINGWINDOW:4:2 \ - MINLEN:25 - -You should see output that looks like this:: - - ... - Input Read Pairs: 100000 Both Surviving: 95236 (95.24%) Forward Only Surviving: 4764 (4.76%) Reverse Only Surviving: 0 (0.00%) Dropped: 0 (0.00%) - TrimmomaticPE: Completed successfully - -Questions: - -* How do you figure out what the parameters mean? -* How do you figure out what parameters to use? -* What adapters do you use? -* What version of Trimmomatic are we using here? (And FastQC?) -* Are parameters different for RNAseq and genomic? -* What's with these annoyingly long and complicated filenames? -* What do we do with the single-ended files (s1_se and s2_se?) - -Links: - -* `Trimmomatic `__ - -4. FastQC again +Now we're going to do some trimming to remove any leftover adaptor sequences and low quality trailing edges. +We'll be using `Trimmomatic `__. +For a discussion of optimal RNAseq trimming strategies, see `MacManes, 2014 +`__. + +To run Trimmomatic +:: + + for f in *R1*.gz; do + input_R1="$f" + input_R2=$(echo "$f" | sed s/R1/R2/) + output_R1_PE=$(echo "$f" | sed s/R1/R1.pe/) + output_R1_SE=$(echo "$f" | sed s/R1/R1.se/) + output_R2_PE=$(echo "$f" | sed s/R1/R2.pe/) + output_R2_SE=$(echo "$f" | sed s/R1/R2.se/) + java -jar $TRIM/trimmomatic PE $input_R1 $input_R2 \ + $output_R1_PE $output_R1_SE $output_R2_PE $output_R2_SE \ + ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 \ + LEADING:2 TRAILING:2 \ + SLIDINGWINDOW:4:15 \ + MINLEN:25 &>> log.txt + # Titus workshop: ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 LEADING:2 TRAILING:2 SLIDINGWINDOW:4:2 + # my parameters: ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:30:10:1 SLIDINGWINDOW:4:15 + done + +You can check the log file to see how much the trimming affected the data + +For better understanding of the Trimmomatic parameters, you need to check +`Trimmomatic web page `__ . +You might need to teak these parameters for you own data + +.. * Are parameters different for RNAseq and genomic? +.. * What do we do with the single-ended files (s1_se and s2_se?) + +3. FastQC again --------------- -Run FastQC again:: - - fastqc salivary_repl1_R1.qc.fq.gz - fastqc salivary_repl1_R2.qc.fq.gz - -(Note that you don't need to load the module again.) - -Copy them to your laptop and open them, OR you can view mine: `salivary_repl1_R1.qc_fastqc.html -`__ -and `salivary_repl1_R2.qc_fastqc.html -`__ - -Let's take a look at the output files:: - - less salivary_repl1_R1.qc.fq.gz +Run FastQC again +:: -(again, use 'q' to exit less). + mkdir trimmedData_fastqc + fastqc *.pe.fq.gz --outdir=./trimmedData_fastqc -Questions: +.. * Does it matter that you still have adapters!? -* Why are some of the reads shorter than others? -* is the quality trimmed data "better" than before? -* Does it matter that you still have adapters!? +.. shell stop Next: :doc:`m-tophat` diff --git a/refTrans/model-orgs.txt b/refTrans/model-orgs.txt index a760bf4..b134ea8 100644 --- a/refTrans/model-orgs.txt +++ b/refTrans/model-orgs.txt @@ -1,16 +1,14 @@ -2014 / December / mRNAseq on model organisms -============================================ +Reference based RNAseq assembly in model organisms +================================================== -This workshop was given on December 4th and 5th, 2014, by -C. Titus Brown and Matt Scholz. - -.. figure:: _static/model-rnaseq-pipeline.png +.. figure:: files/model-rnaseq-pipeline.png Tutorials: .. toctree:: :maxdepth: 2 + m-dataNsoftware_HPC m-quality m-tophat m-htseq @@ -18,4 +16,4 @@ Tutorials: m-data-analysis m-hpc-qsub m-advice - more-resources +.. more-resources From 6bef24e14c65b43c2601d8eb9703744a5e898721 Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Tue, 3 Feb 2015 03:21:52 -0600 Subject: [PATCH 05/14] 1. modify and apply literate-resting code to m-tophat. 2.fix typos and design issues in the other files --- refTrans/m-dataNsoftware_HPC.txt | 27 +++--- refTrans/m-quality.txt | 42 ++++++--- refTrans/m-tophat.txt | 156 +++++++++++++++++++++++-------- refTrans/model-orgs.txt | 2 +- 4 files changed, 164 insertions(+), 63 deletions(-) diff --git a/refTrans/m-dataNsoftware_HPC.txt b/refTrans/m-dataNsoftware_HPC.txt index ad70fcc..ebff038 100644 --- a/refTrans/m-dataNsoftware_HPC.txt +++ b/refTrans/m-dataNsoftware_HPC.txt @@ -1,10 +1,10 @@ -Get an example data and install/load software on HPC +Get example data and install/load software on HPC ==================================================== .. shell start -1. Log into HPC ---------------- +Log into HPC +------------ Log into the HPC with SSH:: @@ -23,8 +23,8 @@ Assing the path of your working directory to a ageneric variable:: .. :: workingPath=$"/mnt/ls12/Tamer/refTransProject/" -2. Getting the data -------------------- +Getting the data +---------------- We'll be using a few RNAseq data sets from Fagerberg et al., `Analysis of the Human Tissue-specific Expression by Genome-wide Integration of @@ -86,8 +86,8 @@ which tells you that this file is about 900 MB. You can go a head and use these for the reset of the protocl. But for the sake of time (& memory), we will run our demo on a subset of this dat -3. Copying in some data to work with. -------------------------------------- +Copying in some data to work with +--------------------------------- .. check if the subset data directory exists .. :: @@ -120,15 +120,15 @@ and do the same for the lung samples fi -4. Link your data into a working directory ------------------------------------------- +Link your data into a working directory +--------------------------------------- Rather than copying the files into the working directory, let’s just link them in. this creates a reference so that UNIX knows where to find them but doesn’t need to actually move them around. :: ln -fs $workingPath/refTransData/dataSubset/*.fq.gz $workingPath -4. install/load software ------------------------- +install/load software +--------------------- Most of softwares used through this protocol require to be installed on the root dircteory of HPC. This action requires adminstrative privilages. The IT supervisors of the HPC install software packages in the form modules that you can simply load by the command: module load @@ -138,9 +138,14 @@ Now let us load all the modules reuired to go through this protocol module load FastQC/0.11.2 module load Trimmomatic/0.32 + module load TopHat2/2.0.12 +.. note:: Loading TopHat2 include loading Bowtie/2.2.3.0 and Samtools/0.1.19.0 + We will elaborate more on each software when we use it .. shell stop +---- + Next: :doc:`m-quality` diff --git a/refTrans/m-quality.txt b/refTrans/m-quality.txt index 8337ce1..6652268 100644 --- a/refTrans/m-quality.txt +++ b/refTrans/m-quality.txt @@ -3,8 +3,8 @@ Short read quality and trimming .. shell start -1. FastQC ---------- +Assessment of sequence quality (FastQC) +--------------------------------------- We're going to use `FastQC `__ to summarize the data. FastQC cans on Fastq files as well as GZip compressed Fastq (.gz) files. The program operates on one file at a time @@ -21,8 +21,7 @@ files have been processed. If your files are big in size, it might be more wise to run then in parallel. -Now type 'ls' in the output directory to explore the result files -:: +Now type 'ls' in the output directory to explore the result files:: ls ./rawData_fastqc @@ -32,9 +31,9 @@ and you will see :: salivary_repl1_R1_fastqc.zip ... -Copy these to your laptop and open them in a browser:: +In the root of your terminal, you can copy these files to your laptop and open them in a browser:: - scp :/rawData_fastqc/salivary*fastqc . + scp :/rawData_fastqc/salivary*fastqc.* /tmp and then open the html files in your browser. @@ -42,8 +41,8 @@ Links: * `FastQC tutorial video `__ -2. Trimmomatic --------------- +Trimmomatic +----------- Now we're going to do some trimming to remove any leftover adaptor sequences and low quality trailing edges. We'll be using `Trimmomatic `__. @@ -70,7 +69,26 @@ To run Trimmomatic # my parameters: ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:30:10:1 SLIDINGWINDOW:4:15 done -You can check the log file to see how much the trimming affected the data +Now type 'ls' in the output directory to explore the result files:: + + ls *.fq.gz + +and you will see :: + + salivary_repl1_R1.fq.gz + salivary_repl1_R1.pe.fq.gz + salivary_repl1_R1.se.fq.gz + salivary_repl1_R2.fq.gz + salivary_repl1_R2.pe.fq.gz + salivary_repl1_R2.se.fq.gz + ... + +Trimming of paired end reads in four output files for each pair input files. +The two output files labeled with (.pe) represnt the reads that are still in the pair end format. +The other two files labled with (.se) represent the unpaired reads remaining from one input file +after dropping the corresponding read from the second input file. + +You can check the log file to see how much the trimming affected the data. For better understanding of the Trimmomatic parameters, you need to check `Trimmomatic web page `__ . @@ -79,8 +97,8 @@ You might need to teak these parameters for you own data .. * Are parameters different for RNAseq and genomic? .. * What do we do with the single-ended files (s1_se and s2_se?) -3. FastQC again ---------------- +FastQC again +------------ Run FastQC again :: @@ -92,4 +110,6 @@ Run FastQC again .. shell stop +---- + Next: :doc:`m-tophat` diff --git a/refTrans/m-tophat.txt b/refTrans/m-tophat.txt index e3bb17b..9c68764 100644 --- a/refTrans/m-tophat.txt +++ b/refTrans/m-tophat.txt @@ -1,65 +1,141 @@ Mapping reads to the transcriptome with TopHat ============================================== +.. shell start + Now that we have some quality-controlled reads, we're going to *map* the reads to the reference gene set, for the purpose of counting how many reads have come from each gene. We'll be using the `TopHat software `__ -For this purpose, we've already installed the human reference gene set -on the HPC (as part of the data you loaded at the beginning). In this -case we've loaded in the `Illumina iGenomes project -`__ into the RNAseq-model -data set. +Download the Bowtie2 index +-------------------------- + +TopHat is built on the short read mapping program `Bowtie +`__. Bowtie indexes the genome +with a Burrows-Wheeler index to keep its memory footprint small. +We can build those indexes for our genome using `bowtie2-build +`__. +However for most if not all model organisms, you can download the Bowtie index for +a given genome assembly. + +The `Illumina iGenomes project `__ +provided the RNA-Seq user community with a set of genome sequence indexes (including Bowtie indexes). + +We will download the iGenome index for the UCSC hg19 assembly version and +define the path of the genome and genome_index_base which the name of any +of the index files up to but not including the first period. + +.. check if the data directory exists +.. :: + + if [ ! -d "$workingPath/Homo_sapiens" ]; then + +:: + + wget ftp://igenome:G3nom3s4u@ussd-ftp.illumina.com/Homo_sapiens/UCSC/hg19/Homo_sapiens_UCSC_hg19.tar.gz + tar -zxvf Homo_sapiens_UCSC_hg19.tar.gz + homo_sapiens_genome="$workingPath/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa" + Bowtie2Index="$workingPath/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome" + +.. :: + + fi + +Making use of gene model annotations +------------------------------------ +Model organisms usually have a good annotation of their transcriptome. Annotation +information is typically shared in GTF files. Illumina iGenomes indexes come +with GTF transcript annotation files (Congratualtion, this means you have downloaded them already) +Let us define the path of the required GTF in the downloaded iGenome package +:: + GTF_file="$workingPath/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf" + +We can pass the annotation file to TopHat through the optional argument (--GTF). +TopHat will first extract the transcript sequences and use Bowtie to align reads to this virtual transcriptome first. +Only the reads that do not fully map to the transcriptome will then be mapped on the genome. + +For TopHat to make use of GTF file, it has to create a Bowtie index for the transcriptome +sequences which could be time consuming. If you are aligning multiple samples with TopHat, +we can create these indexes once and reused for multiple and even *simultaneous* TopHat runs. +To create the transcriptome-index, we can inovke TopHat without input reads +:: + + tophat \ + --num-threads 4 \ + --GTF "$GTF_file" \ + --transcriptome-index "$workingPath/trans_Index/transcriptome" \ + $Bowtie2Index + + +Prepare the input reads +----------------------- +You can pass these two PE files directly to the TopHat to do the pair end orinted alignment. +TopHat allows the use of additional unpaired reads but it has to pass as single file. +So we need to merge every two SE files belonging to the same sample +:: + + for f in *R1.se*; do + input_R1="$f" + input_R2=$(echo "$f" | sed s/R1/R2/) + output=$(echo "$f" | sed s/R1/R/) + printf $input_R1" "$input_R2" "$output"\n" + cat $input_R1 $input_R2 > $output + done Mapping reads ------------- -Load the TopHat software:: - - module load TopHat2/2.0.12 - -And now run TopHat:: - - cd ~/rnaseq - tophat -p 4 \ - -G ~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Annotation/Genes/genes.gtf \ - --transcriptome-index=$HOME/RNAseq-model/transcriptome \ - -o tophat_salivary_repl1 \ - ~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Sequence/Bowtie2Index/genome \ - salivary_repl1_R1.qc.fq.gz salivary_repl1_R2.qc.fq.gz - -This will take about 15 minutes. +And now run TopHat on one sample:: + + tophat \ + --num-threads 4 \ + --GTF "$GTF_file" \ + --transcriptome-index "$workingPath/trans_Index/transcriptome" \ + --output-dir salivary_repl1 \ + $Bowtie2Index \ + salivary_repl1_R1.pe.fq.gz salivary_repl1_R2.pe.fq.gz,salivary_repl1_R.se.fq.gz + +This will take about ~10-15 minutes. We can run all the samples sequencially using a for loop +:: + + 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 -d "_") + printf $input_R1"\n"$input_R2"\n"$input_R"\n"$output"\n" + tophat \ + --num-threads 4 \ + --GTF "$GTF_file" \ + --transcriptome-index "$workingPath/trans_Index/transcriptome" \ + --output-dir $output \ + $Bowtie2Index \ + $input_R1 $input_R2,$input_R + done -Questions: - -* What are all these parameters?! -* How do we pick the transcriptome/genome? -* Why is it so slow? -* What is different about mapping RNAseq reads vs mapping genomic reads. - -Links: - -* `TopHat manual `__ -* `Illumina iGenomes project `__ Counting mapped reads percentage -------------------------------- Let's ask samtools for the total number of reads that mapped:: - samtools view -c -F 4 tophat_salivary_repl1/accepted_hits.bam + samtools view -c -F 4 salivary_repl1/accepted_hits.bam + +You should get around 186,256. + +We can count the reads in our input files:: + + R1_pe="$(zcat salivary_repl1_R1.pe.fq.gz | echo $((`wc -l`/4)))" + R2_pe="$(zcat salivary_repl1_R2.pe.fq.gz | echo $((`wc -l`/4)))" + R_se="$(zcat salivary_repl1_R.se.fq.gz | echo $((`wc -l`/4)))" + echo $((R1_pe + R2_pe + R_se)) -You should get around 179,312. +You should get around 189,150 -If we look at the `FastQC report -`__, -we can see (at the top) that the total number of reads in the R1 file -is 95,236 (which is the same number of reads as in the R2 file, because -we only looked at the paired reads that came out of Trimmomatic). +So you have 186,256 mapped reads out of 189,150 -- about 98%. That's pretty good! -In total, that's 190,472 reads in your trimmed (QC) data, and 179,312 -mapped -- about 94%. That's pretty good! +.. shell stop ---- diff --git a/refTrans/model-orgs.txt b/refTrans/model-orgs.txt index b134ea8..625d4d1 100644 --- a/refTrans/model-orgs.txt +++ b/refTrans/model-orgs.txt @@ -6,7 +6,7 @@ Reference based RNAseq assembly in model organisms Tutorials: .. toctree:: - :maxdepth: 2 + :maxdepth: 1 m-dataNsoftware_HPC m-quality From 39df3845baea98554239523f84bd2fe0522d3263 Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Tue, 3 Feb 2015 04:34:41 -0600 Subject: [PATCH 06/14] 1.fix the cut function in m-tophat. 2.make the para_HPC.txt for parallel qsub of TopHat jobs --- refTrans/m-tophat.txt | 7 ++++++- refTrans/model-orgs.txt | 1 - refTrans/para_HPC.txt | 37 +++++++++++++++++++++++++++++++++++++ 3 files changed, 43 insertions(+), 2 deletions(-) create mode 100644 refTrans/para_HPC.txt diff --git a/refTrans/m-tophat.txt b/refTrans/m-tophat.txt index 9c68764..a05b714 100644 --- a/refTrans/m-tophat.txt +++ b/refTrans/m-tophat.txt @@ -103,7 +103,7 @@ This will take about ~10-15 minutes. We can run all the samples sequencially usi 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 -d "_") + output=$(basename "$f" | cut -f 1,2 -d "_") printf $input_R1"\n"$input_R2"\n"$input_R"\n"$output"\n" tophat \ --num-threads 4 \ @@ -114,6 +114,11 @@ This will take about ~10-15 minutes. We can run all the samples sequencially usi $input_R1 $input_R2,$input_R done +With large input files, TopHat will take hours and may be days. It will be more +efficient to run TopHat on all samples in parallel jobs + +:doc:`./para_HPC + Counting mapped reads percentage -------------------------------- diff --git a/refTrans/model-orgs.txt b/refTrans/model-orgs.txt index 625d4d1..68204f8 100644 --- a/refTrans/model-orgs.txt +++ b/refTrans/model-orgs.txt @@ -14,6 +14,5 @@ Tutorials: m-htseq m-more-tophat m-data-analysis - m-hpc-qsub m-advice .. more-resources diff --git a/refTrans/para_HPC.txt b/refTrans/para_HPC.txt new file mode 100644 index 0000000..7841b24 --- /dev/null +++ b/refTrans/para_HPC.txt @@ -0,0 +1,37 @@ +Submitting jobs to the MSU HPC queue +==================================== + +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 "_") + script="./qsub_scripts/"$output"_script.sh" + printf $input_R1"\n"$input_R2"\n"$input_R"\n"$output"\n" + 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 + + +You can use 'qstat | grep ' to check on your jobs' status. + +Back: :doc:`m-tophat` From ef8fb3a087d7353d8588f8cc3ae6cb83a38d2d26 Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Fri, 6 Feb 2015 16:22:18 -0600 Subject: [PATCH 07/14] draft version of refTrans for model organisms --- refTrans/index.txt | 5 +- refTrans/m-HPC-cuffdiff.txt | 62 +++++++++ refTrans/{para_HPC.txt => m-HPC-tophat.txt} | 19 ++- refTrans/m-advice.txt | 56 +++----- refTrans/m-count.txt | 58 ++++++++ refTrans/m-cuffdiff.txt | 56 ++++++++ refTrans/m-cufflinks.txt | 80 ++++++++++++ refTrans/m-data-analysis.txt | 113 ++-------------- refTrans/m-dataNsoftware_HPC.txt | 14 +- refTrans/m-edgeR.txt | 138 ++++++++++++++++++++ refTrans/m-edgeR_difExp.txt | 63 +++++++++ refTrans/m-func-analysis.txt | 50 +++++++ refTrans/m-hpc-qsub.txt | 41 ------ refTrans/m-htseq.txt | 81 +++++++++--- refTrans/m-more-tophat.txt | 67 ---------- refTrans/m-quality.txt | 17 ++- refTrans/m-tophat.txt | 18 +-- refTrans/model-orgs.txt | 8 +- refTrans/more-resources.txt | 18 --- refTrans/semi-orgs.txt | 10 +- 20 files changed, 650 insertions(+), 324 deletions(-) create mode 100644 refTrans/m-HPC-cuffdiff.txt rename refTrans/{para_HPC.txt => m-HPC-tophat.txt} (79%) create mode 100644 refTrans/m-count.txt create mode 100644 refTrans/m-cuffdiff.txt create mode 100644 refTrans/m-cufflinks.txt create mode 100644 refTrans/m-edgeR.txt create mode 100644 refTrans/m-edgeR_difExp.txt create mode 100644 refTrans/m-func-analysis.txt delete mode 100644 refTrans/m-hpc-qsub.txt delete mode 100644 refTrans/m-more-tophat.txt delete mode 100644 refTrans/more-resources.txt diff --git a/refTrans/index.txt b/refTrans/index.txt index a203b42..61faa8d 100644 --- a/refTrans/index.txt +++ b/refTrans/index.txt @@ -2,10 +2,10 @@ Reference based RNAseq assembly (refTrans) protocol =================================================== -:author: Tamer Mansour and Titus Brown +:author: Tamer A. Mansour and Titus Brown 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. +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 `__ and made use of the `MSU High Performance Compute Center `__ . @@ -18,4 +18,5 @@ and made use of the `MSU High Performance Compute Center ` Reference material ------------------ :doc:`../docs/command-line` + .. :doc:`../amazon/index` \ No newline at end of file diff --git a/refTrans/m-HPC-cuffdiff.txt b/refTrans/m-HPC-cuffdiff.txt new file mode 100644 index 0000000..81f25ec --- /dev/null +++ b/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 ' to check on your jobs' status. + +.. shell stop + +---- + +Back: :doc:`m-cuffdiff` diff --git a/refTrans/para_HPC.txt b/refTrans/m-HPC-tophat.txt similarity index 79% rename from refTrans/para_HPC.txt rename to refTrans/m-HPC-tophat.txt index 7841b24..956bb59 100644 --- a/refTrans/para_HPC.txt +++ b/refTrans/m-HPC-tophat.txt @@ -1,5 +1,7 @@ -Submitting jobs to the MSU HPC queue -==================================== +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 @@ -14,24 +16,29 @@ 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:: +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 "_") + output=$(basename "$f" | cut -f 1,2 -d "_")_tophat script="./qsub_scripts/"$output"_script.sh" - printf $input_R1"\n"$input_R2"\n"$input_R"\n"$output"\n" 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 ' to check on your jobs' status. +.. shell stop + +---- + Back: :doc:`m-tophat` diff --git a/refTrans/m-advice.txt b/refTrans/m-advice.txt index 3bf6f32..4c22420 100644 --- a/refTrans/m-advice.txt +++ b/refTrans/m-advice.txt @@ -18,44 +18,6 @@ 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. -Downloading your data ---------------------- - -If you do your sequencing at the MSU Core Facility, you'll get an e-mail -from them when you're samples are ready. The e-mail will give you an -FTP site, a username, and a password, as well as a URL. You can -use these to download your data. For example, if you get:: - - hostname: titan.bch.msu.edu - username: rnaseqmodel - password: QecheJa6 - - URI: ftp://rnaseqmodel:QecheJa6@titan.bch.msu.edu - -you can go to ftp://rnaseqmodel:QecheJa6@titan.bch.msu.edu in your -Web browser; that is, it lets you combine your username and password -to open that link. - -In this case, you will see a 'testdata' directory. If you click on that, -you'll see a bunch of fastq.gz files. These are the files that you -want to get onto the HPC. - -To download these files onto the HPC, log into the HPC, go to the -directory on the HPC you want to put the files in, and run a 'wget' -- -for example, on the HPC:: - - mkdir ~/testdata - cd ~/testdata - - wget -r -np -nH ftp://rnaseqmodel:QecheJa6@titan.bch.msu.edu/testdata/ - -This will download _all_ of the files in that directory. You can also do -them one at a time, e.g. to get 'Ath_Mut_1_R1.fastq.gz', you would do - - wget ftp://rnaseqmodel:QecheJa6@titan.bch.msu.edu/testdata/Ath_Mut_1_R1.fastq.gz - -Tada! - Developing your own pipeline ---------------------------- @@ -90,4 +52,20 @@ you'll need to develop your own analysis pipeline. Here are some tips. 8. Provide your scripts and raw counts files as part of any publication or thesis, perhaps via `figshare `__. -Next: :doc:`more-resources` + +Informational resources +---------------------- + +`UT (Austin) Sequencing Core prices `__ - costs and yields for sequencing. + +`ANGUS - summer NGS course `__ - lots of resources and materials and book reference + +`Data Carpentry `__ - intro to R, etc. + +`Software Carpentry `__ - more scripting, Python, etc. + + +Places to share data, scripts, and results files +------------------------------------------------ + +`Figshare `__. diff --git a/refTrans/m-count.txt b/refTrans/m-count.txt new file mode 100644 index 0000000..1fcd9dc --- /dev/null +++ b/refTrans/m-count.txt @@ -0,0 +1,58 @@ +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 `__. +and `Featurecount `__. +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 `__, +`UCSC `__, and +`Ensemble `__). +Also they are avaliable through other genomic software projects e.g. +`Illumina iGenomes project `__. + +---- + +:doc:`./m-cufflinks` + +Mapping the reads to transcripts (gene isoforms) rather than the whole gene will enable +the identification of changes in alterative splicaing machinary. +`Cufflinks `__ is the canonical software +package to do that. Cufflinks can do the mapping based on avaliable annoataion files +but also it can make its own gene model. This feature enables Cufflinks to predict denovo transcripts. +On the other hand, mapping to transcripts requires the implementation of error prone modeling +approaches and thus requires rigorous laboratory confirmation. 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. +There are other softwares for analysis of differential splicing events (exon skipping, intron retention, etc.) +e.g. `MATS `__, and +`MISO `__. + + +---- + +Links: + +* `Trapnell et al. 2013 Differential analysis of gene regulation at transcript resolution with RNA-seq `__ + +---- + +Next: :doc:`m-m-htseq` diff --git a/refTrans/m-cuffdiff.txt b/refTrans/m-cuffdiff.txt new file mode 100644 index 0000000..a5a7dcb --- /dev/null +++ b/refTrans/m-cuffdiff.txt @@ -0,0 +1,56 @@ +Differential expression by Cuffdiff +=================================== + +.. shell start + +We'll be using `Cufdiff +`__ +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 ' '` diff --git a/refTrans/m-cufflinks.txt b/refTrans/m-cufflinks.txt new file mode 100644 index 0000000..f195149 --- /dev/null +++ b/refTrans/m-cufflinks.txt @@ -0,0 +1,80 @@ +Counting the reads that map to each transcript +============================================== + +.. shell start + +Now that we know which reads go with which gene, we'll use +`Cufflinks `__ to allocate those reads +to the possible transcripts of every gene + +.. note:: This section reuires cufflinks/2.1.1 installed + +You can ask Cufflinks to use only annotated transcripts by the argument (--GTF). +However there is another argument (--GTF-guide) which Tells Cufflinks to use the +supplied reference annotation to guide reference annotation based transcript assembly. +Output will include all reference transcripts as well as any novel genes and isoforms that are assembled. + +As we discussed before, Cufflinks uses the ambiguously aligned reads in its calculation. +To do this in a more conservative approach, there are several arguments that you can pass to Cufflinks: +1. --multi-read-correct Tells Cufflinks to do an initial estimation procedure to more accurately weight reads mapping to multiple locations in the genome +2. --frag-bias-correct Providing Cufflinks with a genome reference to run a bias detection and correction algorithm which improve accuracy of transcript abundance estimates. + +:: + + for f in *_tophat; do + cufflinks \ + --GTF-guide "$GTF_file" \ + --multi-read-correct \ + --frag-bias-correct ${homo_sapiens_genome} \ + --num-threads 4 \ + --output-dir $f/cufflink_out \ + $f/accepted_hits.bam + done + +Out folder will have four files:: + + transcripts.gtf : This GTF file contains Cufflinks’ assembled isoforms. + isoforms.fpkm_tracking : This file contains the estimated isoform-level expression values + genes.fpkm_tracking : This file contains the estimated gene-level expression values + skipped.gtf : List of skipped loci. Cufflinks skip a locus if more than 1000000 fragment + mapped to it. You can change this through the (–max-bundle-frags) + +Note, Cufflinks express the expression values in **FPKM** (Fragments per kilobase of exon per million reads mapped). +This unit corrects for the gene size and the depth of sequencing. So let us say that the expression of a given gene was 100, +this means that for every one million fragments(read if single ended or pair if paired ended) +mapping to the whole transcriptome, 100 fragments of this million will map to one kilobase of this transcript. + +---- + +Because we have multiple RNA-Seq libraries, It is recommend that you merge these assemblies into a master transcriptome. +This step is required for a differential expression analysis of the new transcripts we’ve assembled. +`Cuffmerge `__ performs this merge step. + +First step is to specify all the input GTF files in a **manifest** file listing full paths to the files. +:: + + for f in *_tophat; do + echo $workingPath/$f/cufflink_out/transcripts.gtf >> assembly_GTF_list.txt + done + +Now we can run the cuffmerge software +:: + + cuffmerge \ + --ref-sequence ${homo_sapiens_genome} \ + --ref-gtf $GTF_file \ + --num-threads 4 \ + assembly_GTF_list.txt + +The cuffmerge will output a folder called "merged_asm" which contain "merged.gtf" file +containing all the transcripts that cufflinks had assembled. Cuffmrege has an optional argument +(--ref-gtf) which allow you to pass the reference annotation. Therefore, the merged.gtf has cufflinks +transcripts with their matching reference transcripts. + + +.. shell stop + +---- + +Next: :doc:`m-Cuffdif` + diff --git a/refTrans/m-data-analysis.txt b/refTrans/m-data-analysis.txt index c249c47..60773f0 100644 --- a/refTrans/m-data-analysis.txt +++ b/refTrans/m-data-analysis.txt @@ -1,108 +1,17 @@ -Data analysis & differential expression -======================================= +Differential expression +======================= -We'll be using `edgeR -`__ -to do the basic differential expression analysis of our counts. +Now we have tried two different methods to summarize our expression: +* The HTseq has summarized the raw reads per genes and now this ready for further analysis by Packages like DEseq and edgeR. +* Cufflinks has summarized the expression values as FPKM and now this is ready for further analysis by software like Cuffdif. -To run edgeR, you need to write a data loading and manipulation script -in R. In this case, I've provided one -- `lung_saliva.R -`__. -This script will load in two samples with two replicates, execute an -MA plot, and provide a spreadsheet with differential expression -information in it. To download it, `click here -`__. +So now we are going to disscuss the differential expression analysis by these softwares: -Links: +.. toctree:: + :maxdepth: 2 -* `False Discovery Rate `__ -* `Learn R with Swirl `__ -* `Data Carpentry `__ + m-edgeR + m-cuffdiff -Running edgeR on a data subset ------------------------------- -To run the script on the HPC, download the script and data files for -the 100k read subsets. At the command line, do:: - - cd ~/rnaseq - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/lung_saliva.R - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/lung_repl1_counts.txt - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/lung_repl2_counts.txt - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/salivary_repl1_counts.txt - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/salivary_repl2_counts.txt - -Note: the ``saliva_repl1_counts.txt`` and ``lung_repl1_counts.txt`` -are the files produced by :doc:`m-tophat` and :doc:`m-more-tophat`, -respectively. I've also run them on the replicate data sets, which -produced the *repl2_counts.txt files. - -Next, to run the R script, do:: - - module load R - Rscript lung_saliva.R - -This will produce two files, `edgeR-MA-plot.pdf -`__ -and `edgeR-lung-vs-salivary.csv -`__; -they will be in your ``rnaseq`` folder in your home directory -on the HPC. The CSV file can be opened directly in Excel; you can -also look at it `here -`__. -It consists of five columns: gene name, log fold change, P-value, and -FDR-adjusted P-value. - -Links: - -* `edgeR tutorial from UT Austin `__ - -Functional and network analysis on differentially expressed genes ------------------------------------------------------------------ - -There are a number of sites that let you analyze gene lists, including -`DAVID `__. I have gotten a number of -recommendations for `PANTHER `__. PANTHER -lets you upload gene lists and explore their functional categories -interactively or in a gene list/annotation format. - -More specifically, you can explore your RNAseq data with - -* functional classifications, in pie chart or in list; -* tests for statistical overrepresentation; -* tests for statistical enrichment based on associated fold change. - -You need to crop and transform the data a little bit before using the -functional classification. The steps are: - -1. Download `the CSV file `__. If you've produced your own, copy it over from the HPC. - -2. Open it in Excel. - -3. Choose an FDR cutoff (suggest FDR < 0.05 or lower) and delete all the rows after that (or, copy the rows into a new spreadsheet -- might be quicker). - -4. Save as a "Tab-delimited text." (Note, on Mac OS X, you may need to save this as "Windows formatted text" instead.) - -This is now a file you can upload to PANTHER. - -(`You can download my copy of this here `__) - -To do the analysis, go to http://www.pantherdb.org/. In box 1, select "Choose file" -and find the CSV file you want to upload to PANTHER. Nothing else in box -one should be changed. - -In box 2, select "Homo sapiens." - -In box 3, select either "Functional analysis classification viewed in -gene list" or "Functional analysis classification viewed in pie -chart." - -Click submit. - -Now you can explore the results! - -Links: - -* `PANTHER database `__ - -Next: :doc:`m-hpc-qsub` +---- \ No newline at end of file diff --git a/refTrans/m-dataNsoftware_HPC.txt b/refTrans/m-dataNsoftware_HPC.txt index ebff038..7f6dee2 100644 --- a/refTrans/m-dataNsoftware_HPC.txt +++ b/refTrans/m-dataNsoftware_HPC.txt @@ -115,7 +115,6 @@ and do the same for the lung samples gunzip -c ERR315424_1.fastq.gz | head -400000 | gzip > dataSubset/lung_repl2_R1.fq.gz gunzip -c ERR315424_2.fastq.gz | head -400000 | gzip > dataSubset/lung_repl2_R2.fq.gz - .. :: fi @@ -126,6 +125,12 @@ Rather than copying the files into the working directory, let this creates a reference so that UNIX knows where to find them but doesn’t need to actually move them around. :: ln -fs $workingPath/refTransData/dataSubset/*.fq.gz $workingPath + +.. note:: + + You can use your own data instead. The code can be applied to any paired end data set whatever the number of the studied + groups and whatever the number of replicates. To minimize the code manipulation, you would rename your samples according to match + our conventions: _repl#_R(1or2).fq.qz install/load software --------------------- @@ -139,8 +144,13 @@ Now let us load all the modules reuired to go through this protocol module load FastQC/0.11.2 module load Trimmomatic/0.32 module load TopHat2/2.0.12 + module load PySAM/0.6 + module load HTSeq/0.6.1 + module load R + module load cufflinks/2.1.1 + -.. note:: Loading TopHat2 include loading Bowtie/2.2.3.0 and Samtools/0.1.19.0 +.. note:: Loading TopHat2 on MSU HPC includes loading Bowtie/2.2.3.0 and Samtools/0.1.19.0 We will elaborate more on each software when we use it diff --git a/refTrans/m-edgeR.txt b/refTrans/m-edgeR.txt new file mode 100644 index 0000000..5f80db9 --- /dev/null +++ b/refTrans/m-edgeR.txt @@ -0,0 +1,138 @@ +Differential expression by EdgeR +================================ + +.. shell start + +We'll be using `edgeR +`__ +to do the basic differential expression analysis of our counts. + +Create R script to run edgeR +---------------------------- + +To run edgeR, you need to write a data loading and manipulation script in R. +This script will load in all the studied samples, execute an MA plot, +and provide a spreadsheet with differential expression of the studied tissues. +We will go togther step by step to understand and write this scrip through the command line + +First, We need to combine all the count files of the studied samples into one dataset. +Also we need to classfy these samples into groups to compare them. In this case we will have +two replicates for each studied tissue. This code this going to deals with your +samples whatever the no of replicates and the studied groups. + +:: + + echo " + First_sample <- TRUE + columnID <- 2 + group <- c() + for (tophat_folder in dir(pattern='*_tophat')){ + file.name <- paste(tophat_folder, '/total_count.txt', sep='') + sample.name <- sub('_tophat','',tophat_folder) + if (First_sample){ + all.data <- read.delim(file.name, col.names=c('gene', sample.name), sep='\t', colClasses=c('character', 'numeric'), row.names=1) + First_sample <- FALSE + } else { + temp <- read.delim(file.name, col.names=c('gene', sample.name), sep='\t', colClasses=c('character', 'numeric'), row.names=1) + #all.data$sample.name <- temp$sample.name + all.data[,columnID] <- temp[,1] + colnames(all.data)[columnID] = sample.name + columnID = columnID + 1 + } + group_label <- sub('_repl[1-2]','', sample.name) + group <- c(group,group_label) + }" > edgeR_difExp.R + +The count output of HTseqcount ends with five lines of summary statistics about the mapped reads. +So we need to the eliminate the last five rows in our dataset. + +:: + + echo "all.data <- all.data[1:(nrow(all.data)-5),]" >> edgeR_difExp.R + +We now have a data frame with all the data in it. let us have an additional +line to show us the first few lines of our dataset:: + + echo "head(all.data)" >> edgeR_difExp.R + +It is time to load edgR package into memory o start working with + +:: + + echo "library('edgeR')" >> edgeR_difExp.R + +Now, it is time to calculte the differention expression. +`edgeR `__ +has different functions for calculation of differential expression. +We are using the 2 group comparison model. If you have more than 2 groups to compare, +you will have to select the appropriate functions according to the edgR manual + +:: + + echo " + dge = DGEList(counts=all.data, group=group) + dge <- estimateCommonDisp(dge) + dge <- estimateTagwiseDisp(dge) + et <- exactTest(dge) + etp <- topTags(et, n=100000)" >> edgeR_difExp.R + +Let's add a line to have a quick view of the results :: + + echo "summary(etp\$table)" >> edgeR_difExp.R + +We also can plot our results + +:: + + echo " + etp\$table\$logFC = -etp\$table\$logFC + pdf('edgeR-MA-plot.pdf') + plot( + etp\$table\$logCPM, + etp\$table\$logFC, + xlim=c(-3, 20), ylim=c(-12, 12), pch=20, cex=.3, + col = ifelse( etp\$table\$FDR < .1, 'red', 'black' ) ) + dev.off()" >> edgeR_difExp.R + +Finally we will write a spreadsheet with differential expression of the studied tissues + +:: + + echo "write.csv(etp\$table, 'edgeR-lung-vs-salivary.csv')" >> edgeR_difExp.R + + +.. Note:: + + you can copy the script from this link and save it in the working directory + under the name edgeR_difExp.R + + :doc:`./m-edgeR_difExp` + + +Running edgeR on a data subset +------------------------------ + +Next, to run the R script, do + +:: + + Rscript edgeR_difExp.R + +This will produce two files, `edgeR-MA-plot.pdf +`__ +and `edgeR-lung-vs-salivary.csv +`__; +The last file consists of five columns: gene name, log fold change, P-value, and +FDR-adjusted P-value. + +.. shell stop + +---- + +Links: + +* `edgeR tutorial from UT Austin `__ + +---- + +Next: :doc:`m-func-analysis` diff --git a/refTrans/m-edgeR_difExp.txt b/refTrans/m-edgeR_difExp.txt new file mode 100644 index 0000000..b6d7caf --- /dev/null +++ b/refTrans/m-edgeR_difExp.txt @@ -0,0 +1,63 @@ +edgeR_difExp.R +============== +You can copy the script from here and save it in the working directory under the name edgeR_difExp.R:: + + # First, combine all the count files into one dataset + # and create a list of group labels where all the replicates from the same group have the same label + First_sample <- TRUE + columnID <- 2 + group <- c() + for (tophat_folder in dir(pattern="*_tophat")){ + file.name <- paste(tophat_folder, "/total_count.txt", sep="") + sample.name <- sub("_tophat","",tophat_folder) + if (First_sample){ + all.data <- read.delim(file.name, col.names=c("gene", sample.name), sep="\t", colClasses=c("character", "numeric"), row.names=1) + First_sample <- FALSE + } else { + temp <- read.delim(file.name, col.names=c("gene", sample.name), sep="\t", colClasses=c("character", "numeric"), row.names=1) + #all.data$sample.name <- temp$sample.name + all.data[,columnID] <- temp[,1] + colnames(all.data)[columnID] = sample.name + columnID = columnID + 1 + } + group_label <- sub("_repl[1-2]","", sample.name) + group <- c(group,group_label) + } + + # eliminate the last five rows + all.data <- all.data[1:(nrow(all.data)-5),] + + # We now have a data frame with all the data in it! + head(all.data) + + ### + + # now, do differential expression w/edgeR + library("edgeR") + + # calculate differential expression and dispersion + dge = DGEList(counts=all.data, group=group) + dge <- estimateCommonDisp(dge) + dge <- estimateTagwiseDisp(dge) + + # plot! + et <- exactTest(dge) + etp <- topTags(et, n=100000) + + summary(etp$table) + + etp$table$logFC = -etp$table$logFC + pdf("edgeR-MA-plot.pdf") + plot( + etp$table$logCPM, + etp$table$logFC, + xlim=c(-3, 20), ylim=c(-12, 12), pch=20, cex=.3, + col = ifelse( etp$table$FDR < .1, "red", "black" ) ) + dev.off() + + # output CSV + write.csv(etp$table, "edgeR-lung-vs-salivary.csv") + +---- + +Back: :doc:`m-data-analysis` diff --git a/refTrans/m-func-analysis.txt b/refTrans/m-func-analysis.txt new file mode 100644 index 0000000..9da065a --- /dev/null +++ b/refTrans/m-func-analysis.txt @@ -0,0 +1,50 @@ +Functional and network analysis on differentially expressed genes +================================================================= + +There are a number of sites that let you analyze gene lists, including +`DAVID `__. I have gotten a number of +recommendations for `PANTHER `__. PANTHER +lets you upload gene lists and explore their functional categories +interactively or in a gene list/annotation format. + +More specifically, you can explore your RNAseq data with + +* functional classifications, in pie chart or in list; +* tests for statistical overrepresentation; +* tests for statistical enrichment based on associated fold change. + +You need to crop and transform the data a little bit before using the +functional classification. The steps are: + +1. Download `the CSV file `__. If you've produced your own, copy it over from the HPC. + +2. Open it in Excel. + +3. Choose an FDR cutoff (suggest FDR < 0.05 or lower) and delete all the rows after that (or, copy the rows into a new spreadsheet -- might be quicker). + +4. Save as a "Tab-delimited text." (Note, on Mac OS X, you may need to save this as "Windows formatted text" instead.) + +This is now a file you can upload to PANTHER. + +(`You can download my copy of this here `__) + +To do the analysis, go to http://www.pantherdb.org/. In box 1, select "Choose file" +and find the CSV file you want to upload to PANTHER. Nothing else in box +one should be changed. + +In box 2, select "Homo sapiens." + +In box 3, select either "Functional analysis classification viewed in +gene list" or "Functional analysis classification viewed in pie +chart." + +Click submit. + +Now you can explore the results! + +---- + +Links: + +* `PANTHER database `__ + diff --git a/refTrans/m-hpc-qsub.txt b/refTrans/m-hpc-qsub.txt deleted file mode 100644 index 4fd0d51..0000000 --- a/refTrans/m-hpc-qsub.txt +++ /dev/null @@ -1,41 +0,0 @@ -Submitting jobs to the MSU HPC queue -==================================== - -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'. - -Here are four scripts (plus a common library script) that you could use -to run some analyses on the HPC. To run, do something like this:: - - module load powertools - getexample RNAseq-model - - mkdir ~/rnaseq - mkdir ~/rnaseq/script - cd ~/rnaseq/script - ln -fs ~/RNAseq-model/data/*.fastq . - - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/env.sh - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/process-1.sh - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/process-2.sh - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/process-3.sh - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/subset/process-4.sh - - qsub process-1.sh - qsub process-2.sh - qsub process-3.sh - qsub process-4.sh - -You can use 'qstat | grep ' to check on your jobs' status. - -Next: :doc:`m-advice` diff --git a/refTrans/m-htseq.txt b/refTrans/m-htseq.txt index cc2f732..2508756 100644 --- a/refTrans/m-htseq.txt +++ b/refTrans/m-htseq.txt @@ -1,25 +1,72 @@ Counting the reads that map to each gene ======================================== +.. shell start + Now that we know which reads go with which gene, we'll use `htseq-count `__. -First, load the PySAM and HTSeq software packages:: +.. note:: This section reuires PySAM/0.6 and HTSeq/0.6.1 installed + +HTSeq can not deal with BAM files containing mix of paired end and single end reads. +So we need to split each BAM file according to the type of reads. `Samtools +`__ can do this job for us using +the view command. Paire end reads are labeled with FLAG=1 so we can separate the +paired reads using the (-f) option which will *keep* only reads labeled with FALG=1. +Also we can get the single end reads in a smilliar way but using (-F) option +which will *exclude* reads labeled with FALG=1. + +Let us do this for all the BAM files +:: + + for f in *_tophat; do + printf $f"\n" + samtools view -f 1 -bh -o "$f"/accepted_hits_PE.bam "$f"/accepted_hits.bam + samtools view -F 1 -bh -o "$f"/accepted_hits_SE.bam "$f"/accepted_hits.bam + done - module load PySAM/0.6 - module load HTSeq/0.6.1 +Second preparatory step is to sort the paired end BAM files by name. You do not have to +but HTSeq might through an error if you are using a huge BAM file. +:: -And next, run HTSeq:: + for f in *_tophat; do + printf $f"\n" + samtools sort -n $f/accepted_hits_PE.bam $f/accepted_hits_PE2 + done - htseq-count --format=bam --stranded=no --order=pos \ - tophat_salivary_repl1/accepted_hits.bam \ - ~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Annotation/Genes/genes.gtf > salivary_repl1_counts.txt +And next, run HTSeq +:: + + for f in *_tophat; do + output_PE=$(basename "$f" | cut -f 1,2 -d "_")_PEcounts.txt + printf $output_PE"\n" + + htseq-count --format=bam --stranded=no --order=name \ + $f/accepted_hits_PE2.bam \ + $GTF_file > $f/$output_PE + + output_SE=$(basename "$f" | cut -f 1,2 -d "_")_SEcounts.txt + printf $output_SE"\n" + htseq-count --format=bam --stranded=no --order=name \ + $f/accepted_hits_SE.bam \ + $GTF_file > $f/$output_SE + done -When this is done, type:: +When this is done, let us view one of the output file:: - less salivary_repl1_counts.txt + less salivary_repl1_tophat/salivary_repl1_counts.txt -(again, use 'q' to exit). These are your gene counts. +(use 'q' to exit). These are your gene counts. + +Final step will be summing the reads from both PE and SE analysis +:: + + for f in *_tophat; do + paste $f/*counts.txt | cut -f1,2,4 > $f/combined_temp.txt + cat $f/combined_temp.txt | awk '{print $2+$3}' > $f/total_temp.txt + paste $f/combined_temp.txt $f/total_temp.txt | cut -f1,4 > $f/total_count.txt + rm $f/*_temp.txt + done Note, these are *raw* gene counts - the number of reads that map to each feature (gene, in this case). They are not normalized by length @@ -27,16 +74,8 @@ of gene. According to `this post on seqanswers `__, both DEseq and edgeR want exactly this kind of information! -Questions: - -* what are the 'ENS*...' names? -* what do these parameters mean? -* what parameters does HTSeq take? -* why are we using so many programs? - -Links: +.. shell stop -* `HTSeq `__ -* `htseq-count documentation `__ +---- -Next: :doc:`m-more-tophat` +Next: :doc:`m-data-analysis` diff --git a/refTrans/m-more-tophat.txt b/refTrans/m-more-tophat.txt deleted file mode 100644 index ea518af..0000000 --- a/refTrans/m-more-tophat.txt +++ /dev/null @@ -1,67 +0,0 @@ -Processing another sample with TopHat and HTSeq -=============================================== - -Let's script this! - -Type:: - - cd ~/rnaseq/ - module load Trimmomatic/0.32 - cat > lung-tophat.sh < lung_repl1_R1.fq.gz - head -400000 ~/RNAseq-model/data/ERR315326_2.fastq | gzip > lung_repl1_R2.fq.gz - - # run Trimmomatic. Here, the inputs are 'lung_repl1_R1.fq.gz' and - # 'lung_repl1_R2.fq.gz', and the outputs are 'lung_repl1_R1.qc.fq.gz' - # and 'lung_repl1_R2.qc.fq.gz'. - java -jar \$TRIM/trimmomatic PE lung_repl1_R1.fq.gz lung_repl1_R2.fq.gz \ - lung_repl1_R1.qc.fq.gz s1_se lung_repl1_R2.qc.fq.gz s2_se \ - ILLUMINACLIP:\$TRIM/adapters/TruSeq2-PE.fa:2:40:15 \ - LEADING:2 TRAILING:2 \ - SLIDINGWINDOW:4:2 \ - MINLEN:25 - - # now run Tophat! - # The inputs are the outputs of the previous Trimmomatic step. - # The outputs are going to be under the 'tophat_lung_repl1' directory. - tophat -p 4 \ - -G ~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Annotation/Genes/genes.gtf \ - --transcriptome-index=$HOME/RNAseq-model/transcriptome \ - -o tophat_lung_repl1 \ - ~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Sequence/Bowtie2Index/genome \ - lung_repl1_R1.qc.fq.gz lung_repl1_R2.qc.fq.gz - - # count the hits by gene -- 'tophat_lung_repl1' is the main output, - # from Tophat. - htseq-count --format=bam --stranded=no --order=pos tophat_lung_repl1/accepted_hits.bam \ - ~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Annotation/Genes/genes.gtf > lung_repl1_counts.txt - EOF - -(Be sure to press the Enter or Return key after pasting this in!) This is -`called a 'heredoc' `__ -and it gives a way to write a shell script via copy-paste ;). - -Next, type:: - - bash lung-tophat.sh - -This will run all of the commands in the file 'lung-tophat.sh'. - -You can use the 'nano' editor to modify this file -- type:: - - nano lung-tophat.sh - -Next: :doc:`m-data-analysis` diff --git a/refTrans/m-quality.txt b/refTrans/m-quality.txt index 6652268..1ac24e5 100644 --- a/refTrans/m-quality.txt +++ b/refTrans/m-quality.txt @@ -7,7 +7,9 @@ Assessment of sequence quality (FastQC) --------------------------------------- We're going to use `FastQC `__ to summarize the data. -FastQC cans on Fastq files as well as GZip compressed Fastq (.gz) files. The program operates on one file at a time +FastQC operates on Fastq files as well as GZip compressed Fastq (.gz) files. + +.. note:: This section reuires FastQC/0.11.2 installed Now run FastQC on your sequencing file :: @@ -15,10 +17,7 @@ Now run FastQC on your sequencing file mkdir rawData_fastqc fastqc *.fq.gz --outdir=./rawData_fastqc -.. note:: When you pass multiple files to FastQC, it operates a queueing system -where only one file is opened at a time, and rest of the files will wait until existing -files have been processed. If your files are big in size, it might be more wise -to run then in parallel. +.. .. note:: When you pass multiple files to FastQC, it operates one file at a time give it a seprate output. Now type 'ls' in the output directory to explore the result files:: @@ -49,6 +48,8 @@ We'll be using `Trimmomatic `__. For a discussion of optimal RNAseq trimming strategies, see `MacManes, 2014 `__. +.. note:: This section reuires Trimmomatic/0.32 installed + To run Trimmomatic :: @@ -65,9 +66,11 @@ To run Trimmomatic LEADING:2 TRAILING:2 \ SLIDINGWINDOW:4:15 \ MINLEN:25 &>> log.txt - # Titus workshop: ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 LEADING:2 TRAILING:2 SLIDINGWINDOW:4:2 - # my parameters: ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:30:10:1 SLIDINGWINDOW:4:15 done + +.. # Titus workshop: ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 LEADING:2 TRAILING:2 SLIDINGWINDOW:4:2 +.. # my parameters: ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:30:10:1 SLIDINGWINDOW:4:15 + Now type 'ls' in the output directory to explore the result files:: diff --git a/refTrans/m-tophat.txt b/refTrans/m-tophat.txt index a05b714..d5d0eb0 100644 --- a/refTrans/m-tophat.txt +++ b/refTrans/m-tophat.txt @@ -8,6 +8,8 @@ reads to the reference gene set, for the purpose of counting how many reads have come from each gene. We'll be using the `TopHat software `__ +.. note:: This section reuires TopHat2/2.0.12, Bowtie/2.2.3.0 and Samtools/0.1.19.0 installed + Download the Bowtie2 index -------------------------- @@ -23,7 +25,7 @@ The `Illumina iGenomes project `__ provided the RNA-Seq user community with a set of genome sequence indexes (including Bowtie indexes). We will download the iGenome index for the UCSC hg19 assembly version and -define the path of the genome and genome_index_base which the name of any +define the path of the genome and genome_index_base which is the name of any of the index files up to but not including the first period. .. check if the data directory exists @@ -79,7 +81,6 @@ So we need to merge every two SE files belonging to the same sample input_R1="$f" input_R2=$(echo "$f" | sed s/R1/R2/) output=$(echo "$f" | sed s/R1/R/) - printf $input_R1" "$input_R2" "$output"\n" cat $input_R1 $input_R2 > $output done @@ -92,7 +93,7 @@ And now run TopHat on one sample:: --num-threads 4 \ --GTF "$GTF_file" \ --transcriptome-index "$workingPath/trans_Index/transcriptome" \ - --output-dir salivary_repl1 \ + --output-dir salivary_repl1_tophat \ $Bowtie2Index \ salivary_repl1_R1.pe.fq.gz salivary_repl1_R2.pe.fq.gz,salivary_repl1_R.se.fq.gz @@ -103,8 +104,7 @@ This will take about ~10-15 minutes. We can run all the samples sequencially usi 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 "_") - printf $input_R1"\n"$input_R2"\n"$input_R"\n"$output"\n" + output=$(basename "$f" | cut -f 1,2 -d "_")_tophat tophat \ --num-threads 4 \ --GTF "$GTF_file" \ @@ -114,10 +114,12 @@ This will take about ~10-15 minutes. We can run all the samples sequencially usi $input_R1 $input_R2,$input_R done -With large input files, TopHat will take hours and may be days. It will be more -efficient to run TopHat on all samples in parallel jobs +.. Note:: -:doc:`./para_HPC + With large input files, TopHat will take hours and may be days. It will be more + efficient to run TopHat on all samples in parallel jobs + + :doc:`./m-HPC-tophat` Counting mapped reads percentage diff --git a/refTrans/model-orgs.txt b/refTrans/model-orgs.txt index 68204f8..8d5e11c 100644 --- a/refTrans/model-orgs.txt +++ b/refTrans/model-orgs.txt @@ -6,13 +6,13 @@ Reference based RNAseq assembly in model organisms Tutorials: .. toctree:: - :maxdepth: 1 + :maxdepth: 2 m-dataNsoftware_HPC m-quality m-tophat - m-htseq - m-more-tophat + m-count m-data-analysis + m-func-analysis m-advice -.. more-resources + diff --git a/refTrans/more-resources.txt b/refTrans/more-resources.txt deleted file mode 100644 index c85627d..0000000 --- a/refTrans/more-resources.txt +++ /dev/null @@ -1,18 +0,0 @@ -More resources -============== - -Informational resources ----------------------- - -`UT (Austin) Sequencing Core prices `__ - costs and yields for sequencing. - -`ANGUS - summer NGS course `__ - lots of resources and materials and book reference - -`Data Carpentry `__ - intro to R, etc. - -`Software Carpentry `__ - more scripting, Python, etc. - -Places to share data, scripts, and results files ------------------------------------------------- - -`Figshare `__. diff --git a/refTrans/semi-orgs.txt b/refTrans/semi-orgs.txt index b27c55f..10d4aa9 100644 --- a/refTrans/semi-orgs.txt +++ b/refTrans/semi-orgs.txt @@ -1,12 +1,9 @@ -2014 / December / mRNAseq on semi-model organisms -================================================= - -This workshop was given on December 10th and 11th, 2014, by C. Titus -Brown and Matt Scholz. +Reference based RNAseq assembly in semi-model organisms +======================================================= .. .. figure:: _static/model-rnaseq-pipeline.png -What are semi-model organisms? `One with bad transcriptomes. `__! +What are semi-model organisms? `One with bad genome. `__! Tutorials: @@ -20,4 +17,3 @@ Tutorials: s-even-more-tophat s-data-analysis m-advice - more-resources From d5cb602504fdfa68dc767564f153e0607e5799e8 Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Sat, 7 Feb 2015 11:13:08 -0600 Subject: [PATCH 08/14] second draft: merging the model and semi-model modules --- refTrans/GETTING-STARTED.txt | 11 -- refTrans/firstCommitment.txt | Bin 132 -> 0 bytes refTrans/index.txt | 13 +- refTrans/m-count.txt | 32 +++- refTrans/m-edgeR.txt | 18 +- refTrans/m-tophat.txt | 8 + refTrans/model-orgs.txt | 18 -- refTrans/s-building-a-reference.txt | 121 ------------- refTrans/s-data-analysis.txt | 60 ------- refTrans/s-even-more-tophat.txt | 65 ------- refTrans/s-more-tophat.txt | 48 ------ refTrans/s-quality.txt | 252 ---------------------------- refTrans/s-tophat-round2.txt | 80 --------- refTrans/semi-orgs.txt | 19 --- 14 files changed, 44 insertions(+), 701 deletions(-) delete mode 100644 refTrans/GETTING-STARTED.txt delete mode 100644 refTrans/firstCommitment.txt delete mode 100644 refTrans/model-orgs.txt delete mode 100644 refTrans/s-building-a-reference.txt delete mode 100644 refTrans/s-data-analysis.txt delete mode 100644 refTrans/s-even-more-tophat.txt delete mode 100644 refTrans/s-more-tophat.txt delete mode 100644 refTrans/s-quality.txt delete mode 100644 refTrans/s-tophat-round2.txt delete mode 100644 refTrans/semi-orgs.txt diff --git a/refTrans/GETTING-STARTED.txt b/refTrans/GETTING-STARTED.txt deleted file mode 100644 index 70c85bc..0000000 --- a/refTrans/GETTING-STARTED.txt +++ /dev/null @@ -1,11 +0,0 @@ -1. Edit conf.py and replace 'labibi' in the 'project' line. - -2. Either replace or eliminate the Google Analytics ID, the disqus name, - and the github information. :: - - html_context = { - "google_analytics_id" : 'UA-36028965-1', - "disqus_shortname" : 'labibi', - "github_base_account" : 'ctb', - "github_project" : 'labibi', - } diff --git a/refTrans/firstCommitment.txt b/refTrans/firstCommitment.txt deleted file mode 100644 index 59d9296100432b9afcc9613e8a4e7ab6afa4a5b1..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 132 zcmXAh!3}^Q6hnJ9aR=sb4i`X#ghasr7nfHba%lR}{@?d0@(NmIl~lFJKt`8}vKn&e oAPIQE$G!hO(R#9tp9<}8qZwNciY)b&6LPe@vFzR4M>K2F2Sx=NAOHXW diff --git a/refTrans/index.txt b/refTrans/index.txt index 61faa8d..a8ab136 100644 --- a/refTrans/index.txt +++ b/refTrans/index.txt @@ -9,11 +9,20 @@ Materials are based on and extending two workshops given at MSU by Titus Brown a Both workshops were sponsored by `iCER `__ and made use of the `MSU High Performance Compute Center `__ . +.. figure:: files/model-rnaseq-pipeline.png + +Tutorials: + .. toctree:: :maxdepth: 2 - semi-orgs - model-orgs + m-dataNsoftware_HPC + m-quality + m-tophat + m-count + m-data-analysis + m-func-analysis + m-advice Reference material ------------------ diff --git a/refTrans/m-count.txt b/refTrans/m-count.txt index 1fcd9dc..95609aa 100644 --- a/refTrans/m-count.txt +++ b/refTrans/m-count.txt @@ -31,20 +31,34 @@ Also they are avaliable through other genomic software projects e.g. :doc:`./m-cufflinks` -Mapping the reads to transcripts (gene isoforms) rather than the whole gene will enable -the identification of changes in alterative splicaing machinary. +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 `__, +`development `__, +`immmunology `__, +`fertility `__, and many others. + `Cufflinks `__ is the canonical software -package to do that. Cufflinks can do the mapping based on avaliable annoataion files -but also it can make its own gene model. This feature enables Cufflinks to predict denovo transcripts. -On the other hand, mapping to transcripts requires the implementation of error prone modeling -approaches and thus requires rigorous laboratory confirmation. Cufflinks also +package to do that. But there are other softwares for analysis of differential +splicing events (exon skipping, intron retention, etc.) +e.g. `MATS `__, and +`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. -There are other softwares for analysis of differential splicing events (exon skipping, intron retention, etc.) -e.g. `MATS `__, and -`MISO `__. + +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. + ---- diff --git a/refTrans/m-edgeR.txt b/refTrans/m-edgeR.txt index 5f80db9..ba3fc24 100644 --- a/refTrans/m-edgeR.txt +++ b/refTrans/m-edgeR.txt @@ -7,9 +7,6 @@ We'll be using `edgeR `__ to do the basic differential expression analysis of our counts. -Create R script to run edgeR ----------------------------- - To run edgeR, you need to write a data loading and manipulation script in R. This script will load in all the studied samples, execute an MA plot, and provide a spreadsheet with differential expression of the studied tissues. @@ -19,7 +16,6 @@ First, We need to combine all the count files of the studied samples into one da Also we need to classfy these samples into groups to compare them. In this case we will have two replicates for each studied tissue. This code this going to deals with your samples whatever the no of replicates and the studied groups. - :: echo " @@ -45,7 +41,6 @@ samples whatever the no of replicates and the studied groups. The count output of HTseqcount ends with five lines of summary statistics about the mapped reads. So we need to the eliminate the last five rows in our dataset. - :: echo "all.data <- all.data[1:(nrow(all.data)-5),]" >> edgeR_difExp.R @@ -56,7 +51,6 @@ line to show us the first few lines of our dataset:: echo "head(all.data)" >> edgeR_difExp.R It is time to load edgR package into memory o start working with - :: echo "library('edgeR')" >> edgeR_difExp.R @@ -66,7 +60,6 @@ Now, it is time to calculte the differention expression. has different functions for calculation of differential expression. We are using the 2 group comparison model. If you have more than 2 groups to compare, you will have to select the appropriate functions according to the edgR manual - :: echo " @@ -81,7 +74,6 @@ Let's add a line to have a quick view of the results :: echo "summary(etp\$table)" >> edgeR_difExp.R We also can plot our results - :: echo " @@ -95,11 +87,9 @@ We also can plot our results dev.off()" >> edgeR_difExp.R Finally we will write a spreadsheet with differential expression of the studied tissues - :: - echo "write.csv(etp\$table, 'edgeR-lung-vs-salivary.csv')" >> edgeR_difExp.R - + echo "write.csv(etp\$table, 'edgeR-lung-vs-salivary.csv')" >> edgeR_difExp.R .. Note:: @@ -109,11 +99,7 @@ Finally we will write a spreadsheet with differential expression of the studied :doc:`./m-edgeR_difExp` -Running edgeR on a data subset ------------------------------- - -Next, to run the R script, do - +Now we can run the R script to do the analysis :: Rscript edgeR_difExp.R diff --git a/refTrans/m-tophat.txt b/refTrans/m-tophat.txt index d5d0eb0..d4be2ff 100644 --- a/refTrans/m-tophat.txt +++ b/refTrans/m-tophat.txt @@ -51,7 +51,15 @@ information is typically shared in GTF files. Illumina iGenomes indexes come with GTF transcript annotation files (Congratualtion, this means you have downloaded them already) Let us define the path of the required GTF in the downloaded iGenome package :: + GTF_file="$workingPath/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf" + +.. Note:: + The human genome is compased of multi fasta entries for chr1-22, M, X,and Y. + genes.gtf has some non canonical chromosomal headers e.g. chr17_ctg5_hap1, chr17_gl000205_random, + chrUn_gl000228, and many more. You can learn more about these header from this + `discussion `__ . These header will cause some software to + through some warning messages. Do not worry about it :) We can pass the annotation file to TopHat through the optional argument (--GTF). TopHat will first extract the transcript sequences and use Bowtie to align reads to this virtual transcriptome first. diff --git a/refTrans/model-orgs.txt b/refTrans/model-orgs.txt deleted file mode 100644 index 8d5e11c..0000000 --- a/refTrans/model-orgs.txt +++ /dev/null @@ -1,18 +0,0 @@ -Reference based RNAseq assembly in model organisms -================================================== - -.. 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 - diff --git a/refTrans/s-building-a-reference.txt b/refTrans/s-building-a-reference.txt deleted file mode 100644 index ae8f9e1..0000000 --- a/refTrans/s-building-a-reference.txt +++ /dev/null @@ -1,121 +0,0 @@ -Building a new reference transcriptome -====================================== - -The distinguishing feature of (what I call) semi-model organisms is that -while they may have a decent genome reference, their transcriptome -annotation is poor. There can be several reasons for this, but generally -it boils down to lack of resources and/or attention -- it takes a *lot* -of effort to build a high quality transcriptome! - -For this purpose, we've already installed the chicken reference genome -set on the HPC (as part of the data you loaded at the beginning). In -this case we've loaded in the `Illumina iGenomes project -`__ into the RNAseq-semimodel -location. - -See paper: - - http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3334321/ - -Map all the reads to the genome with TopHat -------------------------------------------- - -.. @@ add links etc. - -We'll be using the `TopHat software -`__. - -Do:: - - module load TopHat2/2.0.12 - - tophat -p 4 \ - -G $HOME/RNAseq-semimodel/reference/Gallus_gallus/UCSC/galGal3/Annotation/Genes/genes.gtf \ - --transcriptome-index=$HOME/RNAseq-semimodel/tophat/transcriptome \ - -o tophat_all \ - $HOME/RNAseq-semimodel/reference/Gallus_gallus/UCSC/galGal3/Sequence/Bowtie2Index/genome \ - female_repl1_R1.qc.fq.gz,male_repl1_R1.qc.fq.gz,female_repl2_R1.qc.fq.gz,male_repl2_R1.qc.fq.gz \ - female_repl1_R2.qc.fq.gz,male_repl1_R2.qc.fq.gz,female_repl2_R2.qc.fq.gz,male_repl2_R2.qc.fq.gz - -Questions: - -* What are all these parameters?! -* How do we pick the transcriptome/genome? -* Why is it so slow? -* What is different about mapping RNAseq reads vs mapping genomic reads? - -Links: - -* `TopHat manual `__ -* `Illumina iGenomes project `__ - -Evaluating the mapping ----------------------- - -Check out the details:: - - less tophat_all/align_summary.txt - -Build a new transcriptome ("ab initio") from the combined reads using Cufflinks -------------------------------------------------------------------------------- - -Now that we've mapped the reads, let's put them all together into exons -and gene models:: - - module load cufflinks/2.2.1 - - cufflinks -o cuff_all tophat_all/accepted_hits.bam - -.. @@ cufflinks diagram - -Questions: - -* What exactly is Cufflinks doing? - -Merge the new transcriptome with the existing reference transcriptome ---------------------------------------------------------------------- - -We already have some decent gene models; let's merge our new and the old ones:: - - ls -1 cuff_all/transcripts.gtf > cuff_list.txt - - cuffmerge -g $HOME/RNAseq-semimodel/reference/Gallus_gallus/UCSC/galGal3/Annotation/Genes/genes.gtf \ - -o cuffmerge_all \ - -s $HOME/RNAseq-semimodel/reference/Gallus_gallus/UCSC/galGal3/Sequence/WholeGenomeFasta/genome.fa \ - cuff_list.txt - -Do some cleanup:: - - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/remove-nostrand.py - python remove-nostrand.py cuffmerge_all/merged.gtf > cuffmerge_all/nostrand.gtf - -Questions: - -* why do you want to merge? -* why would you have a list of more than one thing in list.txt? -* come to think of it, why aren't you (re)mapping all your reads every time? -* what's with the 'remove nostrand' script? - -Extracting your new transcriptome sequences -------------------------------------------- - -To get a look at the actual DNA sequences, do:: - - gffread -w cuffmerge_all.fa \ - -g $HOME/RNAseq-semimodel/reference/Gallus_gallus/UCSC/galGal3/Sequence/WholeGenomeFasta/genome.fa \ - cuffmerge_all/nostrand.gtf - -Questions: - -* What's the difference between a GTF file and the FA file? - -Checking out your new transcriptome ------------------------------------ - -Take a look at the top of your FASTA file:: - - head -30 cuffmerge_all.fa - -Head on over to `the chicken genome browser `__ and try BLATing the sequence! - -:doc:`s-tophat-round2` diff --git a/refTrans/s-data-analysis.txt b/refTrans/s-data-analysis.txt deleted file mode 100644 index 9928247..0000000 --- a/refTrans/s-data-analysis.txt +++ /dev/null @@ -1,60 +0,0 @@ -Data analysis & differential expression -======================================= - -We'll be using `edgeR -`__ -to do the basic differential expression analysis of our counts. - -To run edgeR, you need to write a data loading and manipulation script -in R. In this case, I've provided one -- `chick.R -`__. -This script will load in two samples with two replicates, execute an -MA plot, do an MDS analysis/plot, and provide a spreadsheet with -differential expression information in it. To download it, `click -here -`__. - -Links: - -* `False Discovery Rate `__ -* `Learn R with Swirl `__ -* `Data Carpentry `__ - -Running edgeR on a data subset ------------------------------- - -To run the script on the HPC, download the script and data files for -the 100k read subsets. At the command line, do:: - - cd ~/rnaseq - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/chick.R - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/chick-subset/female_repl1_counts.txt - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/chick-subset/female_repl2_counts.txt - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/chick-subset/male_repl1_counts.txt - curl -O http://2014-msu-rnaseq.readthedocs.org/en/latest/_static/chick-subset/male_repl2_counts.txt - -These various .txt files are produced by :doc:`s-tophat-round2`, :doc:`s-more-tophat`, and :doc:`s-even-more-tophat`. - -Next, to run the R script, do:: - - module load R - Rscript chick.R - -This will produce three files, `chick-edgeR-MA-plot.pdf -`__ -`chick-edgeR-MDS.pdf -`__ -and `chick-edgeR.csv -`__; -they will be in your ``rnaseq`` folder in your home directory -on the HPC. The CSV file can be opened directly in Excel; you can -also look at it `here -`__. -It consists of five columns: gene name, log fold change, P-value, and -FDR-adjusted P-value. - -Links: - -* `edgeR tutorial from UT Austin `__ - -Next: :doc:`m-advice` diff --git a/refTrans/s-even-more-tophat.txt b/refTrans/s-even-more-tophat.txt deleted file mode 100644 index 7c8a591..0000000 --- a/refTrans/s-even-more-tophat.txt +++ /dev/null @@ -1,65 +0,0 @@ -Using the HPC to run jobs in parallel -===================================== - -Type:: - - cat > chick-tophat-3.sh < male_repl1_counts.txt - - EOF - -Next, do it for male_repl2:: - - cat > chick-tophat-4.sh < male_repl2_counts.txt - - EOF - -To submit these to the queue, do:: - - qsub chick-tophat-3.sh - qsub chick-tophat-4.sh - -Next: doc:`s-data-analysis` diff --git a/refTrans/s-more-tophat.txt b/refTrans/s-more-tophat.txt deleted file mode 100644 index 2aadb26..0000000 --- a/refTrans/s-more-tophat.txt +++ /dev/null @@ -1,48 +0,0 @@ -Processing another sample with TopHat and HTSeq -=============================================== - -Let's script this! - -Type:: - - cd ~/rnaseq/ - cat > chick-tophat-2.sh < female_repl2_counts.txt - - EOF - -(Be sure to press the Enter or Return key after pasting this in!) This is -`called a 'heredoc' `__ -and it gives a way to write a shell script via copy-paste ;). - -Next, type:: - - bash chick-tophat-2.sh - -This will run all of the commands in the file 'chick-tophat-2.sh'. - -You can use the 'nano' editor to modify this file -- type:: - - nano chick-tophat-2.sh - -:doc:`s-even-more-tophat` diff --git a/refTrans/s-quality.txt b/refTrans/s-quality.txt deleted file mode 100644 index 5e9c35d..0000000 --- a/refTrans/s-quality.txt +++ /dev/null @@ -1,252 +0,0 @@ -Short read quality and trimming -=============================== - -.. note:: - - Reminder: if you're on Windows, you should install `mobaxterm `__. - -Log into the HPC with SSH; use your MSU NetID and log into the machine -'hpc.msu.edu'. There copy/paste:: - - cd - module load powertools - getexample RNAseq-semimodel - -This will put all of the example files for today in your home directory -under the directory 'RNAseq-semimodel'. - -0. Getting the data -------------------- - -http://genomebiology.com/content/14/3/R26 - -http://www.ebi.ac.uk/ena/data/view/SRA055442 - -http://www.ebi.ac.uk/ena/data/view/SAMN01096082 -http://www.ebi.ac.uk/ena/data/view/SAMN01096083 - -http://www.ebi.ac.uk/ena/data/view/SAMN01096084 -http://www.ebi.ac.uk/ena/data/view/SAMN01096085 - -Note that each -sample has two replicates, and each replicate has two files. - -**Don't download them**, but if you were downloading these yourself, -you would want the "Fastq files (ftp)", both File 1 and File 2. (They -take a few hours to download!) - -We've already loaded the data onto the MSU HPC, and you've loaded -it with 'module load powertools'. - -To log into a compute node, type:: - - ~/RNAseq-semimodel/login.sh - -If this doesn't work, do:: - - ssh dev-intel14-phi - -Now do:: - - ls -l ~/RNAseq-semimodel/data/ - -You'll see something like :: - - -r--r--r-- 1 mscholz common-data 6781517200 Dec 9 09:46 SRR534005_1.fastq.gz - -r--r--r-- 1 mscholz common-data 7023515467 Dec 9 09:50 SRR534005_2.fastq.gz - -r--r--r-- 1 mscholz common-data 7285848617 Dec 9 09:41 SRR534006_1.fastq.gz - -r--r--r-- 1 mscholz common-data 7542383700 Dec 9 09:43 SRR534006_2.fastq.gz - -r--r--r-- 1 mscholz common-data 7219923066 Dec 9 09:47 SRR536786_1.fastq.gz - -r--r--r-- 1 mscholz common-data 7467116873 Dec 9 09:49 SRR536786_2.fastq.gz - -r--r--r-- 1 mscholz common-data 7694614208 Dec 9 09:40 SRR536787_1.fastq.gz - -r--r--r-- 1 mscholz common-data 7944043814 Dec 9 09:44 SRR536787_2.fastq.gz - -These files are each approximately 7-8 GB in size! - -1. Copying in some data to work with. -------------------------------------- - -First, make a directory:: - - mkdir ~/rnaseq - cd ~/rnaseq - -Copy in a subset of the data (100,000 reads):: - - gunzip -c ~/RNAseq-semimodel/data/SRR534005_1.fastq.gz | head -400000 | gzip > female_repl1_R1.fq.gz - gunzip -c ~/RNAseq-semimodel/data/SRR534005_2.fastq.gz | head -400000 | gzip > female_repl1_R2.fq.gz - -These are FASTQ files -- let's take a look:: - - less female_repl1_R1.fq.gz - -(type 'q' to exit less) - -Question: - -* why are some files named SRR*? -* why are some files named female*? -* why are there R1 and R2 in the name? - -Links: - -* `FASTQ Format `__ - -2. FastQC ---------- - -We're going to use `FastQC `__ to summarize the data. - -First, we need to load the FastQC software into our account:: - - module load FastQC/0.11.2 - -(You have to do this each time you log in and want to use FastQC.) - -Now, run FastQC on both of the female files:: - - fastqc female_repl1_R1.fq.gz - fastqc female_repl1_R2.fq.gz - -Now type 'ls':: - - ls - -and you will see :: - - female_repl1_R1_fastqc.html - female_repl1_R1_fastqc.zip - female_repl1_R2_fastqc.html - female_repl1_R2_fastqc.zip - -Copy these to your laptop and open them in a browser. If you're on a -Mac or Linux machine, you can type:: - - scp username@hpc.msu.edu:rnaseq/female*fastqc.* /tmp - -and then open the html files in your browser. For Windows, if you're using -mobaxterm, most of you should have a file transfer window on the left. -Click 'refresh' (green circle icon fourth from the left) and then navigate -into the 'rnaseq' folder; you should see the 'female_repl...' files there. -Drag and drop those onto your Windows machine. - -You can also view my versions: `female_repl1_R1_fastqc.html -`__ -and `female_repl1_R2_fastqc.html -`__ - -Questions: - -* What should you pay attention to in the FastQC report? -* Which is "better", R1 or R2? - -Links: - -* `FastQC `__ -* `FastQC tutorial video `__ - -3. Trimmomatic --------------- - -Now we're going to do some trimming! We'll be using -`Trimmomatic `__. For -a discussion of optimal RNAseq trimming strategies, see `MacManes, 2014 `__. - -First, load the Trimmomatic software:: - - module load Trimmomatic/0.32 - -Next, run Trimmomatic:: - - java -jar $TRIM/trimmomatic PE female_repl1_R1.fq.gz female_repl1_R2.fq.gz\ - female_repl1_R1.qc.fq.gz s1_se female_repl1_R2.qc.fq.gz s2_se \ - ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 \ - LEADING:2 TRAILING:2 \ - SLIDINGWINDOW:4:2 \ - MINLEN:25 - -You should see output that looks like this:: - - ... - Quality encoding detected as phred33 - Input Read Pairs: 100000 Both Surviving: 95583 (95.58%) Forward Only Surviving: 4262 (4.26%) Reverse Only Surviving: 86 (0.09%) Dropped: 69 (0.07%) - ... - -Questions: - -* How do you figure out what the parameters mean? -* How do you figure out what parameters to use? -* What adapters do you use? -* What version of Trimmomatic are we using here? (And FastQC?) -* Are parameters different for RNAseq and genomic? -* What's with these annoyingly long and complicated filenames? -* What do we do with the single-ended files (s1_se and s2_se?) - -Links: - -* `Trimmomatic `__ - -4. FastQC again ---------------- - -Run FastQC again:: - - fastqc female_repl1_R1.qc.fq.gz - fastqc female_repl1_R2.qc.fq.gz - -(Note that you don't need to load the module again.) - -Copy them to your laptop and open them, OR you can view mine: `female_repl1_R1.qc_fastqc.html -`__ -and `female_repl1_R2.qc_fastqc.html -`__ - -Let's take a look at the output files:: - - less female_repl1_R1.qc.fq.gz - -(again, use 'q' to exit less). - -Questions: - -* Why are some of the reads shorter than others? -* is the quality trimmed data "better" than before? -* Does it matter that you still have adapters!? - -5. Subset and trim the rest of the sequences --------------------------------------------- - -Copy and paste all of the below at once:: - - gunzip -c ~/RNAseq-semimodel/data/SRR534006_1.fastq.gz | head -400000 | gzip > female_repl2_R1.fq.gz - gunzip -c ~/RNAseq-semimodel/data/SRR534006_2.fastq.gz | head -400000 | gzip > female_repl2_R2.fq.gz - - gunzip -c ~/RNAseq-semimodel/data/SRR536786_1.fastq.gz | head -400000 | gzip > male_repl1_R1.fq.gz - gunzip -c ~/RNAseq-semimodel/data/SRR536786_2.fastq.gz | head -400000 | gzip > male_repl1_R2.fq.gz - - gunzip -c ~/RNAseq-semimodel/data/SRR536787_1.fastq.gz | head -400000 | gzip > male_repl2_R1.fq.gz - gunzip -c ~/RNAseq-semimodel/data/SRR536787_2.fastq.gz | head -400000 | gzip > male_repl2_R2.fq.gz - - java -jar $TRIM/trimmomatic PE female_repl2_R1.fq.gz female_repl2_R2.fq.gz\ - female_repl2_R1.qc.fq.gz s1_se female_repl2_R2.qc.fq.gz s2_se \ - ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 \ - LEADING:2 TRAILING:2 \ - SLIDINGWINDOW:4:2 \ - MINLEN:25 - - java -jar $TRIM/trimmomatic PE male_repl1_R1.fq.gz male_repl1_R2.fq.gz\ - male_repl1_R1.qc.fq.gz s1_se male_repl1_R2.qc.fq.gz s2_se \ - ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 \ - LEADING:2 TRAILING:2 \ - SLIDINGWINDOW:4:2 \ - MINLEN:25 - - java -jar $TRIM/trimmomatic PE male_repl2_R1.fq.gz male_repl2_R2.fq.gz\ - male_repl2_R1.qc.fq.gz s1_se male_repl2_R2.qc.fq.gz s2_se \ - ILLUMINACLIP:$TRIM/adapters/TruSeq3-PE.fa:2:40:15 \ - LEADING:2 TRAILING:2 \ - SLIDINGWINDOW:4:2 \ - MINLEN:25 - - -Next: :doc:`s-building-a-reference` diff --git a/refTrans/s-tophat-round2.txt b/refTrans/s-tophat-round2.txt deleted file mode 100644 index 4627dcc..0000000 --- a/refTrans/s-tophat-round2.txt +++ /dev/null @@ -1,80 +0,0 @@ -Mapping reads to the transcriptome with TopHat -============================================== - -Now that we have some quality-controlled reads and a new reference -transcriptome, we're going to map the -reads to the reference genome, **using the new reference transcriptome**. -We'll again be using the `TopHat software -`__ - -Mapping reads -------------- - -Load the TopHat software, if you haven't already:: - - module load TopHat2/2.0.12 - -And now run TopHat:: - - cd ~/rnaseq - tophat -p 4 \ - -G cuffmerge_all/nostrand.gtf \ - --transcriptome-index=$HOME/RNAseq-semimodel/tophat/merged \ - -o tophat_female_repl1 \ - ~/RNAseq-semimodel/reference/Gallus_gallus/UCSC/galGal3/Sequence/Bowtie2Index/genome \ - female_repl1_R1.qc.fq.gz female_repl1_R2.qc.fq.gz - -This will take about 15 minutes. - -Questions: - -* How do we pick the transcriptome/genome? - -Viewing the mapped reads percentage ------------------------------------ - -Let's look at these numbers specifically:: - - less tophat_female_repl1/align_summary.txt - -Making gene counts ------------------- - -Now that we know which reads go with which gene, we'll use -`htseq-count `__. - -First, load the PySAM and HTSeq software packages:: - - module load HTSeq/0.6.1 - -And next, run HTSeq:: - - htseq-count --format=bam --stranded=no --order=pos \ - tophat_female_repl1/accepted_hits.bam \ - cuffmerge_all/nostrand.gtf > female_repl1_counts.txt - -When this is done, type:: - - less female_repl1_counts.txt - -(again, use 'q' to exit). These are your gene counts. - -Note, these are *raw* gene counts - the number of reads that map to -each feature (gene, in this case). They are not normalized by length -of gene. According to `this post on seqanswers -`__, both -DEseq and edgeR want exactly this kind of information! - -Questions: - -* what are the 'TCONS...' names? -* what do these parameters mean? -* what parameters does HTSeq take? -* why are we using so many programs? - -Links: - -* `HTSeq `__ -* `htseq-count documentation `__ - -Next: :doc:`s-more-tophat` diff --git a/refTrans/semi-orgs.txt b/refTrans/semi-orgs.txt deleted file mode 100644 index 10d4aa9..0000000 --- a/refTrans/semi-orgs.txt +++ /dev/null @@ -1,19 +0,0 @@ -Reference based RNAseq assembly in semi-model organisms -======================================================= - -.. .. figure:: _static/model-rnaseq-pipeline.png - -What are semi-model organisms? `One with bad genome. `__! - -Tutorials: - -.. toctree:: - :maxdepth: 2 - - s-quality - s-building-a-reference.rst - s-tophat-round2.rst - s-more-tophat - s-even-more-tophat - s-data-analysis - m-advice From 83e55ee933bc10290b293e46dbe00a763b5fac4c Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Mon, 9 Feb 2015 11:37:27 -0600 Subject: [PATCH 09/14] trial to correct for what appeats as non-ASCII --- refTrans/m-cufflinks.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/refTrans/m-cufflinks.txt b/refTrans/m-cufflinks.txt index f195149..bbd0548 100644 --- a/refTrans/m-cufflinks.txt +++ b/refTrans/m-cufflinks.txt @@ -47,7 +47,7 @@ mapping to the whole transcriptome, 100 fragments of this million will map to on ---- Because we have multiple RNA-Seq libraries, It is recommend that you merge these assemblies into a master transcriptome. -This step is required for a differential expression analysis of the new transcripts we’ve assembled. +This step is required for a differential expression analysis of the new transcripts we've assembled. `Cuffmerge `__ performs this merge step. First step is to specify all the input GTF files in a **manifest** file listing full paths to the files. From 6e3537e71b968fe2cccf1f0424cafd4bb10a998e Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Mon, 9 Feb 2015 12:14:22 -0600 Subject: [PATCH 10/14] Fix the 2 files found with abnormal characters --- refTrans/m-cufflinks.txt | 4 ++-- refTrans/m-dataNsoftware_HPC.txt | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/refTrans/m-cufflinks.txt b/refTrans/m-cufflinks.txt index bbd0548..9fae366 100644 --- a/refTrans/m-cufflinks.txt +++ b/refTrans/m-cufflinks.txt @@ -33,11 +33,11 @@ To do this in a more conservative approach, there are several arguments that you Out folder will have four files:: - transcripts.gtf : This GTF file contains Cufflinks’ assembled isoforms. + transcripts.gtf : This GTF file contains Cufflinks' assembled isoforms. isoforms.fpkm_tracking : This file contains the estimated isoform-level expression values genes.fpkm_tracking : This file contains the estimated gene-level expression values skipped.gtf : List of skipped loci. Cufflinks skip a locus if more than 1000000 fragment - mapped to it. You can change this through the (–max-bundle-frags) + mapped to it. You can change this through the (-max-bundle-frags) Note, Cufflinks express the expression values in **FPKM** (Fragments per kilobase of exon per million reads mapped). This unit corrects for the gene size and the depth of sequencing. So let us say that the expression of a given gene was 100, diff --git a/refTrans/m-dataNsoftware_HPC.txt b/refTrans/m-dataNsoftware_HPC.txt index 7f6dee2..324338d 100644 --- a/refTrans/m-dataNsoftware_HPC.txt +++ b/refTrans/m-dataNsoftware_HPC.txt @@ -121,8 +121,8 @@ and do the same for the lung samples Link your data into a working directory --------------------------------------- -Rather than copying the files into the working directory, let’s just link them in. -this creates a reference so that UNIX knows where to find them but doesn’t need to actually move them around. +Rather than copying the files into the working directory, let's just link them in. +this creates a reference so that UNIX knows where to find them but doesn't need to actually move them around. :: ln -fs $workingPath/refTransData/dataSubset/*.fq.gz $workingPath From ef9656ae6803dce3f05cefbaea1479653ea09650 Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Sat, 14 Feb 2015 11:18:20 -0600 Subject: [PATCH 11/14] amazon implementation --- refTrans/data-snapshot.txt | 186 ++++++++++++++++++++ refTrans/igenome-backup.txt | 70 ++++++++ refTrans/index.txt | 2 +- refTrans/m-count.txt | 12 +- refTrans/m-cuffdiff.txt | 7 +- refTrans/m-cufflinks.txt | 8 +- refTrans/m-dataNsoftware_HPC.txt | 13 +- refTrans/m-dataNsoftware_amazon.txt | 253 ++++++++++++++++++++++++++++ refTrans/m-edgeR.txt | 2 + refTrans/m-htseq.txt | 4 +- refTrans/m-quality.txt | 29 +++- refTrans/m-resources.txt | 26 +++ refTrans/m-tophat.txt | 72 +++++--- 13 files changed, 649 insertions(+), 35 deletions(-) create mode 100644 refTrans/data-snapshot.txt create mode 100644 refTrans/igenome-backup.txt create mode 100644 refTrans/m-dataNsoftware_amazon.txt create mode 100644 refTrans/m-resources.txt diff --git a/refTrans/data-snapshot.txt b/refTrans/data-snapshot.txt new file mode 100644 index 0000000..85f0f6a --- /dev/null +++ b/refTrans/data-snapshot.txt @@ -0,0 +1,186 @@ +Create a data snapshot on Amazon +================================ +Essentially you create a disk; attach it; format it; +and then copy things to and from it. + + +Getting started: lunch EC2 instance with extra EBS volume +--------------------------------------------------------- +We will follow up the main events and highlight the critical points: + +1. log into your AWS account + +2. In the AWS managmenet consol: + 1. click EC2 (Amazon Elastic Compute Cloud) to open the EC2 Dashboard + 2. Change your location on the top right corner of the page to be US East (N. Virginia) + 3. press “Launch Instance” (midway down the page) + 4. Choose an Amazon Machine Image (AMI): Ubuntu Server 14.04 LTS (HVM), SSD Volume Type - ami-9a562df2 + 5. Choose an Instance Type: for now I am using the free tier t2.micro instance + 6. Next: Configure Instance Details: Keep the defaults + 7. **Next: Add Storage**: Keep the defualt storage with 8 GiB as it is and add new EBS volume of 1 GiB general purpose SSD type (to dreate data snapshot on it). Make sure it does not delete after termination + 8. Next: Tag instance: give a name to your instance (I called mine refTrans_instance) + 9. Next: Configure Security Group: Create a new security group and adjust your security rules to enable ports 22, 80, and 443 (SSH, HTTP, and HTTPS). + 10. click "review and lunch" the click "lunch" + 11. from the top drop-down menu, choose an existing key pair or create a new one and download to your disk + 12. click lunch instance + 13. scroll down and click view instance + 14. copy the Public DNS in your clipboard + +Logging into your new instance "in the cloud" (Windows version) +-------------------------------------------------------------- + 1. Download PuTTY and PuTTYgen from: http://www.chiark.greenend.org.uk/~sgtatham/putty/download.html + 2. Run PuTTYgen, Find and load your '.pem' file, and Now, "save private key" + 3. Run PuTTY, paste the public DNS in the host name & expand SSH, click Auth (do not expand), find your private key and click open + 4. click yes if probmbted then Log in as "ubuntu" + +Getting the data +---------------- + +We'll be using a few RNAseq data sets from Fagerberg et al., `Analysis +of the Human Tissue-specific Expression by Genome-wide Integration of +Transcriptomics and Antibody-based Proteomics +`__. + +You can get this data from the European Nucleotide Archive under +`ERP003613 `. +All samples in this project are `paired end +`__ . +So each sample is represented by 2 files. These files are in +`FASTQ Format `__ . + + +In this tutorial we will work with two tissues: `salivary gland +`__ and `lung +`__. Note that each +tissue has two replicates, and each replicate has two files for the paired end reads + +make a directory to download the data to the Amazon cloud +:: + + mkdir refTransData + cd refTransData + +Download the sequencing files of the salivary gland +:: + + wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR315/ERR315325/* + wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR315/ERR315382/* + +Download the sequencing files of the lung tissue +:: + + wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR315/ERR315326/* + wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR315/ERR315424/* + + +Now do:: + + ls -la refTransData/ + +You'll see something like :: + + -r--r--r-- 1 mscholz common-data 3714262571 Dec 4 08:44 ERR315325_1.fastq.gz + -r--r--r-- 1 mscholz common-data 3714262571 Dec 4 08:44 ERR315325_2.fastq.gz + ... + +which tells you that this file is about 900 MB. You can go a head and use these files +for the reset of the protocl. But for the sake of time (& memory), we will run our demo +on a subset of this dat + + +Prepare the working data on the EBS volume +------------------------------------------ + +1. Use the lsblk command to view your available disk devices and their mount points + :: + + lsblk + + + You should see something like this :: + + NAME MAJ:MIN RM SIZE RO TYPE MOUNTPOINT + xvda 202:0 0 8G 0 disk + └─xvda1 202:1 0 8G 0 part / + xvdb 202:16 0 1G 0 disk + + + | /dev/xvda1 is mounted as the root device (note the MOUNTPOINT is listed as /, the root of the Linux file system hierarchy), + | and /dev/xvdb is attached, but it has not been mounted yet. + +2. New volume does not have a file system. Use the sudo file -s device command to double check :: + + sudo file -s /dev/xvdb + + + If the output of the previous command shows simply data for the device, then there is no file system on the device and you need to create one. + +3. Create an ext4 file system on the volume :: + + sudo mkfs -t ext4 /dev/xvdb + +4. Create a mounting point directory for the volume. The mount point is where the volume is located in the file system tree and where you read and write files to after you mount the volume:: + + sudo mkdir refTransData/dataSubset + +5. Mount the volume at the location of the project folder:: + + sudo mount /dev/xvdb refTransData/dataSubset + +6. Copy in a subset of the data (100,000 reads):: + + sudo sh -c 'gunzip -c ERR315325_1.fastq.gz | head -400000 | gzip > dataSubset/salivary_repl1_R1.fq.gz' + sudo sh -c 'gunzip -c ERR315325_2.fastq.gz | head -400000 | gzip > dataSubset/salivary_repl1_R2.fq.gz' + sudo sh -c 'gunzip -c ERR315382_1.fastq.gz | head -400000 | gzip > dataSubset/salivary_repl2_R1.fq.gz' + sudo sh -c 'gunzip -c ERR315382_2.fastq.gz | head -400000 | gzip > dataSubset/salivary_repl2_R2.fq.gz' + + + and do the same for the lung samples:: + + sudo sh -c 'gunzip -c ERR315326_1.fastq.gz | head -400000 | gzip > dataSubset/lung_repl1_R1.fq.gz' + sudo sh -c 'gunzip -c ERR315326_2.fastq.gz | head -400000 | gzip > dataSubset/lung_repl1_R2.fq.gz' + sudo sh -c 'gunzip -c ERR315424_1.fastq.gz | head -400000 | gzip > dataSubset/lung_repl2_R1.fq.gz' + sudo sh -c 'gunzip -c ERR315424_2.fastq.gz | head -400000 | gzip > dataSubset/lung_repl2_R2.fq.gz' + + +Getting Set Up with the AWS Command Line Interface and make the snapshot +------------------------------------------------------------------------ +The AWS Command Line Interface is a unified tool to manage your AWS services. It will help us making a snapshot of our data and a lot more +Follow the `Amazon documenation `__. Briefly: + +1. Get your access key ID and secret access key through `IAM console `__. Do not forget to add adminstrative access role + +2. install AWS Command Line Interface:: + + cd ~ + sudo apt-get install awscli + +3. configure credentials by running:: + + aws configure + + The AWS CLI will prompt you for four pieces of information:: + + AWS Access Key ID [None]: You got it from IAM console, it should be something like AKIAIOSFODNN7EXAMPLE + AWS Secret Access Key [None]: You got it from IAM console, it should be something like wJalrXUtnFEMI/K7MDENG/bPxRfiCYEXAMPLEKEY + Default region name [None]: This is the region of your bucket (in our case it is us-east-1) + Default output format [None]: json + + +4. Create a snapshot (get the volume-id from the EC2 dashboard in the volums tab):: + + sudo umount -d /dev/xvdb + aws ec2 create-snapshot --volume-id vol-90ab2b8b --description "refTrans sample data" + + The snapshot has an ID: snap-6a5ddae5 + +5. Now we can share this snapshot but we need to modify its the permissions:: + + aws ec2 modify-snapshot-attribute --snapshot-id snap-6a5ddae5 --attribute createVolumePermission --operation-type add + +Now you can stop your instance and even delete the EBS divers. + + +---- + +Back: :doc:`m-dataNsoftware_amazon` diff --git a/refTrans/igenome-backup.txt b/refTrans/igenome-backup.txt new file mode 100644 index 0000000..697153e --- /dev/null +++ b/refTrans/igenome-backup.txt @@ -0,0 +1,70 @@ +Storage of my iGenome tarball in the Amazon S3 +============================================== + +Connect to MSU HPC and get a copy of the Bowtie2Index, genome and gtf annotation +-------------------------------------------------------------------------------- +After I logged into my AWS instance, I opened sftp connection with MSU HPC +and got the file I already prepared there in this way:: + + sftp (username)@hpc.msu.edu + get /path/to/Homo_sapiens_UCSC_hg19_small.tar.gz + exit + +Create S3 bucket +---------------- +In the AWS managmenet consol: + 1. Click S3 (Amazon simple storage service) + 2. Create Bucket + 3. Define the Bucket Name: reftransdata (Choose the Region: US Standard) + +Getting Set Up with the AWS Command Line Interface and copy a backup to S3 +-------------------------------------------------------------------------- +The AWS Command Line Interface is a unified tool to manage your AWS services. It will help us making a snapshot of our data and a lot more +Follow the `Amazon documenation `__. Briefly: + +1. Get your access key ID and secret access key through `IAM console `__. +Do not forget to add adminstrative access role + +2. install AWS Command Line Interface +:: + + cd ~ + sudo apt-get install awscli + +3. configure credentials by running +:: + + aws configure + +| The AWS CLI will prompt you for four pieces of information:: + + AWS Access Key ID [None]: You got it from IAM console, it should be something like AKIAIOSFODNN7EXAMPLE + AWS Secret Access Key [None]: You got it from IAM console, it should be something like wJalrXUtnFEMI/K7MDENG/bPxRfiCYEXAMPLEKEY + Default region name [None]: This is the region of your bucket (in our case it is us-east-1) + Default output format [None]: json + + +4. copy the tarball to your S3 space +:: + + aws s3 cp $workingPath/Homo_sapiens_UCSC_hg19_small.tar.gz s3://reftransdata/Homo_sapiens_UCSC_hg19_small.tar.gz + +.. Go to S3 consol and right click the file and choose "make it public" + + +Get a copy back from your S3 to the current EC2 +----------------------------------------------- +Install and configure the AWS Command Line Interface then copy the file to the current working directory + +:: + + aws s3 cp s3://reftransdata/Homo_sapiens_UCSC_hg19_small.tar.gz + + +open the file and change name to match the one obtained from igenome +-------------------------------------------------------------------- +:: + + tar -zxvf Homo_sapiens_UCSC_hg19_small.tar.gz + mv Homo_sapiens2 Homo_sapiens + diff --git a/refTrans/index.txt b/refTrans/index.txt index a8ab136..75bfeea 100644 --- a/refTrans/index.txt +++ b/refTrans/index.txt @@ -16,7 +16,7 @@ Tutorials: .. toctree:: :maxdepth: 2 - m-dataNsoftware_HPC + m-resources m-quality m-tophat m-count diff --git a/refTrans/m-count.txt b/refTrans/m-count.txt index 95609aa..96e6e71 100644 --- a/refTrans/m-count.txt +++ b/refTrans/m-count.txt @@ -2,7 +2,7 @@ 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: +the expression level of each expression unit. There are two main flavours of this process: .. toctree:: :maxdepth: 2 @@ -54,8 +54,10 @@ This feature might explain a lot of the differences in the output of the 2 softw 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. @@ -69,4 +71,10 @@ Links: ---- -Next: :doc:`m-m-htseq` +Next, choose one way to go: + +.. toctree:: + :maxdepth: 1 + + m-htseq + m-cufflinks diff --git a/refTrans/m-cuffdiff.txt b/refTrans/m-cuffdiff.txt index a5a7dcb..0c0612c 100644 --- a/refTrans/m-cuffdiff.txt +++ b/refTrans/m-cuffdiff.txt @@ -44,7 +44,12 @@ Now we can run the Cuffdiff --labels ${groups[0]},${groups[1]} \ $GTF_file ${gp_replicates[0]} ${gp_replicates[1]} -On HPCC, most probably you need to schadule this job +Cuffdiff usually requires high memory and would be more faster if you increased the no of threads +Check your computer resources before running it + +Cuffdif will output a folder cuffdiff_out in the working directory. The folder caintains a lot of files. +You need to read about it `here `__ + .. shell stop diff --git a/refTrans/m-cufflinks.txt b/refTrans/m-cufflinks.txt index 9fae366..0e3af71 100644 --- a/refTrans/m-cufflinks.txt +++ b/refTrans/m-cufflinks.txt @@ -18,6 +18,8 @@ As we discussed before, Cufflinks uses the ambiguously aligned reads in its calc To do this in a more conservative approach, there are several arguments that you can pass to Cufflinks: 1. --multi-read-correct Tells Cufflinks to do an initial estimation procedure to more accurately weight reads mapping to multiple locations in the genome 2. --frag-bias-correct Providing Cufflinks with a genome reference to run a bias detection and correction algorithm which improve accuracy of transcript abundance estimates. +However it worth to mention that applying these arguments will cause cufflinks to tke much more longer time to finish the job. +Forexample, with our set of data, cuffinks will need ~ 10 hours to finish compared to less than 2 hours without these arguments :: @@ -26,12 +28,12 @@ To do this in a more conservative approach, there are several arguments that you --GTF-guide "$GTF_file" \ --multi-read-correct \ --frag-bias-correct ${homo_sapiens_genome} \ - --num-threads 4 \ + --num-threads 2 \ --output-dir $f/cufflink_out \ $f/accepted_hits.bam done -Out folder will have four files:: +For each sample, We directed the outout to a folder called cufflink_out to be inside the tophat output folde. cufflink_out folder will have four files:: transcripts.gtf : This GTF file contains Cufflinks' assembled isoforms. isoforms.fpkm_tracking : This file contains the estimated isoform-level expression values @@ -63,7 +65,7 @@ Now we can run the cuffmerge software cuffmerge \ --ref-sequence ${homo_sapiens_genome} \ --ref-gtf $GTF_file \ - --num-threads 4 \ + --num-threads 2 \ assembly_GTF_list.txt The cuffmerge will output a folder called "merged_asm" which contain "merged.gtf" file diff --git a/refTrans/m-dataNsoftware_HPC.txt b/refTrans/m-dataNsoftware_HPC.txt index 324338d..fe20aae 100644 --- a/refTrans/m-dataNsoftware_HPC.txt +++ b/refTrans/m-dataNsoftware_HPC.txt @@ -16,12 +16,19 @@ Then log into a compute node e.g. dev-intel14-phi:: ssh dev-intel14-phi -Assing the path of your working directory to a ageneric variable:: + +Assign the working directory +---------------------------- +Assign the path of your working directory to a ageneric variable:: workingPath=$"" + export workingPath .. :: + workingPath=$"/mnt/ls12/Tamer/refTransProject/" + export workingPath + Getting the data ---------------- @@ -88,6 +95,8 @@ on a subset of this dat Copying in some data to work with --------------------------------- +These data files are too big and will take us hours and hourse for analysis. So let us make it simple. +We will create smaller files but getting only the first 100000 reads. .. check if the subset data directory exists .. :: @@ -132,6 +141,7 @@ this creates a reference so that UNIX knows where to find them but doesn't need groups and whatever the number of replicates. To minimize the code manipulation, you would rename your samples according to match our conventions: _repl#_R(1or2).fq.qz + install/load software --------------------- Most of softwares used through this protocol require to be installed on the root dircteory of HPC. @@ -154,6 +164,7 @@ Now let us load all the modules reuired to go through this protocol We will elaborate more on each software when we use it + .. shell stop ---- diff --git a/refTrans/m-dataNsoftware_amazon.txt b/refTrans/m-dataNsoftware_amazon.txt new file mode 100644 index 0000000..a79fdbe --- /dev/null +++ b/refTrans/m-dataNsoftware_amazon.txt @@ -0,0 +1,253 @@ +Get example data and install/load software on Amazon +==================================================== + +.. shell start + +Cloud computing involves deploying groups of remote servers and software +networks that allow centralized data storage and online access to computer +services or resources `wikipedia `__ . + +`Amazon Web Services `__ offers a broad set of cloud +comutational services. In this protocol we will implement our whole pipeline on an +Amazon Elastic Compute Cloud (EC2). + +This does mean that the first thing you need to do is get your data +over to the cloud. + +The basics +---------- + +... Amazon is happy to rent disk space to you, in addition to compute time. +They'll rent you disk space in a few different ways, but the way that's +most useful for us is through what's called Elastic Block Store. This +is essentially a hard-disk rental service. + +There are two basic concepts -- "volume" and "snapshot". A "volume" can +be thought of as a pluggable-in hard drive: you create an empty volume of +a given size, attach it to a running instance, and voila! You have extra +hard disk space. Volume-based hard disks have two problems, however: +first, they cannot be used outside of the "availability zone" they've +been created in, which means that you need to be careful to put them +in the same zone that your instance is running in; and they can't be shared +amongst people. + +Snapshots, the second concept, are the solution to transporting and +sharing the data on volumes. A "snapshot" is essentially a frozen +copy of your volume; you can copy a volume into a snapshot, and a +snapshot into a volume. + +Run through :doc:`../amazon/index` once, to get the hang of the mechanics. +Also you can read :doc:`./data-snapshot` to see how I created my data snapshot + +Launch a data snapshot +---------------------- +Go ahead and open the EC2 service to lunch a suitable instance(s). +Change your location on the top right corner of the page to be US East (N. Virginia). +Consider these notes: + +1. On "Step 1: Choose an Amazon Machine Image (AMI)": The protocol was prepared to run on "Ubuntu Server 14.04 LTS (HVM)" +2. On "Step 2: Choose an Instance Type": You can run this protocol on t2.medium instance (2 vCPU and 4 GiB memory). But when you start running your own data you can scale this up or down according to your needs +3. On "Step 3: Configure Instance Details": I just used the default settings +4. On "Step 4: Add storage": + a. start the default with 20 GiB + b. add a new new EBS volume of 1 GiB general purpose SSD type. + c. Use our snapshot (snap-6a5ddae5) to create the volume. You can search for it by the snapshot description "refTrans sample data". + d. Set the device name as /dev/sdf (be carefull about naming your device. For Linux instance, recommended device names are '/dev/sd[f-p]'. `Source: AWS Block Device Mapping Documentation `__) + e. Mark the volume with delete after termination (you always have the snapshot to make new volumes). +5. On "Step 5: Tag Instance": Give your instance a name. +6. On "Step 6: Configure Security Group": Create a new security group and adjust your security rules to enable ports 22, 80, and 443 (SSH, HTTP, and HTTPS). +7. Make a new key-paur or use yours if you already have one + +Logging into your new cloud instance (Windows version) +------------------------------------------------------ +You need an SSH client to conect from you own computer to the Amazon instance you just started in the cloud: + +1. Download PuTTY and PuTTYgen from: http://www.chiark.greenend.org.uk/~sgtatham/putty/download.html +2. Run PuTTYgen, Find and load your '.pem' file, and Now, "save private key" +3. Run PuTTY, paste the public DNS in the host name & expand SSH, click Auth (do not expand), find your private key and click open +4. click yes if probmbted then Log in as "ubuntu" +5. create a folder to contain all the files of this expermint and define this path + :: + + mkdir refTransProject + workingPath=$"/home/ubuntu/refTransProject" + export workingPath + +6. Check for the directory structure :: + + lsblk + + + You should see something like this:: + + NAME MAJ:MIN RM SIZE RO TYPE MOUNTPOINT + xvda 202:0 0 8G 0 disk + └─xvda1 202:1 0 8G 0 part / + xvdb 202:16 0 1G 0 disk + + + | /dev/xvda1 is mounted as the root device (note the MOUNTPOINT is listed as /, the root of the Linux file system hierarchy), and + | /dev/xvdf is attached, but it has not been mounted yet. + + +7. Create a mounting point and mount the new volume + :: + + mkdir $workingPath/refTransData + sudo mount /dev/xvdf $workingPath/refTransData + + + | Check using lsblk again to make sure everything is ok. + +Now, link the data to the working directory. +This creates a reference so that UNIX knows where to find them but doesn't need to actually move them around. +:: + + ln -fs $workingPath/refTransData/*.fq.gz $workingPath + + +Now if you typed:: + + ls refTransProject/ + +You should see the 8 fastq files in your folder + +.. note:: + + You can use your own data instead. The code can be applied to any paired end data set whatever the number of the studied + groups and whatever the number of replicates. To minimize the code manipulation, you would rename your samples according to match + our conventions: _repl#_R(1or2).fq.qz + +install software and link to the working directory +-------------------------------------------------- + +.. note:: + During the installation, copy the commands to your prompt one line + at a time to make sure that every line pass safelly and because some + commands are interactive requiring your confirmation to continue + +install prerequisites +:: + + cd ~ + sudo apt-get install gcc g++ pkg-config wget make + sudo add-apt-repository ppa:webupd8team/java + sudo apt-get update + sudo apt-get install oracle-java7-installer + +.. Source http://www.webupd8.org/2012/01/install-oracle-java-jdk-7-in-ubuntu-via.html + +install FastQC/0.11.2 +:: + + cd ~ + wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.2.zip + sudo apt-get install unzip + unzip fastqc_v0.11.2.zip + cd FastQC + chmod 755 fastqc + sudo ln -s /home/ubuntu/FastQC/fastqc /usr/local/bin/fastqc + +install Trimmomatic/0.32 +:: + + cd ~ + wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.32.zip + unzip Trimmomatic-0.32.zip + mv /home/ubuntu/Trimmomatic-0.32/trimmomatic-0.32.jar /home/ubuntu/Trimmomatic-0.32/trimmomatic + TRIM=$"/home/ubuntu/Trimmomatic-0.32" + export TRIM + + +install Bowtie/2.2.3.0 +:: + + cd ~ + wget http://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.3/bowtie2-2.2.3-source.zip + unzip bowtie2-2.2.3-source.zip + cd bowtie2-2.2.3 + make + PATH=$PATH:/home/ubuntu/bowtie2-2.2.3 + export PATH + +install TopHat2/2.0.13 +:: + + cd /usr/src/ + sudo wget http://ccb.jhu.edu/software/tophat/downloads/tophat-2.0.13.Linux_x86_64.tar.gz + sudo tar -zxvf tophat-2.0.13.Linux_x86_64.tar.gz + sudo rm tophat-2.0.13.Linux_x86_64.tar.gz + cd /usr/bin + sudo ln -s /usr/src/tophat-2.0.13.Linux_x86_64/tophat2 ./tophat + +install samtools/0.1.19-1 +:: + + cd ~ + sudo apt-get install samtools + +install HTSeq/0.6.1 +:: + + cd ~ + sudo apt-get install build-essential python2.7-dev python-numpy python-matplotlib + wget https://pypi.python.org/packages/source/H/HTSeq/HTSeq-0.6.1.tar.gz + tar -zxvf HTSeq-0.6.1.tar.gz + cd HTSeq-0.6.1 + python setup.py install --user + sudo cp .local/bin/* /usr/bin/. + +install PySam +:: + + sudo apt-get install python-pip + sudo pip install pysam + +Install GNU R core and the edgeR package +:: + + cd ~ + sudo apt-get install r-base-core + echo " + source('http://bioconductor.org/biocLite.R') + biocLite('edgeR')" > install_edgeR.R + sudo Rscript install_edgeR.R + +install cufflinks/2.1.1 +:: + + cd /usr/src/ + sudo wget http://cole-trapnell-lab.github.io/cufflinks/assets/downloads/cufflinks-2.1.1.Linux_x86_64.tar.gz + sudo tar -zxvf cufflinks-2.1.1.Linux_x86_64.tar.gz + sudo rm cufflinks-2.1.1.Linux_x86_64.tar.gz + cd /usr/bin + sudo ln -s /usr/src/cufflinks-2.1.1.Linux_x86_64/cufflinks . + sudo ln -s /usr/src/cufflinks-2.1.1.Linux_x86_64/cuffmerge . + sudo ln -s /usr/src/cufflinks-2.1.1.Linux_x86_64/cuffcompare . + sudo ln -s /usr/src/cufflinks-2.1.1.Linux_x86_64/cuffdiff . + sudo ln -s /usr/src/cufflinks-2.1.1.Linux_x86_64/gffread . + sudo ln -s /usr/src/cufflinks-2.1.1.Linux_x86_64/gtf_to_sam . + +install tmux + sudo apt-get install tmux + +.. Install FastX/0.0.13.2:: + cd ~ + sudo apt-get install gcc g++ pkg-config wget make + wget http://hannonlab.cshl.edu/fastx_toolkit/libgtextutils-0.6.1.tar.bz2 + tar -xjf libgtextutils-0.6.1.tar.bz2 + cd libgtextutils-0.6.1 + sudo sh -c "./configure && make && make install" + cd ~ + wget http://hannonlab.cshl.edu/fastx_toolkit/fastx_toolkit-0.0.13.2.tar.bz2 + tar xjf fastx_toolkit-0.0.13.2.tar.bz2 + cd fastx_toolkit-0.0.13.2/ + sudo sh -c "./configure && make && make install" + export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH + + +.. shell stop + +---- + +Next: :doc:`m-quality` diff --git a/refTrans/m-edgeR.txt b/refTrans/m-edgeR.txt index ba3fc24..d5bba7c 100644 --- a/refTrans/m-edgeR.txt +++ b/refTrans/m-edgeR.txt @@ -7,6 +7,8 @@ We'll be using `edgeR `__ to do the basic differential expression analysis of our counts. +.. note:: This section reuires GNU R core installed with edgeR package installed + To run edgeR, you need to write a data loading and manipulation script in R. This script will load in all the studied samples, execute an MA plot, and provide a spreadsheet with differential expression of the studied tissues. diff --git a/refTrans/m-htseq.txt b/refTrans/m-htseq.txt index 2508756..9900808 100644 --- a/refTrans/m-htseq.txt +++ b/refTrans/m-htseq.txt @@ -54,7 +54,7 @@ And next, run HTSeq When this is done, let us view one of the output file:: - less salivary_repl1_tophat/salivary_repl1_counts.txt + less salivary_repl1_tophat/salivary_repl1_PEcounts.txt (use 'q' to exit). These are your gene counts. @@ -78,4 +78,4 @@ DEseq and edgeR want exactly this kind of information! ---- -Next: :doc:`m-data-analysis` +Next: :doc:`m-edgeR` diff --git a/refTrans/m-quality.txt b/refTrans/m-quality.txt index 1ac24e5..15a7933 100644 --- a/refTrans/m-quality.txt +++ b/refTrans/m-quality.txt @@ -3,6 +3,15 @@ Short read quality and trimming .. shell start +It is time to run your Tmux session and it is a good trend to do this from your root directory:: + + cd ~ + tmux + +Now, the remote computer will sustain your jobs even your ssh connection was interrupted. +You can always restore your tmux session by going to your root directory and type "tmux attach" + + Assessment of sequence quality (FastQC) --------------------------------------- @@ -14,15 +23,16 @@ FastQC operates on Fastq files as well as GZip compressed Fastq (.gz) files. Now run FastQC on your sequencing file :: + cd $workingPath mkdir rawData_fastqc - fastqc *.fq.gz --outdir=./rawData_fastqc + fastqc *.fq.gz --outdir=$workingPath/rawData_fastqc .. .. note:: When you pass multiple files to FastQC, it operates one file at a time give it a seprate output. Now type 'ls' in the output directory to explore the result files:: - ls ./rawData_fastqc + ls $workingPath/rawData_fastqc and you will see :: @@ -30,9 +40,15 @@ and you will see :: salivary_repl1_R1_fastqc.zip ... -In the root of your terminal, you can copy these files to your laptop and open them in a browser:: +In the root of your terminal, you can copy these files to your laptop and open them in a browser + +For HPCC:: scp :/rawData_fastqc/salivary*fastqc.* /tmp + +For Amazon:: + + scp -i ubuntu@:/path/to/file . and then open the html files in your browser. @@ -40,8 +56,8 @@ Links: * `FastQC tutorial video `__ -Trimmomatic ------------ +Trimming of low quality bases (Trimmomatic) +------------------------------------------- Now we're going to do some trimming to remove any leftover adaptor sequences and low quality trailing edges. We'll be using `Trimmomatic `__. @@ -53,6 +69,7 @@ For a discussion of optimal RNAseq trimming strategies, see `MacManes, 2014 To run Trimmomatic :: + cd $workingPath for f in *R1*.gz; do input_R1="$f" input_R2=$(echo "$f" | sed s/R1/R2/) @@ -107,7 +124,7 @@ Run FastQC again :: mkdir trimmedData_fastqc - fastqc *.pe.fq.gz --outdir=./trimmedData_fastqc + fastqc *.pe.fq.gz --outdir=$workingPath/trimmedData_fastqc .. * Does it matter that you still have adapters!? diff --git a/refTrans/m-resources.txt b/refTrans/m-resources.txt new file mode 100644 index 0000000..b582a4f --- /dev/null +++ b/refTrans/m-resources.txt @@ -0,0 +1,26 @@ +High Performance Computation +============================ + +During the last few years, Next generation sequencing (NGS) has become widely available and cost effective. +The tremendous amounts of data generated by NGS are fundamentally changing how research is done. +However, these huge amounts of data requires dedicated software and special computation resources. +High Performance Computaters are tyically used for this kind of data analysis. + +We adapted this protocol to run on 2 examples of those computers + +.. toctree:: + :maxdepth: 2 + + m-dataNsoftware_HPC + m-dataNsoftware_amazon + + +Using either system, you need to communicate from your local machine to run your job on a distant server. +We use `SSH client `__ to do this job. +A problem that you have to be aware of is that when you log out of the communication session, +whatever the job you are running on the distant computer will stop. Some of our jobs should take hours and may be days +and it will be too frustrating to start over each time you lose the connection with the server. +So. We need a way to trick the server to continue our jobs even if we are not logged in. +Couple softwares do that, for example `Tmux `__. +Tmux will let us begin a session and detach from it, letting it run by itself. + diff --git a/refTrans/m-tophat.txt b/refTrans/m-tophat.txt index d4be2ff..e662681 100644 --- a/refTrans/m-tophat.txt +++ b/refTrans/m-tophat.txt @@ -8,7 +8,7 @@ reads to the reference gene set, for the purpose of counting how many reads have come from each gene. We'll be using the `TopHat software `__ -.. note:: This section reuires TopHat2/2.0.12, Bowtie/2.2.3.0 and Samtools/0.1.19.0 installed +.. note:: This section reuires TopHat2/2.0.13, Bowtie/2.2.3.0 and Samtools/0.1.19.0 installed Download the Bowtie2 index -------------------------- @@ -24,9 +24,7 @@ a given genome assembly. The `Illumina iGenomes project `__ provided the RNA-Seq user community with a set of genome sequence indexes (including Bowtie indexes). -We will download the iGenome index for the UCSC hg19 assembly version and -define the path of the genome and genome_index_base which is the name of any -of the index files up to but not including the first period. +We will download the iGenome index for the UCSC hg19 assembly version. .. check if the data directory exists .. :: @@ -35,15 +33,36 @@ of the index files up to but not including the first period. :: + cd $workingPath wget ftp://igenome:G3nom3s4u@ussd-ftp.illumina.com/Homo_sapiens/UCSC/hg19/Homo_sapiens_UCSC_hg19.tar.gz tar -zxvf Homo_sapiens_UCSC_hg19.tar.gz - homo_sapiens_genome="$workingPath/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa" - Bowtie2Index="$workingPath/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome" - + .. :: fi +.. note:: + + The homo sapians iGenome tarball is huge (~20 GiB). I have prepared a smaller version containing the minimal required files. + This file is stored in Amazon S3. To use this file replace the last 3 lines of code with:: + + cd $workingPath + wget https://s3.amazonaws.com/reftransdata/Homo_sapiens_UCSC_hg19_small.tar.gz + tar -zxvf Homo_sapiens_UCSC_hg19_small.tar.gz + mv Homo_sapiens2 Homo_sapiens + +.. for more details about creating and using this file, check this file: :doc:`./igenome-backup` + + +Define the path of the genome and genome_index_base which is the name of any +of the index files up to but not including the first period. + +:: + + homo_sapiens_genome="$workingPath/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa" + Bowtie2Index="$workingPath/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome" + + Making use of gene model annotations ------------------------------------ Model organisms usually have a good annotation of their transcriptome. Annotation @@ -72,7 +91,7 @@ To create the transcriptome-index, we can inovke TopHat without input reads :: tophat \ - --num-threads 4 \ + --num-threads 1 \ --GTF "$GTF_file" \ --transcriptome-index "$workingPath/trans_Index/transcriptome" \ $Bowtie2Index @@ -81,7 +100,7 @@ To create the transcriptome-index, we can inovke TopHat without input reads Prepare the input reads ----------------------- You can pass these two PE files directly to the TopHat to do the pair end orinted alignment. -TopHat allows the use of additional unpaired reads but it has to pass as single file. +TopHat allows the use of additional unpaired reads but it has to pass as a single file. So we need to merge every two SE files belonging to the same sample :: @@ -98,14 +117,14 @@ Mapping reads And now run TopHat on one sample:: tophat \ - --num-threads 4 \ + --num-threads 2 \ --GTF "$GTF_file" \ --transcriptome-index "$workingPath/trans_Index/transcriptome" \ --output-dir salivary_repl1_tophat \ $Bowtie2Index \ salivary_repl1_R1.pe.fq.gz salivary_repl1_R2.pe.fq.gz,salivary_repl1_R.se.fq.gz -This will take about ~10-15 minutes. We can run all the samples sequencially using a for loop +This will take about ~10-15 minutes. We can run all the samples sequencially using a for loop. :: for f in *R1.pe*; do @@ -114,7 +133,7 @@ This will take about ~10-15 minutes. We can run all the samples sequencially usi input_R=$(echo "$f" | sed s/R1.pe/R.se/) output=$(basename "$f" | cut -f 1,2 -d "_")_tophat tophat \ - --num-threads 4 \ + --num-threads 1 \ --GTF "$GTF_file" \ --transcriptome-index "$workingPath/trans_Index/transcriptome" \ --output-dir $output \ @@ -124,20 +143,35 @@ This will take about ~10-15 minutes. We can run all the samples sequencially usi .. Note:: - With large input files, TopHat will take hours and may be days. It will be more - efficient to run TopHat on all samples in parallel jobs + With large input files, TopHat will take hours and may be days. You can increase the + number of threads but make sure that you have the required underlying resources. + It will be more efficient to run TopHat on all samples in parallel jobs. On Amazon AWS, + you can run multiple instances. For HPC, you can submet multiple jobs simultaneously. :doc:`./m-HPC-tophat` +The output will be a folder for each sample names like _tophat. Each folder contain multiple files + +1. accepted_hits.bam: The main output file with contains the read alignment in `SAM `__ format. +2. unmapped.bam: Another file in the same format but contains only the unmapped reads +3. junctions.bed, insertions.bed, deletions.bed: A `UCSC BED track `__ of junctions, insertions and deletions reported by TopHat +4. align_summary.txt: a text file that contains summry statistics about the mapped reads + + Counting mapped reads percentage -------------------------------- -Let's ask samtools for the total number of reads that mapped:: +You can open the align_summary file in any out folder :: + + vim salivary_repl1_tophat/align_summary.txt + +We can calculte these statistics ourselves. Let's ask samtools for the total number of reads that mapped. +We count all the alignment in the BAM file but exclude unmapped reads (-F 4) and also exclude non primary alignments (-F 256):: - samtools view -c -F 4 salivary_repl1/accepted_hits.bam + samtools view -c -F 4 -F 256 salivary_repl1_tophat/accepted_hits.bam -You should get around 186,256. +You should get around 170,849. We can count the reads in our input files:: @@ -148,10 +182,10 @@ We can count the reads in our input files:: You should get around 189,150 -So you have 186,256 mapped reads out of 189,150 -- about 98%. That's pretty good! +So you have 170,849 mapped reads out of 189,150 -- about 90.3%. That's pretty good! .. shell stop ---- -Next: :doc:`m-htseq` +Next: :doc:`m-count` From fff3452ff36644fa0730c7f0adb524ef27f2b272 Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Sat, 14 Feb 2015 17:10:38 -0600 Subject: [PATCH 12/14] fix the m-count page hierarchy --- refTrans/m-count.txt | 31 +++++++++---------------------- 1 file changed, 9 insertions(+), 22 deletions(-) diff --git a/refTrans/m-count.txt b/refTrans/m-count.txt index 96e6e71..fdff133 100644 --- a/refTrans/m-count.txt +++ b/refTrans/m-count.txt @@ -4,17 +4,9 @@ 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 are two main flavours of this process: -.. toctree:: - :maxdepth: 2 - - m-htseq - m-cufflinks - - ----- - -:doc:`./m-htseq` +Counting the reads that map to each gene +---------------------------------------- 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 `__. @@ -27,10 +19,12 @@ Annotaions file of model oraganisms are avaliable on all gene banks (e.g. Also they are avaliable through other genomic software projects e.g. `Illumina iGenomes project `__. ----- +Go for: +:doc:`./m-htseq` -:doc:`./m-cufflinks` +Counting the reads that map to each transcript +---------------------------------------------- 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 @@ -62,19 +56,12 @@ On the other hand, mapping to transcripts requires the implementation of error p approaches and thus requires rigorous laboratory confirmation. +Go for: +:doc:`./m-cufflinks` + ---- Links: * `Trapnell et al. 2013 Differential analysis of gene regulation at transcript resolution with RNA-seq `__ - ----- - -Next, choose one way to go: - -.. toctree:: - :maxdepth: 1 - - m-htseq - m-cufflinks From a2e24f7bd98c73e4c7db3e43339680a049f4f1b0 Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Tue, 3 Mar 2015 00:58:53 -0500 Subject: [PATCH 13/14] update typos and fix installation on amazon --- index.txt | 3 - refTrans/data-snapshot.txt | 34 ++++++++---- refTrans/m-dataNsoftware_HPC.txt | 4 +- refTrans/m-dataNsoftware_amazon.txt | 85 +++++++++++++++-------------- refTrans/m-quality.txt | 3 + refTrans/m-resources.txt | 24 ++++---- 6 files changed, 86 insertions(+), 67 deletions(-) diff --git a/index.txt b/index.txt index bc54d8e..4b2e804 100644 --- a/index.txt +++ b/index.txt @@ -58,9 +58,6 @@ Reference based RNAseq assembly: the refTrans Protocol 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 `__ -and made use of the `MSU High Performance Compute Center `__ . - Additional information ---------------------- diff --git a/refTrans/data-snapshot.txt b/refTrans/data-snapshot.txt index 85f0f6a..0b4f21e 100644 --- a/refTrans/data-snapshot.txt +++ b/refTrans/data-snapshot.txt @@ -1,5 +1,8 @@ Create a data snapshot on Amazon ================================ + +.. shell start + Essentially you create a disk; attach it; format it; and then copy things to and from it. @@ -108,26 +111,31 @@ Prepare the working data on the EBS volume | /dev/xvda1 is mounted as the root device (note the MOUNTPOINT is listed as /, the root of the Linux file system hierarchy), | and /dev/xvdb is attached, but it has not been mounted yet. -2. New volume does not have a file system. Use the sudo file -s device command to double check :: +2. New volume does not have a file system. Use the sudo file -s device command to double check + :: sudo file -s /dev/xvdb If the output of the previous command shows simply data for the device, then there is no file system on the device and you need to create one. -3. Create an ext4 file system on the volume :: +3. Create an ext4 file system on the volume + :: sudo mkfs -t ext4 /dev/xvdb -4. Create a mounting point directory for the volume. The mount point is where the volume is located in the file system tree and where you read and write files to after you mount the volume:: +4. Create a mounting point directory for the volume. The mount point is where the volume is located in the file system tree and where you read and write files to after you mount the volume + :: sudo mkdir refTransData/dataSubset -5. Mount the volume at the location of the project folder:: +5. Mount the volume at the location of the project folder + :: sudo mount /dev/xvdb refTransData/dataSubset -6. Copy in a subset of the data (100,000 reads):: +6. Copy in a subset of the data (100,000 reads) + :: sudo sh -c 'gunzip -c ERR315325_1.fastq.gz | head -400000 | gzip > dataSubset/salivary_repl1_R1.fq.gz' sudo sh -c 'gunzip -c ERR315325_2.fastq.gz | head -400000 | gzip > dataSubset/salivary_repl1_R2.fq.gz' @@ -135,7 +143,8 @@ Prepare the working data on the EBS volume sudo sh -c 'gunzip -c ERR315382_2.fastq.gz | head -400000 | gzip > dataSubset/salivary_repl2_R2.fq.gz' - and do the same for the lung samples:: + and do the same for the lung samples + :: sudo sh -c 'gunzip -c ERR315326_1.fastq.gz | head -400000 | gzip > dataSubset/lung_repl1_R1.fq.gz' sudo sh -c 'gunzip -c ERR315326_2.fastq.gz | head -400000 | gzip > dataSubset/lung_repl1_R2.fq.gz' @@ -150,12 +159,14 @@ Follow the `Amazon documenation `__. Do not forget to add adminstrative access role -2. install AWS Command Line Interface:: +2. install AWS Command Line Interface + :: cd ~ sudo apt-get install awscli -3. configure credentials by running:: +3. configure credentials by running + :: aws configure @@ -167,19 +178,22 @@ Follow the `Amazon documenation `__ . +services or resources (`wikipedia `__) . `Amazon Web Services `__ offers a broad set of cloud comutational services. In this protocol we will implement our whole pipeline on an @@ -43,20 +43,21 @@ Launch a data snapshot ---------------------- Go ahead and open the EC2 service to lunch a suitable instance(s). Change your location on the top right corner of the page to be US East (N. Virginia). -Consider these notes: +Push the **Launch Instance** button then follow these instructions: 1. On "Step 1: Choose an Amazon Machine Image (AMI)": The protocol was prepared to run on "Ubuntu Server 14.04 LTS (HVM)" -2. On "Step 2: Choose an Instance Type": You can run this protocol on t2.medium instance (2 vCPU and 4 GiB memory). But when you start running your own data you can scale this up or down according to your needs -3. On "Step 3: Configure Instance Details": I just used the default settings +2. On "Step 2: Choose an Instance Type": You can run this protocol on **t2.medium instance** (2 vCPU and 4 GiB memory). But when you start running your own data you can scale this up or down according to your needs +3. On "Step 3: Configure Instance Details": We will accept the default settings 4. On "Step 4: Add storage": - a. start the default with 20 GiB - b. add a new new EBS volume of 1 GiB general purpose SSD type. - c. Use our snapshot (snap-6a5ddae5) to create the volume. You can search for it by the snapshot description "refTrans sample data". + a. Start the default volume with 20 GiB + b. Add a new new EBS volume of 1 GiB general purpose SSD type. + c. Use our snapshot (snap-6a5ddae5) to create the additional volume. You can search for it by the snapshot description "refTrans sample data". d. Set the device name as /dev/sdf (be carefull about naming your device. For Linux instance, recommended device names are '/dev/sd[f-p]'. `Source: AWS Block Device Mapping Documentation `__) - e. Mark the volume with delete after termination (you always have the snapshot to make new volumes). + e. Mark the volume with **delete after termination** (you always have the snapshot to make new volumes). 5. On "Step 5: Tag Instance": Give your instance a name. -6. On "Step 6: Configure Security Group": Create a new security group and adjust your security rules to enable ports 22, 80, and 443 (SSH, HTTP, and HTTPS). -7. Make a new key-paur or use yours if you already have one +6. On "Step 6: Configure Security Group": Create a new security group and add security rules to enable ports 22, 80, and 443 (SSH, HTTP, and HTTPS). +7. On "Step 7: Review Instance Launch": Review the information of your instance. You will see an alarm about the security of your instnce. This is ok, once you click launch, you will be able to make a key-paur (or use yours if you already have one) +8. After you launch your instance, the confirmation page will show your instance ID. Click the instance ID to watch your instance status. In the lower half of the page, you will description info of your instance. Copy the puplic IP address to use in the next step. Logging into your new cloud instance (Windows version) ------------------------------------------------------ @@ -64,7 +65,7 @@ You need an SSH client to conect from you own computer to the Amazon instance yo 1. Download PuTTY and PuTTYgen from: http://www.chiark.greenend.org.uk/~sgtatham/putty/download.html 2. Run PuTTYgen, Find and load your '.pem' file, and Now, "save private key" -3. Run PuTTY, paste the public DNS in the host name & expand SSH, click Auth (do not expand), find your private key and click open +3. Run PuTTY, paste the public DNS in the host name. In the left category panel, expand the SSH entry, select Auth, find your private key and click open 4. click yes if probmbted then Log in as "ubuntu" 5. create a folder to contain all the files of this expermint and define this path :: @@ -130,7 +131,9 @@ install prerequisites :: cd ~ + sudo apt-get update sudo apt-get install gcc g++ pkg-config wget make + sudo apt-get install unzip sudo add-apt-repository ppa:webupd8team/java sudo apt-get update sudo apt-get install oracle-java7-installer @@ -140,34 +143,34 @@ install prerequisites install FastQC/0.11.2 :: - cd ~ - wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.2.zip - sudo apt-get install unzip - unzip fastqc_v0.11.2.zip + cd /usr/src/ + sudo wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.2.zip + sudo unzip fastqc_v0.11.2.zip cd FastQC - chmod 755 fastqc - sudo ln -s /home/ubuntu/FastQC/fastqc /usr/local/bin/fastqc + sudo chmod 755 fastqc + sudo ln -s /usr/src/FastQC/fastqc /usr/local/bin/fastqc + install Trimmomatic/0.32 :: - cd ~ - wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.32.zip - unzip Trimmomatic-0.32.zip - mv /home/ubuntu/Trimmomatic-0.32/trimmomatic-0.32.jar /home/ubuntu/Trimmomatic-0.32/trimmomatic - TRIM=$"/home/ubuntu/Trimmomatic-0.32" + cd /usr/src/ + sudo wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.32.zip + sudo unzip Trimmomatic-0.32.zip + sudo mv /usr/src/Trimmomatic-0.32/trimmomatic-0.32.jar /usr/src/Trimmomatic-0.32/trimmomatic + TRIM=$"/usr/src/Trimmomatic-0.32" export TRIM install Bowtie/2.2.3.0 :: - cd ~ - wget http://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.3/bowtie2-2.2.3-source.zip - unzip bowtie2-2.2.3-source.zip + cd /usr/src/ + sudo wget http://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.3/bowtie2-2.2.3-source.zip + sudo unzip bowtie2-2.2.3-source.zip cd bowtie2-2.2.3 - make - PATH=$PATH:/home/ubuntu/bowtie2-2.2.3 + sudo make + PATH=$PATH:/usr/src/bowtie2-2.2.3 export PATH install TopHat2/2.0.13 @@ -176,7 +179,6 @@ install TopHat2/2.0.13 cd /usr/src/ sudo wget http://ccb.jhu.edu/software/tophat/downloads/tophat-2.0.13.Linux_x86_64.tar.gz sudo tar -zxvf tophat-2.0.13.Linux_x86_64.tar.gz - sudo rm tophat-2.0.13.Linux_x86_64.tar.gz cd /usr/bin sudo ln -s /usr/src/tophat-2.0.13.Linux_x86_64/tophat2 ./tophat @@ -189,17 +191,21 @@ install samtools/0.1.19-1 install HTSeq/0.6.1 :: - cd ~ + cd /usr/src/ sudo apt-get install build-essential python2.7-dev python-numpy python-matplotlib - wget https://pypi.python.org/packages/source/H/HTSeq/HTSeq-0.6.1.tar.gz - tar -zxvf HTSeq-0.6.1.tar.gz + sudo wget https://pypi.python.org/packages/source/H/HTSeq/HTSeq-0.6.1.tar.gz + sudo tar -zxvf HTSeq-0.6.1.tar.gz cd HTSeq-0.6.1 - python setup.py install --user - sudo cp .local/bin/* /usr/bin/. + sudo python setup.py install --user + PATH=$PATH:/home/ubuntu/.local/bin + export PATH + sudo chmod -R 755 ~/.local + install PySam :: + cd ~ sudo apt-get install python-pip sudo pip install pysam @@ -208,10 +214,7 @@ Install GNU R core and the edgeR package cd ~ sudo apt-get install r-base-core - echo " - source('http://bioconductor.org/biocLite.R') - biocLite('edgeR')" > install_edgeR.R - sudo Rscript install_edgeR.R + sudo Rscript -e "source('http://bioconductor.org/biocLite.R'); biocLite('edgeR');" install cufflinks/2.1.1 :: @@ -219,7 +222,6 @@ install cufflinks/2.1.1 cd /usr/src/ sudo wget http://cole-trapnell-lab.github.io/cufflinks/assets/downloads/cufflinks-2.1.1.Linux_x86_64.tar.gz sudo tar -zxvf cufflinks-2.1.1.Linux_x86_64.tar.gz - sudo rm cufflinks-2.1.1.Linux_x86_64.tar.gz cd /usr/bin sudo ln -s /usr/src/cufflinks-2.1.1.Linux_x86_64/cufflinks . sudo ln -s /usr/src/cufflinks-2.1.1.Linux_x86_64/cuffmerge . @@ -229,6 +231,9 @@ install cufflinks/2.1.1 sudo ln -s /usr/src/cufflinks-2.1.1.Linux_x86_64/gtf_to_sam . install tmux +:: + + cd ~ sudo apt-get install tmux .. Install FastX/0.0.13.2:: diff --git a/refTrans/m-quality.txt b/refTrans/m-quality.txt index 15a7933..f71cc1b 100644 --- a/refTrans/m-quality.txt +++ b/refTrans/m-quality.txt @@ -126,6 +126,9 @@ Run FastQC again mkdir trimmedData_fastqc fastqc *.pe.fq.gz --outdir=$workingPath/trimmedData_fastqc + +Now you can dowload the output files and check for changes after Trimmomatic treatment + .. * Does it matter that you still have adapters!? .. shell stop diff --git a/refTrans/m-resources.txt b/refTrans/m-resources.txt index b582a4f..aad6649 100644 --- a/refTrans/m-resources.txt +++ b/refTrans/m-resources.txt @@ -3,10 +3,9 @@ High Performance Computation During the last few years, Next generation sequencing (NGS) has become widely available and cost effective. The tremendous amounts of data generated by NGS are fundamentally changing how research is done. -However, these huge amounts of data requires dedicated software and special computation resources. -High Performance Computaters are tyically used for this kind of data analysis. +However, these huge amounts of data require dedicated software and special computational resources. -We adapted this protocol to run on 2 examples of those computers +We adapted this protocol to run on 2 examples of computational plateforms .. toctree:: :maxdepth: 2 @@ -14,13 +13,14 @@ We adapted this protocol to run on 2 examples of those computers m-dataNsoftware_HPC m-dataNsoftware_amazon - -Using either system, you need to communicate from your local machine to run your job on a distant server. -We use `SSH client `__ to do this job. -A problem that you have to be aware of is that when you log out of the communication session, -whatever the job you are running on the distant computer will stop. Some of our jobs should take hours and may be days -and it will be too frustrating to start over each time you lose the connection with the server. -So. We need a way to trick the server to continue our jobs even if we are not logged in. -Couple softwares do that, for example `Tmux `__. -Tmux will let us begin a session and detach from it, letting it run by itself. +.. note:: + + Using either system, you need to communicate from your local machine to run your job on a distant server. + We use `SSH client `__ to do this job. + A problem that you have to be aware of is that when you log out of the communication session, + whatever the job you are running on the distant computer will stop. Some of our jobs should take hours and may be days + and it will be too frustrating to start over each time you lose the connection with the server. + So. We need a way to trick the server to continue our jobs even if we are not logged in. + Couple softwares do that, for example `Tmux `__. + Tmux will let us begin a session and detach from it, letting it run by itself. From deac8aca94efdead66cd9291385a280c65e79ae9 Mon Sep 17 00:00:00 2001 From: drtamermansour Date: Mon, 2 Nov 2015 05:41:06 -0500 Subject: [PATCH 14/14] fix the name of amazon drive --- refTrans/data-snapshot.txt | 15 ++++++++------- refTrans/m-dataNsoftware_amazon.txt | 2 +- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/refTrans/data-snapshot.txt b/refTrans/data-snapshot.txt index 0b4f21e..430b9bd 100644 --- a/refTrans/data-snapshot.txt +++ b/refTrans/data-snapshot.txt @@ -20,7 +20,7 @@ We will follow up the main events and highlight the critical points: 4. Choose an Amazon Machine Image (AMI): Ubuntu Server 14.04 LTS (HVM), SSD Volume Type - ami-9a562df2 5. Choose an Instance Type: for now I am using the free tier t2.micro instance 6. Next: Configure Instance Details: Keep the defaults - 7. **Next: Add Storage**: Keep the defualt storage with 8 GiB as it is and add new EBS volume of 1 GiB general purpose SSD type (to dreate data snapshot on it). Make sure it does not delete after termination + 7. **Next: Add Storage**: Keep the defualt storage with 8 GiB as it is and add new EBS volume of 1 GiB general purpose SSD type (to create data snapshot on it). Make sure it does not delete after termination 8. Next: Tag instance: give a name to your instance (I called mine refTrans_instance) 9. Next: Configure Security Group: Create a new security group and adjust your security rules to enable ports 22, 80, and 443 (SSH, HTTP, and HTTPS). 10. click "review and lunch" the click "lunch" @@ -105,16 +105,16 @@ Prepare the working data on the EBS volume NAME MAJ:MIN RM SIZE RO TYPE MOUNTPOINT xvda 202:0 0 8G 0 disk └─xvda1 202:1 0 8G 0 part / - xvdb 202:16 0 1G 0 disk + xvdf 202:16 0 1G 0 disk | /dev/xvda1 is mounted as the root device (note the MOUNTPOINT is listed as /, the root of the Linux file system hierarchy), - | and /dev/xvdb is attached, but it has not been mounted yet. + | and /dev/xvdf is attached, but it has not been mounted yet. 2. New volume does not have a file system. Use the sudo file -s device command to double check :: - sudo file -s /dev/xvdb + sudo file -s /dev/xvdf If the output of the previous command shows simply data for the device, then there is no file system on the device and you need to create one. @@ -122,7 +122,7 @@ Prepare the working data on the EBS volume 3. Create an ext4 file system on the volume :: - sudo mkfs -t ext4 /dev/xvdb + sudo mkfs -t ext4 /dev/xvdf 4. Create a mounting point directory for the volume. The mount point is where the volume is located in the file system tree and where you read and write files to after you mount the volume :: @@ -132,7 +132,7 @@ Prepare the working data on the EBS volume 5. Mount the volume at the location of the project folder :: - sudo mount /dev/xvdb refTransData/dataSubset + sudo mount /dev/xvdf refTransData/dataSubset 6. Copy in a subset of the data (100,000 reads) :: @@ -163,6 +163,7 @@ Follow the `Amazon documenation