Alignment of RNA Seq data with STAR

The goal of this exercise is to align FASTQ files to the mouse genome using STAR. The summary for this project describes the background information.

The STAR manual has a description of how to run the program and all the options you can supply to STAR. It's a good idea to look this over.

To run STAR on your FASTQ files you generally need three things:

• the location of your FASTQ files
• a STAR index of the genome
• a GTF file of transcript descriptions

The FASTQ files are located in the class Data directory:

/home/cws/CompGenomics/Data/RNA-Seq/

What are these files? How can we view them?

At the Stowers, we have many genome indexes for various organisms that have already been created for you. But since we're not on the internal network, our STAR alignment index is in:

/home/cws/CompGenomics/dm6/

You can see that this directory points to resources around the dm6 version of the fly genome, and these indices have been constructed using version 98 of the Ensembl annotations for mouse genes. There is a GTF file of transcript descriptions there. Can you find it?

The general structure of the STAR command to align some data looks as follows:

STAR NumberOfProcessors IndexDirectory FASTQ_files

We will be specifying these options, and a few others to the program using command-line arguments that begin with a double-dash, as described in the following table:

When working out a complex command, it's a good idea to test it on some sample data. Let's use our unix skills to grab some data and make a toy data file to play with as input.

There's a data file in the class data directory, we can non-destructively uncompress it and take the top 1000 lines, and then re-direct that to a local file in our current directory:

# grab 1000 lines of fastq data from a file
gunzip -c /home/cws/CompGenomics/Data/RNA-Seq/gd7_1.fastq.gz | head -1000 > sample.fastq

For now let's just begin with the simplest example possible: align our toy data using STAR

STAR --runThreadN 2 --genomeDir /home/cws/CompGenomics/dm6/Ens_98/STAR_51bp --readFilesIn sample.fastq

The command should complete within a minute or two. When it's done, explore your directory: ls -al
What do you see?

Once you explore the files, remember to remove them. They are simply a test result. You don't really need to keep them, unless you want to compare them to some other results generated with different parameters.

# clean up
rm *.log Aligned.out.sam SJ.out.tab

The command above is good for aligning one file. But what if you have many files? If you feel ambitious, you can write your command in a script so that it gets executed over and over again and is customize to run on each input file.

Writing a script can be as simple as putting your commands into a file, one after the other. You can then run the script like a program, and your commands will all be executed in order.

However, putting things into a script opens up a new world of possibilities. Here is an example of a script that runs STAR for all the gd7 samples:

#!/bin/bash

# define variables
index=/home/cws/CompGenomics/dm6/Ens98_STAR_51
# get our data files
FILES=/home/cws/CompGenomics/Data/RNA-Seq/gd7_*.fastq.gz

for f in $FILES
do
    echo $f
    base=$(basename $f .fastq.gz)
    echo $base
    STAR --runThreadN 3 --genomeDir $index --readFilesIn $f --outSAMtype BAM SortedByCoordinate \
         --quantMode GeneCounts --readFilesCommand zcat --outFileNamePrefix $base"_"
done

echo "done!"

Some things to notice:

• The first line is special. It tells the computer what program will understand and execute your code. the bash program is the same program that gives you your command line.
• We can use a for loop to do something over and over again, in this case, do something for each input data file.
• We use a program called basename to extract pieces of our input file name, and assign them to a new variable.
• We use echo to check that our variables contain the values we think they do.
• You can build this command up, line by line. You can run it by commenting out the STAR part to make sure everything makes sense before running that part.
• You can also place echo in front of the STAR command, which will print out the command, rather than run it. This way you can check to see if it looks like it's doing the right thing.

At the end of the process above, you should have 2 sets of files for the gd7 mutants. You will have sorted BAM files that can be used for further analysis or visualization, and you will have files of gene counts that you can use for differential expression analysis.

To complete your data set, you need the toll10b data as well. The toll10b files have already been generated for you. Copy them from the class directory to your project directory:

cp /home/cws/CompGenomics/Data/RNA-Seq/STAR/toll10b_*ReadsPerGene.out.tab ./

With a full data set, now you're ready to take the next step of consolidating the data to a single table, and running EdgeR.

Running EdgeR