Alignment to Read Counts & Visualization in IGV

1. Initial Setup
2. Mapping vs Assembly
3. Aligners/Mappers
4. Mapping against the genome vs transcriptome
5. Counting reads as a proxy for gene expression
6. Alignments
7. Running STAR on the experiment
8. Quality Assurance - Mapping statistics as QA/QC.

Initial Setup

This document assumes preproc htstream has been completed. To catch up to where we are:

mkdir -p /share/workshop/mrnaseq_workshop/$USER/rnaseq_example
cd /share/workshop/mrnaseq_workshop/$USER/rnaseq_example
wget https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2022-August-RNA-Seq-Analysis/master/software_scripts/scripts/star_index.slurm
mkdir References
[cp -r /share/workshop/mrnaseq_workshop/jli/rnaseq_example/References/star.overlap100.gencode.M29 References/]
cp -r /share/workshop/mrnaseq_workshop/jli/rnaseq_example/HTS_testing .
ln -s /share/workshop/mrnaseq_workshop/jli/rnaseq_example/01-HTS_Preproc .

Mapping vs Assembly

Assembly seeks to put together the puzzle without knowing what the picture is.

• The focus is on the pieces, how they fit together.

Mapping (or alignment to a reference) tries to put together the puzzle pieces directly onto an image of the picture._

• The focus is on the puzzle, regions of the puzzle that contain certain characteristics (ex. what background) that will help you place the piece onto the puzzle.
• In mapping the question is more, given a small chunk of sequence, where in the genome did this sequence most likely come from.
• The goal then is to find the match(es) with either the “best” edit distance (smallest difference), or all matches with edit distance less than max edit distance. Main issues are:
• Large search space
• Regions of similarity (aka repeats)
• Gaps (INDELS)
• Complexity (RNA, splicing, transcripts)

Alignment concepts

• Multimappers:
  • Reads that align equally well to more than one reference location.
  • Generally, multimappers are discounted in variant detection, and are often discounted in counting applications (like RNA-Seq … would “cancel” out anyway).
  • Note: multimapper “rescue” in some algorithms (RSEM, eXpress).
• Duplicates:
  • Reads or read pairs arising from the same original library fragment, either during library preparation (PCR duplicates) or during sequencing (optical duplicates).
  • Generally, duplicates can only be detected reliably with paired-end sequencing. If PE, they’re discounted in variant detection, and discounted in counting applications (like RNA-Seq).
• Clipping vs Splicing
  • soft-clipped: bases in 5’ and 3’ of the read are NOT part of the alignment.
  • hard-clipped: bases in 5’ and 3’ of the read are NOT part of the alignment AND those bases have been removed from the read sequence in the BAM file. The ‘real’ sequence length would be length(SEQ)+ count-of-hard-clipped-bases
  • From biostars

• Inner length, insert size, fragment length
  
  From This Biostars answer

Considerations when mapping

• Placing reads in regions that do not exist in the reference genome (reads extend off the end of linearized fragments) [ mitochondrial, plasmids, structural variants, etc.].
• Sequencing errors and genetic variation: alignment between read and true source in genome may have more differences than alignment with some other copy of repeat.
• What if the closest fully sequenced genome is too divergent?
• Placing reads in repetitive regions: Some algorithms only return 1 mapping; If multiple: map quality = 0
• Algorithms that use paired-end information => might prefer correct distance over correct alignment.

In RNAseq data, you must also consider effect of splice junctions, reads may span an intron.

Aligners/Mappers

Many alignment algorithms to choose from. Examples include:

• Spliced Aligners
  • STAR
  • HiSAT2 (formerly Tophat [Bowtie2])
  • GMAP/GSNAP
  • SOAPsplice
  • MapSplice
• Aligners that can ’clip’
  • bwa-mem
  • Bowtie2 in local mode

Pseudo-aligners (salmon and kalisto)

• Quasi-mapping
• Probabilistic
• Map to transcripts, not genome
• Does transcript quantifications (or gene)
• Blazing FAST and can run on most laptops
• Experience suggests differences between “traditional” mappers are in the low abundance genes.

Mapping against the genome vs transcriptome

May seem intuitive to map RNAseq data to transcriptome, but it is not that simple.

• Transcriptomes are rarely complete.
• Which transcript of a gene should you map to? canonical transcript (which is that)?
• Shouldn’t map to all splice variants as these would show up as ‘multi-mappers’.

More so, a aligner will try to map every read, somewhere, provided the alignment meets its minimum requirements.

• Need to provide a mapper with all possible places the read could have arisen from, which is best represented by the genome. Otherwise, you get mis-mapping because its close enough.

Genome and Genome Annotation

Genome sequence fasta files and annotation (gff, gtf) files go together! These should be identified at the beginning of analysis.

• Genome fasta files should include all primary chromosomes, unplaced sequences and un-localized sequences, as well as any organelles. Should not contain any contigs that represent patches, or alternative haplotypes.
• If you expect contamination, or the presence of additional sequence/genome, add the sequence(s) to the genome fasta file.
• Annotation file should be GTF (preferred), and should be the most comprehensive you can find.
  • Chromosome names in the GTF must match those in the fasta file (they don’t always do).

Counting reads as a proxy for gene expression

The more you can count (and HTS sequencing systems can count a lot) the better the measure of copy number for even rare transcripts in a population.

• Most RNA-seq techniques deal with count data. Reads are mapped to a reference genome, transcripts are detected, and the number of reads that map to a transcript (or gene) are counted.
• Read counts for a transcript are roughly proportional to the gene’s length and transcript abundance (in whole transcript methods).

Technical artifacts should be considered during counting

• Mapping quality
• Map-ability (uniqueness), the read is not ambiguous

Options are (STAR, HTSEQ, featureCounts)

The HTSEQ way

• Given a sam/bam file with aligned sequence reads and a list of genomic feature (genes locations), we wish to count the number of reads (fragments) than overlap each feature.
  • Features are defined by intervals, they have a start and stop position on a chromosome.
  • For this workshop and analysis, features are genes which are the union of all its exons. You could consider each exon as a feature, for alternative splicing.
• Htseq-count has three overlapping modes
  • union
  • intersection-strict
  • intersection-nonempty

  from the HTSeq Paper

Star Implementation

Counts coincide with Htseq-counts under default parameters (union and tests all orientations). Need to specify GTF file at genome generation step or during mapping.

• Output, 4 columns
  • GeneID
  • Counts for unstranded
  • Counts for first read strand
  • Counts for second read strand

Choose the appropriate column given the library preparation characteristics and generate a matrix expression table, columns are samples, rows are genes.

What does stranded and unstranded mean? Which is better and why? Stranded vs Unstraded

Alignments

1. We are now ready to try an alignment:

    cd /share/workshop/mrnaseq_workshop/$USER/rnaseq_example/HTS_testing
   

   and let’s run STAR (via srun) on the pair of streamed test files we created earlier:

    srun --time=15:00:00 -n 8 --mem=32g --reservation=mrnaseq_workshop --account=workshop --pty /bin/bash
   

   Once you’ve been given an interactive session we can run STAR. You can ignore the two warnings/errors and you know your on a cluster node because your server will change. Here you see I’m on tadpole, then after the srun command is successful, I am now on drove-13.

   jli@ganesh:/share/workshop/mrnaseq_workshop/jli/rnaseq_example/HTS_testing$ srun --time=15:00:00 -n 8 --mem=32g --reservation=mrnaseq_workshop --account=workshop --pty /bin/bash srun: error: spank-auks: cred forwarding failed : auks api : connection failed srun: job 49856359 queued and waiting for resources srun: job 49856359 has been allocated resources bash: /home/jli/.bashrc: Permission denied jli@fleet-29:/share/workshop/mrnaseq_workshop/jli/rnaseq_example/HTS_testing$

2. Then run the star commands

    module load star/2.7.10a
    STAR \
    --runThreadN 12 \
       --genomeDir ../References/star.overlap100.gencode.M29 \
       --outSAMtype BAM SortedByCoordinate \
       --quantMode GeneCounts \
       --outFileNamePrefix mouse_110_WT_C.htstream_ \
       --readFilesCommand zcat \
       --readFilesIn mouse_110_WT_C.htstream_R1.fastq.gz mouse_110_WT_C.htstream_R2.fastq.gz
   

   In the command, we are telling star to count reads on a gene level (‘–quantMode GeneCounts’), the prefix for all the output files will be mouse_110WT_C.htstream, the command to unzip the files (zcat), and finally, the input file pair.

   Once finished please ‘exit’ the srun session. You’ll know you were successful when your back on tadpole

   jli@drove-13:/share/workshop/jli/rnaseq_example/HTS_testing$ exit exit jli@tadpole:/share/workshop/jli/rnaseq_example/HTS_testing$

Now let’s take a look at an alignment in IGV.

1. We first need to index the bam file, will use ‘samtools’ for this step, which is a program to manipulate SAM/BAM files. Take a look at the options for samtools and ‘samtools index’.

   module load samtools
   samtools
   samtools index
   

   We need to index the BAM file:

   cd /share/workshop/mrnaseq_workshop/$USER/rnaseq_example/HTS_testing
   samtools index mouse_110_WT_C.htstream_Aligned.sortedByCoord.out.bam
   

   IF for some reason it didn’t finish, is corrupted or you missed the session, you can copy over a completed copy.

   cp /share/biocore/workshops/2022_mRNASeq/HTS_testing/mouse_110_WT_C.htstream_Aligned.sortedByCoord.out.bam* /share/workshop/mrnaseq_workshop/$USER/rnaseq_example/HTS_testing
   

2. Transfer mouse_110_WT_C.htstream_Aligned.sortedByCoord.out.bam and mouse_110_WT_C.htstream_Aligned.sortedByCoord.out.bam.bai (the index file) to your computer using scp,FileZilla or winSCP.

   Windows users can use WinSCP or FileZilla, both of which are GUI based.

   Mac/Linux users can use scp. In a new shell session on my laptop. NOT logged into tadpole. Replace [your_username] with your username.

    mkdir ~/rnaseq_workshop
    cd ~/rnaseq_workshop
    scp [your_username]@tadpole.genomecenter.ucdavis.edu:/share/workshop/mrnaseq_workshop/[your_username]/rnaseq_example/HTS_testing/mouse_110_WT_C.htstream_Aligned.sortedByCoord.out.bam* .
   

   Its ok of the mkdir command fails (“File exists”) because we aleady created the directory earlier.

3. Now we are ready to use IGV.

   Go to the IGV page at the Broad Institute.

   

   And then navigate to the download page, IGV download

   

   Here you can download IGV for your respective platform (Window, Mac OSX, Linux), but we are going to use the web application they supply, IGV web app.

   

4. The first thing we want to do is load the Human genome. Click on “Genomes” in the menu and choose “Human (GRCm38/mm10)”.

   

5. Now let’s load the alignment bam and index files. Click on “Tracks” and choose “Local File …”.

   

   Navigate to where you transferred the bam and index file and select them both.

   

   Now your alignment is loaded. Any loaded bam file aligned to a genome is called a “track”.

   

6. Lets take a look at the alignment associated with the gene Fn1, and if for some reason it doesn’t find HBB (web IGV can be fickle) go to position chr1:71,610,633-71,628,073. If you don’t see any reads, this likely means your in the wrong genome, double check that it says mm10 in the top left.

   

   

   You can zoom in by clicking on the plus sign (top right) or zoom out by clicking the negative sign (click it twice). You also may have to move around by clicking and dragging in the BAM track window.

   You can also zoom in by clicking and dragging across the number line at the top. That section will highlight, and when you release the button, it will zoom into that section.

   

   

   Reset the window by searching for Fn1 again, and then zoom in 2 steps.

   

7. See that the reads should be aligning within the exons in the gene. This makes sense, since RNA-Seq reads are from exons. Play with the settings on the right hand side a bit and selecting reads.

   

   

   

   

   

Running STAR on the experiment

Start Group Exercise:

• work through Running STAR on all samples
• work through QA/QC of the experiment
• complete the questions at the end

1. We can now run STAR across all samples on the real data using a SLURM script, star.slurm, that we should take a look at now.

    cd /share/workshop/mrnaseq_workshop/$USER/rnaseq_example  # We'll run this from the main directory
    wget https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2022-August-RNA-Seq-Analysis/master/software_scripts/scripts/star.slurm
    less star.slurm
   

   #!/bin/bash
   
    #SBATCH --job-name=star # Job name
    #SBATCH --nodes=1
    #SBATCH --ntasks=8
    #SBATCH --time=60
    #SBATCH --mem=32000 # Memory pool for all cores (see also --mem-per-cpu)
    #SBATCH --partition=production
    #SBATCH --reservation=mrnaseq_workshop
    #SBATCH --account=workshop
    #SBATCH --array=1-22
    #SBATCH --output=slurmout/star_%A_%a.out # File to which STDOUT will be written
    #SBATCH --error=slurmout/star_%A_%a.err # File to which STDERR will be written
   
    start=`date +%s`
    echo $HOSTNAME
    echo "My SLURM_ARRAY_TASK_ID: " $SLURM_ARRAY_TASK_ID
   
    sample=`sed "${SLURM_ARRAY_TASK_ID}q;d" samples.txt`
    REF="References/star.overlap100.gencode.M29"
   
    outpath='02-STAR_alignment'
    [[ -d ${outpath} ]] || mkdir ${outpath}
    [[ -d ${outpath}/${sample} ]] || mkdir ${outpath}/${sample}
   
    echo "SAMPLE: ${sample}"
   
    module load star
   
    call="STAR
         --runThreadN ${SLURM_NTASKS} \
         --genomeDir $REF \
         --outSAMtype BAM SortedByCoordinate \
         --readFilesCommand zcat \
         --readFilesIn 01-HTS_Preproc/${sample}/${sample}_R1.fastq.gz 01-HTS_Preproc/${sample}/${sample}_R2.fastq.gz \
         --quantMode GeneCounts \
         --outFileNamePrefix ${outpath}/${sample}/${sample}_ \
         > ${outpath}/${sample}/${sample}-STAR.stdout 2> ${outpath}/${sample}/${sample}-STAR.stderr"
   
    echo $call
    eval $call
   
    end=`date +%s`
    runtime=$((end-start))
    echo $runtime
    

When you are done, type “q” to exit.

2. After looking at the script, lets run it.

    sbatch star.slurm  # moment of truth!
   

   We can watch the progress of our task array using the ‘squeue’ command. Takes about 30 minutes to process each sample.

    squeue -u $USER 
   

Quality Assurance - Mapping statistics as QA/QC.

1. Once your jobs have finished successfully, check the error and out logs like we did in the previous exercise.

   Use a script of ours, star_stats.sh to collect the alignment stats. Don’t worry about the script’s contents at the moment; you’ll use very similar commands to create a counts table in the next section. For now:

    cd /share/workshop/mrnaseq_workshop/$USER/rnaseq_example # We'll run this from the main directory
    wget https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2022-August-RNA-Seq-Analysis/master/software_scripts/scripts/star_stats.sh
    bash star_stats.sh
   

   #!/bin/bash
   
    echo -en "sample_names" > names.txt
    echo -en "total_in_feature\t" > totals.txt
    cat 02-STAR_alignment/*/*ReadsPerGene.out.tab | head -4 | cut -f1 > stats.txt
    cat samples.txt | while read sample; do
        echo ${sample}
        echo -en "\t${sample}" >> names.txt
        head -4 02-STAR_alignment/${sample}/${sample}_ReadsPerGene.out.tab | cut -f4 > temp1
        paste stats.txt temp1 > temp2
        mv temp2 stats.txt
        tail -n +5 02-STAR_alignment/${sample}/${sample}_ReadsPerGene.out.tab | cut -f4 | \
            perl -ne '$tot+=$_ }{ print "$tot\t"' >> totals.txt
    done
    echo -en "\n" >> names.txt
    cat names.txt stats.txt totals.txt > temp1
    mv temp1 summary_star_alignments.txt
    rm stats.txt
    rm names.txt
    rm totals.txt
    

2. Transfer summary_star_alignments.txt to your computer using scp or winSCP, or copy/paste from cat [sometimes doesn’t work],

   In Mac/Linux, users can use scp. Windows users can use WinSCP. In a new shell session on my laptop. NOT logged into tadpole. Replace my [your_username] with your username.

    mkdir -p ~/rnaseq_workshop
    cd ~/rnaseq_workshop
    scp [your_username]@tadpole.genomecenter.ucdavis.edu:/share/workshop/mrnaseq_workshop/[your_username]/rnaseq_example/summary_star_alignments.txt .
   

Questions:

1. Look at the script star.slurm. What does the array=1-22 mean, why is it used, and what is the usage of it in the script itself?
2. Look through the files in an output directory and check out what is present and discuss what each of them mean. (for example: cd /share/workshop/mrnaseq_workshop/$USER/rnaseq_example/02-STAR_alignment/mouse_110_WT_C )
3. Come up with a brief command you might use to check that all of the sample alignments using STAR have a reasonable output and/or did not produce any errors.
4. Open summary_star_alignments.txt in excel (or excel like application), and review. The table that this script creates (“summary_star_alignments.txt”) can be pulled to your laptop via ‘scp’, or WinSCP, etc., and imported into a spreadsheet. Are all samples behaving similarly? Discuss …
5. If time, find some other regions/genes with high expression using IGV with your group. (Looking at genes the paper references is a great place to start)

Stop Group Exercise