MiCM: 01 - Alignment

This workshop will guide you through the basics of ChIP-seq analysis with hands-on exercises. You will learn how to process ChIP-seq data: perform read alignment, peak calling, quality control, visualization through the genome browser, motif finding and gene set enrichment analysis.

In this tutorial, we are going to focus on performing read alignment.

The dataset that we will be using today is ChIP-seq for HIF1a and HIF2a. HIF1a and HIF2a are the two principal isoforms of the Hypoxia-inducible factor (HIF), which is the major transcriptional regulator of cellular responses to hypoxia. This dataset was generated by Smythies et al (2018), and you can find the paper here. Because the original data is very big, we are going to work with a subset of the original reads from the study (reads that map to a segment of the chromosome 6).

Note: All the code in this tutorial will be performed in the terminal.

Let’s move to the scripts directory (or stay here if already there):

cd chip_seq/scripts

01 - Mapping with bowtie2

You can check the help on how to use bowtie2 by typing:

bowtie2 --help 

Bowtie 2 version 2.3.5.1 by Ben Langmead (langmea@cs.jhu.edu, www.cs.jhu.edu/~langmea)
Usage: 
  bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r> | --interleaved <i> | -b <bam>} [-S <sam>]

  <bt2-idx>  Index filename prefix (minus trailing .X.bt2).
             NOTE: Bowtie 1 and Bowtie 2 indexes are not compatible.
  <m1>       Files with #1 mates, paired with files in <m2>.
             Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
  <m2>       Files with #2 mates, paired with files in <m1>.
             Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
  <r>        Files with unpaired reads.
             Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
  <i>        Files with interleaved paired-end FASTQ/FASTA reads
             Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
  <bam>      Files are unaligned BAM sorted by read name.
  <sam>      File for SAM output (default: stdout)

  <m1>, <m2>, <r> can be comma-separated lists (no whitespace) and can be
  specified many times.  E.g. '-U file1.fq,file2.fq -U file3.fq'.

Options (defaults in parentheses):

 Input:
  -q                 query input files are FASTQ .fq/.fastq (default)
  --tab5             query input files are TAB5 .tab5
  --tab6             query input files are TAB6 .tab6
  --qseq             query input files are in Illumina's qseq format
  -f                 query input files are (multi-)FASTA .fa/.mfa
  -r                 query input files are raw one-sequence-per-line
  -F k:<int>,i:<int> query input files are continuous FASTA where reads
                     are substrings (k-mers) extracted from a FASTA file <s>
                     and aligned at offsets 1, 1+i, 1+2i ... end of reference
  -c                 <m1>, <m2>, <r> are sequences themselves, not files
  -s/--skip <int>    skip the first <int> reads/pairs in the input (none)
  -u/--upto <int>    stop after first <int> reads/pairs (no limit)
  -5/--trim5 <int>   trim <int> bases from 5'/left end of reads (0)
  -3/--trim3 <int>   trim <int> bases from 3'/right end of reads (0)
  --trim-to [3:|5:]<int> trim reads exceeding <int> bases from either 3' or 5' end
                     If the read end is not specified then it defaults to 3 (0)
  --phred33          qualities are Phred+33 (default)
  --phred64          qualities are Phred+64
  --int-quals        qualities encoded as space-delimited integers

 Presets:                 Same as:
  For --end-to-end:
   --very-fast            -D 5 -R 1 -N 0 -L 22 -i S,0,2.50
   --fast                 -D 10 -R 2 -N 0 -L 22 -i S,0,2.50
   --sensitive            -D 15 -R 2 -N 0 -L 22 -i S,1,1.15 (default)
   --very-sensitive       -D 20 -R 3 -N 0 -L 20 -i S,1,0.50

  For --local:
   --very-fast-local      -D 5 -R 1 -N 0 -L 25 -i S,1,2.00
   --fast-local           -D 10 -R 2 -N 0 -L 22 -i S,1,1.75
   --sensitive-local      -D 15 -R 2 -N 0 -L 20 -i S,1,0.75 (default)
   --very-sensitive-local -D 20 -R 3 -N 0 -L 20 -i S,1,0.50

 Alignment:
  -N <int>           max # mismatches in seed alignment; can be 0 or 1 (0)
  -L <int>           length of seed substrings; must be >3, <32 (22)
  -i <func>          interval between seed substrings w/r/t read len (S,1,1.15)
  --n-ceil <func>    func for max # non-A/C/G/Ts permitted in aln (L,0,0.15)
  --dpad <int>       include <int> extra ref chars on sides of DP table (15)
  --gbar <int>       disallow gaps within <int> nucs of read extremes (4)
  --ignore-quals     treat all quality values as 30 on Phred scale (off)
  --nofw             do not align forward (original) version of read (off)
  --norc             do not align reverse-complement version of read (off)
  --no-1mm-upfront   do not allow 1 mismatch alignments before attempting to
                     scan for the optimal seeded alignments
  --end-to-end       entire read must align; no clipping (on)
   OR
  --local            local alignment; ends might be soft clipped (off)

 Scoring:
  --ma <int>         match bonus (0 for --end-to-end, 2 for --local) 
  --mp <int>         max penalty for mismatch; lower qual = lower penalty (6)
  --np <int>         penalty for non-A/C/G/Ts in read/ref (1)
  --rdg <int>,<int>  read gap open, extend penalties (5,3)
  --rfg <int>,<int>  reference gap open, extend penalties (5,3)
  --score-min <func> min acceptable alignment score w/r/t read length
                     (G,20,8 for local, L,-0.6,-0.6 for end-to-end)

 Reporting:
  (default)          look for multiple alignments, report best, with MAPQ
   OR
  -k <int>           report up to <int> alns per read; MAPQ not meaningful
   OR
  -a/--all           report all alignments; very slow, MAPQ not meaningful

 Effort:
  -D <int>           give up extending after <int> failed extends in a row (15)
  -R <int>           for reads w/ repetitive seeds, try <int> sets of seeds (2)

 Paired-end:
  -I/--minins <int>  minimum fragment length (0)
  -X/--maxins <int>  maximum fragment length (500)
  --fr/--rf/--ff     -1, -2 mates align fw/rev, rev/fw, fw/fw (--fr)
  --no-mixed         suppress unpaired alignments for paired reads
  --no-discordant    suppress discordant alignments for paired reads
  --dovetail         concordant when mates extend past each other
  --no-contain       not concordant when one mate alignment contains other
  --no-overlap       not concordant when mates overlap at all

 BAM:
  --align-paired-reads Bowtie2 will, by default, attempt to align unpaired BAM reads.
                     Use this option to align paired-end reads instead.
  --preserve-tags    Preserve tags from the original BAM record by
                     appending them to the end of the corresponding SAM output.

 Output:
  -t/--time          print wall-clock time taken by search phases
  --un <path>        write unpaired reads that didn't align to <path>
  --al <path>        write unpaired reads that aligned at least once to <path>
  --un-conc <path>   write pairs that didn't align concordantly to <path>
  --al-conc <path>   write pairs that aligned concordantly at least once to <path>
    (Note: for --un, --al, --un-conc, or --al-conc, add '-gz' to the option name, e.g.
    --un-gz <path>, to gzip compress output, or add '-bz2' to bzip2 compress output.)
  --quiet            print nothing to stderr except serious errors
  --met-file <path>  send metrics to file at <path> (off)
  --met-stderr       send metrics to stderr (off)
  --met <int>        report internal counters & metrics every <int> secs (1)
  --no-unal          suppress SAM records for unaligned reads
  --no-head          suppress header lines, i.e. lines starting with @
  --no-sq            suppress @SQ header lines
  --rg-id <text>     set read group id, reflected in @RG line and RG:Z: opt field
  --rg <text>        add <text> ("lab:value") to @RG line of SAM header.
                     Note: @RG line only printed when --rg-id is set.
  --omit-sec-seq     put '*' in SEQ and QUAL fields for secondary alignments.
  --sam-no-qname-trunc Suppress standard behavior of truncating readname at first whitespace 
                      at the expense of generating non-standard SAM.
  --xeq              Use '='/'X', instead of 'M,' to specify matches/mismatches in SAM record.
  --soft-clipped-unmapped-tlen Exclude soft-clipped bases when reporting TLEN

 Performance:
  -p/--threads <int> number of alignment threads to launch (1)
  --reorder          force SAM output order to match order of input reads
  --mm               use memory-mapped I/O for index; many 'bowtie's can share

 Other:
  --qc-filter        filter out reads that are bad according to QSEQ filter
  --seed <int>       seed for random number generator (0)
  --non-deterministic seed rand. gen. arbitrarily instead of using read attributes
  --version          print version information and quit
  -h/--help          print this usage message

You can also access the bowtie2 manual in here.

The data that we have is paired end data, so we have 2 .fastq files for one sample. They look like this:

ls ../data/fastq

HIF1a_Rep1_R1.fastq.gz
HIF1a_Rep1_R2.fastq.gz
HIF2a_Rep1_R1.fastq.gz
HIF2a_Rep1_R2.fastq.gz

Bowtie2 requires a genome index (compressed files that store the reference sequence). For this time, we’ve already generated the index and can be found here:

ls ../data/reference/GRCh38_chr6

hg38_chr6.1.bt2
hg38_chr6.2.bt2
hg38_chr6.3.bt2
hg38_chr6.4.bt2
hg38_chr6.rev.1.bt2
hg38_chr6.rev.2.bt2

This index was created only using the human chromosome 6. You can create your own index or download pre-built ones from the bowtie2 page.

Let’s create a directory where to store our alignment

mkdir -p ../results

Now we are ready to run bowtie2!

bowtie2 --very-sensitive-local -x ../data/reference/GRCh38_chr6/hg38_chr6 -1 ../data/fastq/HIF1a_Rep1_R1.fastq.gz -2 ../data/fastq/HIF1a_Rep1_R2.fastq.gz -S ../results/HIF1a_Rep1.sam

121159 reads; of these:
  121159 (100.00%) were paired; of these:
    1526 (1.26%) aligned concordantly 0 times
    103901 (85.76%) aligned concordantly exactly 1 time
    15732 (12.98%) aligned concordantly >1 times
    ----
    1526 pairs aligned concordantly 0 times; of these:
      928 (60.81%) aligned discordantly 1 time
    ----
    598 pairs aligned 0 times concordantly or discordantly; of these:
      1196 mates make up the pairs; of these:
        306 (25.59%) aligned 0 times
        456 (38.13%) aligned exactly 1 time
        434 (36.29%) aligned >1 times
99.87% overall alignment rate

Let’s take a look at the parameters that we’ve used for running bowtie2:

• –very-sensitive-local: Bowtie2 comes with some useful combinations of parameters packaged into shorter “preset” parameters. The preset options that come with Bowtie 2 are designed to cover a wide area of the speed/sensitivity/accuracy trade-off space. “Fast” presets are generally faster but less sensitive, while the “sensitive” presets are generally slower but more sensitive and accurate. The local part means that we are doing ‘local’ alignment (remember that the default in Bowtie2 is the end-to-end alignment).
• –x: path to the index
• -1: path to fastq1
• -2: path to fastq2
• -S: path to the output SAM file

02 - Filtering SAM/BAM

For the following steps, we will be using samtools. The manual can be found here. Samtools is a compendium of various tools:

samtools --help

Program: samtools (Tools for alignments in the SAM format)
Version: 1.10 (using htslib 1.10.2-3ubuntu0.1)

Usage:   samtools <command> [options]

Commands:
  -- Indexing
     dict           create a sequence dictionary file
     faidx          index/extract FASTA
     fqidx          index/extract FASTQ
     index          index alignment

  -- Editing
     calmd          recalculate MD/NM tags and '=' bases
     fixmate        fix mate information
     reheader       replace BAM header
     targetcut      cut fosmid regions (for fosmid pool only)
     addreplacerg   adds or replaces RG tags
     markdup        mark duplicates

  -- File operations
     collate        shuffle and group alignments by name
     cat            concatenate BAMs
     merge          merge sorted alignments
     mpileup        multi-way pileup
     sort           sort alignment file
     split          splits a file by read group
     quickcheck     quickly check if SAM/BAM/CRAM file appears intact
     fastq          converts a BAM to a FASTQ
     fasta          converts a BAM to a FASTA

  -- Statistics
     bedcov         read depth per BED region
     coverage       alignment depth and percent coverage
     depth          compute the depth
     flagstat       simple stats
     idxstats       BAM index stats
     phase          phase heterozygotes
     stats          generate stats (former bamcheck)

  -- Viewing
     flags          explain BAM flags
     tview          text alignment viewer
     view           SAM<->BAM<->CRAM conversion
     depad          convert padded BAM to unpadded BAM

We can count the number of reads that are in our SAM

samtools view -c ../results/HIF1a_Rep1.sam

242318

• -c: Instead of printing the alignments, only count them and print the total number

We need to convert the SAM file into a BAM file. In this same step, we will perform some filtering of the alignments:

samtools view -h -F 1804 -q 30 -b ../results/HIF1a_Rep1.sam -o ../results/HIF1a_Rep1.bam

• -h: include header in the output
• -F: Exclude alignments with a given FLAG field. The 1804 FLAG excludes the following reads: read unmapped, mate unmapped, not primary alignment, read quality low, PCR / optical duplicate. You can check the meaning of the SAM FLAGS in here. More information about samtools flags can be found here.
• -q: Skip alignments with MAPQ smaller than the given value. The MAPQ is a non-negative integer Q = -10 log10 p, where p is an estimate of the probability that the alignment does not correspond to the read’s true point of origin. Usually, a good value is between 10-30.
• -b: Output in the BAM format.
• -o: Output to file

We can also remove reads that map blacklisted regions by ENCODE using bedtools.

bedtools intersect -abam ../results/HIF1a_Rep1.bam -b ../data/blacklist/hg38.blacklist.bed.gz -v > ../results/HIF1a_Rep1.rmblacklist.bam

Finally, the next steps deal with marking and removing duplicates. Removing duplicates is usually a good idea, since most likely duplicate reads are originated from PCR, and don’t represent ‘true’ reads.

# collate: group by name
samtools collate ../results/HIF1a_Rep1.rmblacklist.bam -o ../results/HIF1a_Rep1.collate.bam
# fills in mate coordinates and insert size fields
samtools fixmate -m ../results/HIF1a_Rep1.collate.bam ../results/HIF1a_Rep1.fixmate.bam
# sort bam
samtools sort ../results/HIF1a_Rep1.fixmate.bam -o ../results/HIF1a_Rep1.sorted.bam
# mark and remove duplicates
samtools markdup ../results/HIF1a_Rep1.sorted.bam ../results/HIF1a_Rep1.sorted.markdup.bam # with this command we don't remove duplicates
samtools markdup -r ../results/HIF1a_Rep1.sorted.bam ../results/HIF1a_Rep1.sorted.rmdup.bam
# create index
samtools index ../results/HIF1a_Rep1.sorted.rmdup.bam

We can compare the initial number of reads in the fastq file versus the number of reads and fragments after mapping and filtering.

## Function to count the number of reads in a fastq.gz file
function nreadsfastq(){
    zcat ${1} | echo $((`wc -l`/4))
}
nreadsfastq ../data/fastq/HIF1a_Rep1_R1.fastq.gz
nreadsfastq ../data/fastq/HIF1a_Rep1_R2.fastq.gz

121159
121159

samtools flagstat ../results/HIF1a_Rep1.sorted.markdup.bam

234228 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
41550 + 0 duplicates
234228 + 0 mapped (100.00% : N/A)
234152 + 0 paired in sequencing
117076 + 0 read1
117076 + 0 read2
231884 + 0 properly paired (99.03% : N/A)
234152 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

samtools flagstat ../results/HIF1a_Rep1.sorted.rmdup.bam

192678 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
192678 + 0 mapped (100.00% : N/A)
192610 + 0 paired in sequencing
96305 + 0 read1
96305 + 0 read2
190460 + 0 properly paired (98.88% : N/A)
192610 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

And that’s it! Now you are familiar with alignment using bowtie2 and SAM/BAM processing.

03 - Exercises

• 1. Perform the alignment steps for the HIF2a data. The input fastq files can be found in the data folder:

ls ../data/fastq/

HIF1a_Rep1_R1.fastq.gz
HIF1a_Rep1_R2.fastq.gz
HIF2a_Rep1_R1.fastq.gz
HIF2a_Rep1_R2.fastq.gz

• 2. How would you filter for reads that are unmapped and have a mate unmapped? Give it a try using the unfiltered .sam output file from HIF1a. Remember that you can check the SAM FLAGS in here.
• 3. How would you run bowtie2 using the end-to-end mode? How do the results compare to the local mode?

sessionInfo()

R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.6 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_CA.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_CA.UTF-8        LC_COLLATE=en_CA.UTF-8    
 [5] LC_MONETARY=en_CA.UTF-8    LC_MESSAGES=en_CA.UTF-8   
 [7] LC_PAPER=en_CA.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C       

time zone: America/New_York
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] digest_0.6.33   R6_2.5.1        fastmap_1.1.1   xfun_0.40      
 [5] cachem_1.0.8    knitr_1.43      htmltools_0.5.6 rmarkdown_2.24 
 [9] cli_3.6.1       sass_0.4.7      jquerylib_0.1.4 compiler_4.3.1 
[13] tools_4.3.1     mime_0.12       evaluate_0.21   bslib_0.5.1    
[17] yaml_2.3.7      rlang_1.1.1     jsonlite_1.8.7