In this tutorial we will use a splice aware aligner HISAT2 that is the best solution when using a laptop or a personal computer with less than 16GB. If you have a computer or access to computer cluster it is worthwhile to consider using STAR and compare the results.

1 HISAT2

1.1 Data

All FASTQ files that you will need for this exercise can be found in the fastq subfolder of your data folder.

A pre-built mouse genome index for HISAT2 can be found in the references/hisat2_build subfolder of your data folder.

1.2 Mapping

Here, you will map the reads to the mouse reference genome using the RNA-seq aligner HISAT2.

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.

Now 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.

sh 
mkdir outDir

# for paired end reads
hisat2 -p N -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.txt


# for single end reads
hisat2 -p N  -x path/to/index/fileName \
  -U path/to/reads/sample.fastq \
  -S outDir/hisat2.sam \
  --summary-file outDir/hisat2_summary.txt

The flags used are:

  • -p N specifies the number of threads that will be used by the program
  • -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 file of paired end reads that you wish to map
  • -2 /path/to/reads/sample_2.fastq is where you should list the second-read FASTQ file of paired end reads that you wish to map
  • -U /path/to/reads/sample.fastq is where you should list the single end read FASTQ file 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 alignments

This 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:

  • How many RNA-seq read pairs were provided as input to HISAT2?
  • How many of those read pairs were mapped by HISAT2?
  • How many reads were uniquely mapped, i.e. mapped to one genomic location?
  • In general, do the alignments seem to be good? I.e. do they cover the entire reads and contain few mismatches?

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.

2 Convert SAM to BAM

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.mmusculus.HISAT2.bam.

The most commonly used tool for converting from SAM to BAM is Samtools (follow the link for more information about Samtools).

The Samtools command to convert from SAM to a sorted BAM file is:

sh 
samtools sort -o output.bam input.sam

Remember to use an appropriate filename instead of output.bam! Next, we need to index the BAM file.

sh 
samtools index properName.bam

The 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:

sh 
samtools flagstat properName.sorted.bam > properName.sorted.flagstat

Since 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.

sh 
rm input.sam

3 Code to run all the steps above on data on course data

sh 



###############################################################################
###          MAP READS TO REFERENCE                                       #####
###############################################################################


######################################
###  BUILD A HISAT2 REFERENCE     ####
######################################





###########################
###  ALIGN WITH HISAT2  ###
###########################
mkdir hisat2_results

#map all the small files
 for i in $(ls fastq/*.clean.fq.gz.1M_reads.fq.gz); do

    file=${i#fastq/}
    sample=${file%.clean.fq.gz.1M_reads.fq.gz}
   #Single-end read
   hisat2 -p 4 \
    -x references/hisat2_build/mm10 \
    -U $i \
    -S hisat2_results/$sample.1M.sam \
    --summary-file hisat2_results/$sample.hisat2_summary.txt


 done


#map  the large file
 for i in $(ls fastq/*.clean.fq.gz); do

    file=${i#fastq/}
    sample=${file%.clean.fq.gz}
   #Single-end read
   hisat2 -p 4 \
    -x references/hisat2_build/mm10 \
    -U $i \
    -S hisat2_results/$sample.10M.sam \
    --summary-file hisat2_results/$sample.10M.hisat2_summary.txt



 done




############################
####   CONVERT TO BAM   ####

############################


for i in $(ls hisat2_results/*.sam); do

    file=${i#hisat2_results/}
    sample=${file%.sam}
 
 
 # This is how it was done in a older version of samtools  (Version: 0.1.19-44428cd)
  #samtools view -bS -o hisat2_results/$sample.bam hisat2_results/$sample.sam 
  #samtools sort  hisat2_results/$sample.bam hisat2_results/$sample.sorted
  
 # This is how it should be done in the version that is used in our conda enviroment (Version: 1.10 (using htslib 1.10.2))
  samtools sort -o hisat2_results/$sample.sorted.bam hisat2_results/$sample.sam
  samtools index hisat2_results/$sample.sorted.bam
  samtools flagstat hisat2_results/$sample.sorted.bam > $sample.sorted.flagstat

done

############################
####   REMOVE SAM FILES ####
############################


# Remove the intermediate sam and bamfiles.
for i in $(ls hisat2_results/*.sorted.bam); do

    file=${i#hisat2_results/}
    sample=${file%.sorted.bam}

 if test -f "hisat2_results/$sample.sorted.bam"; then
    echo "$sample.sorted.bam exists. Removing $sample.sam and $sample.bam"
    rm hisat2_results/$sample.sam

    # This intermediate step is no longer needed anymore 
    #rm hisat2_results/$sample.bam
 fi

done