FastQ reads are mapped to a reference genome using splice-aware RNA-seq aligner STAR.
All FASTQ files that you will need for this exercise can be found in /proj/uppstore2017171/courses/RNAseqWorkshop/downloads/isoform/referenceBased/data on UPPMAX.
If you want to map more files for practice, you can continue with files for additional samples, found in /proj/uppstore2017171/courses/RNAseqWorkshop/downloads/isoform/RAB11FIP5_fastqFiles on UPPMAX.
A pre-built human genome index for STAR is found here /proj/uppstore2017171/courses/RNAseqWorkshop/downloads/reference/hg19_Gencode14.overhang75 on UPPMAX.
Here we will map the reads to the hg19 reference genome using a popular RNA-seq aligner, STAR. There are many many features that can be tweaked using STAR. For more information concerning different features that can be used see the manual.
Read below for the flags we use for this exercise. Remember to change filenames accordingly so that you can run your program properly and know which files you have used.
To load the STAR module on UPPMAX, execute
module load bioinfo-tools
module load module load star/2.5.3aNow you can map the reads from one of the samples (or several; it’s up to you which ones) using a command such as the one below.
mkdir outDir
STAR --runThreadN N --outSAMstrandField intronMotif --genomeDir /path/to/index \
--readFilesIn /path/to/reads/sample_1.fastq /path/to/reads/sample_2.fastq \
--outFileNamePrefix outDir/Flags used are
--runThreadN N specifies the number of threads that will be used by the program.--outSAMstrandField intronMotif adds information (to the SAM output file) required for downstream analysis with Cufflinks--genomeDir /path/to/index specifies the directory containing the pre-built genome index--readFilesIn /path/to/reads/sample_1.fastq /path/to/reads/sample_2.fastq is where you should list the FASTQ files that you wish to map--outFileNamePrefix outDir specifies the output directoryThis should run fairly quickly and put a file called Aligned.out.sam in the directory that you specified with --outFileNamePrefix.
Try to answer the following:
All FASTQ files that you will need for this exercise can be found in /proj/uppstore2017171/courses/RNAseqWorkshop/downloads/isoform/referenceBased/data on UPPMAX.
A pre-built human genome index for HISAT2 is found here /proj/uppstore2017171/courses/RNAseqWorkshop/downloads/reference/hg19_hisat2 on UPPMAX.
Here, you will map the reads to the hg19 reference genome using the RNA-seq aligner HISAT2. Note that if you are using your own non-human data, you need to use a reference genome for the corresponding species.
There are many features that can be tweaked using HISAT2. For more information on all flags that can be used go here.
Read below for the flags we use for this exercise. Remember to change filenames accordingly so that you can run your program properly and know which files you have used.
To load the HISAT2 module on UPPMAX, execute:
module load bioinfo-tools # This is to get access to all bioinformatics tools available on UPPMAX
module load HISAT2/2.1.0 # The specific mapping program
mkdir outDir
hisat2 -p N --dta-cufflinks -x path/to/index/fileName \
-1 path/to/reads/sample_1.fastq -2 path/to/reads/sample_2.fastq \
-S outDir/hisat2.sam --summary-file outDir/hisat2_summary.txtThe flags used are:
-p N specifies the number of threads that will be used by the program--dta-cufflinks will generate output that is optimal for downstream analysis with Cufflinks-x /path/to/index/fileName specifies the path to the pre-built genome index. Note that the index consists of multiple files ending in .ht2, and only the shared part of the filename should be indicated (e.g. genome if the files are called genome.1.ht2, genome.2.ht2, etc).-1 /path/to/reads/sample_1.fastq is where you should list the first-read FASTQ files that you wish to map-2 /path/to/reads/sample_2.fastq is where you should list the second-read FASTQ files that you wish to map-S outDir/hisat2.sam is the name of the result file that will be created--summary-file outDir/hisat2_summary.txt is the name of a file for summary information about the alignmentsThis should run fairly quickly and create the files you specified with -S and --summary-file.
If everything worked, HISAT2 should report some statistics about how many reads were mapped, on your terminal and in the summary file.
Try to answer the following:
To answer these questions, you should look at the input to and output from HISAT2. You may also need to consult the HISAT2 manual, information about the FASTQ format and the SAM format specification.
If you were able to run HISAT2 and STAR successfully, this should have produced files with mapped reads in SAM format. These files need to be converted to sorted and indexed BAM files for efficient downstream analysis.
You should try to give the BAM files representable names, in order to make it easier to manage your files. A good naming scheme for BAM files is to use names that indicate what you mapped and how. As an example, if you mapped sample 12 using HISAT2 you could create a file named sample12_RAB11FIP5.hg19.HISAT2.bam.
The most commonly used tool for converting from SAM to BAM is Samtools (follow the link for more information about Samtools).
To load the Samtools module on UPPMAX, execute:
module load bioinfo-tools
module load samtools/1.9The Samtools command to convert from SAM to a sorted BAM file is:
samtools sort -o output.bam input.samRemember to use an appropriate filename instead of output.bam! Next, we need to index the BAM file.
samtools index properName.bamThe indexing step should create an index file with the suffix .bai. This sorted, indexed BAM file can be viewed in the Integrative Genomics Viewer (IGV).
You can also get a report on your mapped reads using the samtools command flagstat:
samtools flagstat properName.sorted.bam > properName.sorted.flagstatSince the BAM file contains all the information from the original SAM file, remember to remove the SAM file once you are finished, in order to free up disk space.
End of document