PostAlignment.Rmd

---
title: "ATAC-Seq Data Analysis with Human Data - Perform post-alignment processing"
author: "Anni Liu"
date: '`r format(Sys.Date(), "%B %d, %Y")`'
output: 
  html_document:
    code_folding: show
---

```{r, shorcut, include=FALSE}
## RStudio keyboard shortcut
# Cursor at the beginning of a command line: Ctrl+A
# Cursor at the end of a command line: Ctrl+E
# Clear all the code from your console: Ctrl+L
# Create a pipe operator %>%: Ctrl+Shift+M (Windows) or Cmd+Shift+M (Mac)
# Create an assignment operator <-: Alt+- (Windows) or Option+- (Mac) 
# Knit a document (knitr): Ctrl+Shift+K (Windows) or Cmd+Shift+K (Mac)
# Comment or uncomment current selection: Ctrl+Shift+C (Windows) or Cmd+Shift+C (Mac)
```


# [QC] Perform post-alignment processing
## Select properly paired reads using `ScanBamParam()` and `scanBamFlag()`
```{r}
library(GenomicAlignments)
flags <- scanBamFlag(isProperPair = T)
# Here we only select reads paired in alignment within the preset max fragment length 600 bp
# isProperPair = T: represent reads aligning to identical reference sequences and with a specified distance, 600 bp in our case.

# Recap maxFragLength = 600 # Line 501: Perform the alignment relative to hg38 using the new ATACseq data from female lung
# align(index = "BSgenome.Hsapiens.UCSC.hg38.mainchr",
#       readfile1 = "./data/SRR6870408_1.fastq.gz", 
#       readfile2 = "./data/SRR6870408_2.fastq.gz",
#       output_file = "./data/ATAC_female_lung.bam",
#       nthreads = 4, type = 1,
#       unique = TRUE, maxFragLength = 600)

# Set up parameters for scanning BAM files
chr20_param  <- ScanBamParam(flag = flags,
                             what = c("qname", "mapq", "isize"),
                             which = GRanges("chr20", IRanges(1, 64444167))) 

# https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&chromInfoPage=
# 63025520: length of chromosome 20 (GRCh37/hg19)
# https://genome.ucsc.edu/cgi-bin/hgTracks?hgsid=1795230954_XcEyjdAa9Z968dnsdwUvHaqjwvOg&chromInfoPage=
# 64444167: length of chromosome 20 (GRCh38/hg38)
# qname: reads name
# mapq: quality of alignment
# isize: insert size of our fragment
chr20_param

# class: ScanBamParam
# bamFlag (NA unless specified): isProperPair=TRUE
# bamSimpleCigar: FALSE
# bamReverseComplement: FALSE
# bamTag:  
# bamTagFilter:
# bamWhich: 1 ranges
# bamWhat: qname, mapq, isize
# bamMapqFilter: NA
```

```{r two.chromosomes}
chr1X_param <- ScanBamParam(flag = flags,
                            what = c("qname", "mapq", "isize"),
                            which = IRangesList(chr1=IRanges(1, 248956422),
                                                chrX=IRanges(1, 156040895))) 
```

```{r all.chromosomes}
# Approach 1:
# https://genome.ucsc.edu/cgi-bin/hgTracks?hgsid=1795230954_XcEyjdAa9Z968dnsdwUvHaqjwvOg&chromInfoPage=
scanBamHeader(sorted_bam)[[1]]$targets
mainchr_param <- ScanBamParam(flag = flags,
                              what = c("qname", "mapq", "isize"),
                              which = IRangesList(chr1=IRanges(1, 248956422),
                                                  chr2=IRanges(1, 242193529),
                                                  chr3=IRanges(1, 198295559),
                                                  chr4=IRanges(1, 190214555),
                                                  chr5=IRanges(1, 181538259),
                                                  chr6=IRanges(1, 170805979),
                                                  chr7=IRanges(1, 159345973),
                                                  chr8=IRanges(1, 145138636),
                                                  chr9=IRanges(1, 138394717),
                                                  chr10=IRanges(1, 133797422),
                                                  chr11=IRanges(1, 135086622),
                                                  chr12=IRanges(1, 133275309),
                                                  chr13=IRanges(1, 114364328),
                                                  chr14=IRanges(1, 107043718),
                                                  chr15=IRanges(1, 101991189),
                                                  chr16=IRanges(1, 90338345),
                                                  chr17=IRanges(1, 83257441),
                                                  chr18=IRanges(1, 80373285),
                                                  chr19=IRanges(1, 58617616),
                                                  chr20=IRanges(1, 64444167),
                                                  chr21=IRanges(1, 46709983),
                                                  chr22=IRanges(1, 50818468),
                                                  chrX=IRanges(1, 156040895),
                                                  chrY=IRanges(1, 57227415))) 


# Approach 2:
as(Seqinfo(genome = "hg38"), "GRanges")@seqnames@values
mainchr_param_2 <- ScanBamParam(flag = flags,
                                what = c("qname", "mapq", "isize"),
                                which = as(Seqinfo(genome="hg38"), "GRanges")[1:24]) 
mainchr_param_2@which@partitioning@NAMES <- paste0("chr", c(1:22, "X", "Y"))
mainchr_param_2@which@partitioning@end <- 1:24
```


## Read properly paired reads
```{r}
# Global view of the 2 reads, containing the seqnames, strand, ranges of 2 reads
atac_reads_chr20 <- readGAlignmentPairs(sorted_bam, param = chr20_param)
class(atac_reads_chr20)
# [1] "GAlignmentPairs"
# attr(,"package")
# [1] "GenomicAlignments"
atac_reads_chr20
# GAlignmentPairs object with 686908 pairs, strandMode=1, and 0 metadata columns:
#            seqnames strand   :            ranges  --            ranges
#               <Rle>  <Rle>   :         <IRanges>  --         <IRanges>
#        [1]    chr20      -   :       60100-60196  --       60007-60106
#        [2]    chr20      +   :       60246-60280  --       60248-60288
#        [3]    chr20      +   :       60255-60296  --       60257-60298
#        [4]    chr20      +   :       60260-60296  --       60262-60298
#        [5]    chr20      +   :       60261-60315  --       60263-60323
#        ...      ...    ... ...               ... ...               ...
#   [686904]    chr20      -   : 64332264-64332359  -- 64332180-64332268
#   [686905]    chr20      -   : 64332477-64332521  -- 64332469-64332513
#   [686906]    chr20      +   : 64332967-64333059  -- 64332969-64333061
#   [686907]    chr20      +   : 64333282-64333381  -- 64333312-64333411
#   [686908]    chr20      +   : 64334000-64334042  -- 64334002-64334044
#   -------
#   seqinfo: 25 sequences from an unspecified genome
atac_reads_chr20[1:2, ]

# GAlignmentPairs object with 4348926 pairs, strandMode=1, and 0 metadata columns:
#             seqnames strand   :              ranges  --              ranges
#                <Rle>  <Rle>   :           <IRanges>  --           <IRanges>
#         [1]     chr1      +   :          9996-10084  --         10277-10374
#         [2]     chr1      +   :          9997-10087  --         10082-10179
#         [3]     chr1      +   :         10003-10091  --         10123-10221
#         [4]     chr1      +   :         10004-10099  --         10009-10108
#         [5]     chr1      +   :         10004-10104  --         10037-10136
#         ...      ...    ... ...                 ... ...                 ...
#   [4348922]     chrX      -   : 156030710-156030739  -- 156030677-156030713
#   [4348923]     chrX      +   : 156030720-156030793  -- 156030723-156030795
#   [4348924]     chrX      +   : 156030726-156030759  -- 156030728-156030761
#   [4348925]     chrX      -   : 156030736-156030787  -- 156030725-156030785
#   [4348926]     chrX      -   : 156030748-156030783  -- 156030289-156030331
#   -------
#   seqinfo: 25 sequences from an unspecified genome

atac_reads_chr1X <- readGAlignmentPairs(sorted_bam, param = chr1X_param)
atac_reads_chr1X
# GAlignmentPairs object with 4348926 pairs, strandMode=1, and 0 metadata columns:
#             seqnames strand   :              ranges [read1]  --      ranges [read2]
#                <Rle>  <Rle>   :           <IRanges>  --           <IRanges>
#         [1]     chr1      +   :          9996-10084  --         10277-10374
#         [2]     chr1      +   :          9997-10087  --         10082-10179
#         [3]     chr1      +   :         10003-10091  --         10123-10221
#         [4]     chr1      +   :         10004-10099  --         10009-10108
#         [5]     chr1      +   :         10004-10104  --         10037-10136
#         ...      ...    ... ...                 ... ...                 ...
#   [4348922]     chrX      -   : 156030710-156030739  -- 156030677-156030713
#   [4348923]     chrX      +   : 156030720-156030793  -- 156030723-156030795
#   [4348924]     chrX      +   : 156030726-156030759  -- 156030728-156030761
#   [4348925]     chrX      -   : 156030736-156030787  -- 156030725-156030785
#   [4348926]     chrX      -   : 156030748-156030783  -- 156030289-156030331
#   -------
#   seqinfo: 25 sequences from an unspecified genome

# Gain information on the first or second/last alignment for each alignment pair in atac_reads_chr1X, containing the seqnames, strand, cigar, qwidth, start, end, width, qname, mapq, isize
read1 <- first(atac_reads_chr1X)
read2 <- second(atac_reads_chr1X)
read1[1, ] # This shows the first alignment in the first read
# GAlignments object with 1 alignment and 3 metadata columns:
#       seqnames strand       cigar    qwidth     start       end     width     njunc |              qname
#          <Rle>  <Rle> <character> <integer> <integer> <integer> <integer> <integer> |        <character>
#   [1]     chr1      +    30M2I59M        91      9996     10084        89         0 | SRR6870408.2987101
#            mapq     isize
#       <integer> <integer>
#   [1]         2       379
#   -------
#   seqinfo: 25 sequences from an unspecified genome
read2[1, ] # This shows the first alignment in the paired second read
# GAlignments object with 1 alignment and 3 metadata columns:
#       seqnames strand         cigar    qwidth     start       end     width     njunc |              qname
#          <Rle>  <Rle>   <character> <integer> <integer> <integer> <integer> <integer> |        <character>
#   [1]     chr1      - 15M1I62M1I21M       100     10277     10374        98         0 | SRR6870408.2987101
#            mapq     isize
#       <integer> <integer>
#   [1]         2      -379
#   -------
#   seqinfo: 25 sequences from an unspecified genome
```

```{r}
atac_reads_mainchr <- readGAlignmentPairs(sorted_bam, param = mainchr_param)
```


## Retrieve and summarize MapQ scores
```{r}
# Do our reads have the high-quality mapping?
read1_mapq <- mcols(read1)$mapq # mapq is a metadata column which can be accessed using mcols()
read2_mapq <- mcols(read2)$mapq
read1_mapq
read2_mapq

# *How our MapQ scores distribute across each read?
read1_mapq_dist <- table(read1_mapq); read2_mapq_dist <- table(read2_mapq)
read1_mapq_dist
read2_mapq_dist

# Visualize the distribution of MapQ scores per read
library(ggplot2)
plot_df <- data.frame(MapQ = c(names(read1_mapq_dist), names(read2_mapq_dist)),
                      Frequency = c(read1_mapq_dist, read2_mapq_dist),
                      Read = c(rep("Read1", length(read1_mapq_dist)), 
                              rep("Read2", length(read2_mapq_dist))))
plot_df$MapQ <- factor(plot_df$MapQ, 
                       levels = unique(sort(as.numeric(plot_df$MapQ))))
plot_df |>
  ggplot(mapping = aes(x = MapQ, y = Frequency, fill = MapQ)) + 
    geom_bar(stat = "identity") + 
    facet_grid(~Read) 
```


## Retrieve and summarize insert sizes 
```{r}
# Distance between read1 and read2 = distance between read2 and read1, so we only need to view the insert size of read1
insert_size <- mcols(read1)$isize |> abs() 
# insert_size_2 <- elementMetadata(read1)$isize |> abs()
insert_size
length(insert_size)
# [1] 4348926

# *How our insert sizes distribute?
(insert_size_dist <- table(insert_size))

# Visualize the distribution 
plot_df <- data.frame(Insert_size = as.numeric(names(insert_size_dist)),
                      Count = as.numeric(insert_size_dist)) # as.numeric only extracts the frequency values in a frequency table
(insert_size_plot <- plot_df |> 
  ggplot(mapping = aes(x = Insert_size, y = Count)) + 
  geom_line())

# Use the log2 transformation to the counts on y-axis to show the nucleosome patterning 
(insert_size_plot + 
    scale_y_continuous(trans = "log2"))

# The log2 transformed count of <100 bp is too high -> the concentration of transposases is too high
# Recap: ATAC-seq (Assay for Transposase-Accessible Chromatin with high-throughput sequencing) is a method for determining chromatin accessibility across the genome. It utilizes a hyperactive Tn5 transposase to insert sequencing adapters into open chromatin regions. High-throughput sequencing then yields reads that indicate these regions of increased accessibility.

# Nucleosome free (<100 bp), mononucleosome (180-247 bp), dinucleosome (315-437 bp)
(insert_size_plot + scale_y_continuous(trans = "log2") + 
  geom_vline(xintercept = c(180, 247),colour = "red") + 
  geom_vline(xintercept = c(315, 437),colour="darkblue") +
  geom_vline(xintercept = c(100), colour = "darkgreen"))
```


## Subset ATACseq reads by varying insert sizes
```{r}
# Nucleosome free (<100 bp), mononucleosome (180-247 bp), dinucleosome (315-437 bp)
# Why 100, not 147? Make sure the open region does not overlap the nucleosome location
atac_reads_chr1X_free <- atac_reads_chr1X[insert_size < 100, ]
atac_reads_chr1X_mono <- atac_reads_chr1X[insert_size > 180 & insert_size < 240, ]
atac_reads_chr1X_di <- atac_reads_chr1X[insert_size > 315 & insert_size < 437, ]
```

```{r}
read1 <- first(atac_reads_mainchr)
insert_size <- mcols(read1)$isize |> abs() 
atac_reads_mainchr_free <- atac_reads_mainchr[insert_size < 100, ]
```


## Save the subsetted ATACseq reads as BMA files
```{r}
sorted_bam
# [1] "./Sorted_ATAC_female_lung_bowtie2.bam"
free_bam <- gsub("\\.bam", "_chr1X_free\\.bam", sorted_bam) # "./Sorted_ATAC_female_lung_bowtie2_chr1X_free.bam"
mono_bam <- gsub("\\.bam", "_chr1X_mono\\.bam", sorted_bam)
di_bam <- gsub("\\.bam", "_chr1X_di\\.bam", sorted_bam)
full_bam <- gsub("\\.bam", "_chr1X\\.bam", sorted_bam)

library(rtracklayer)
export(object = atac_reads_chr1X_free, con = free_bam, format = "bam") # Save the object atac_reads_chr1X_free with the filename in `free_bam`
export(object = atac_reads_chr1X_mono, con = mono_bam, format = "bam")
export(object = atac_reads_chr1X_di, con = di_bam, format = "bam")
export(object = atac_reads_chr1X, con = full_bam, format = "bam")
```

```{r}
free_bam <- gsub("\\.bam", "_mainchr_free\\.bam", sorted_bam) 
export(object = atac_reads_mainchr_free, con = free_bam, format = "bam")
```


## Reconstruct one full-length fragment from 2 paired reads using `granges`
```{r}
# Show the GAlignmentPairs object with 1 pair
atac_reads_chr1X[1, ] 
# GAlignmentPairs object with 1 pair, strandMode=1, and 0 metadata columns:
#       seqnames strand :     ranges --      ranges
#          <Rle>  <Rle> :  <IRanges> --   <IRanges>
#   [1]     chr1      + : 9996-10084 -- 10277-10374

# Show the GRange object with 1 range, which collapses read 1 and read2 into one long read
atac_reads_chr1X_long <- granges(atac_reads_chr1X)
atac_reads_chr1X_long[1, ] 
# GRanges object with 1 range and 0 metadata columns:
#       seqnames     ranges strand
#          <Rle>  <IRanges>  <Rle>
#   [1]     chr1 9996-10374      +
#   -------
#   seqinfo: 25 sequences from an unspecified genome
```


## [QC] Estimate the non-redundant fraction (NRF) of the full-length fragments 
```{r}
duplicated(atac_reads_chr1X_long) |> length() # duplicated() returns a logical vector indicating which elements (rows) are duplicates.
# [1] 4348926
duplicated_frag <- sum(duplicated(atac_reads_chr1X_long))
duplicate_rate <- duplicated_frag/length(atac_reads_chr1X_long)
(non_redundant_fraction <- 1 - duplicate_rate) # [1] 0.7760803
```


## Construct an open region bigWig
```{r}
class(atac_reads_chr1X_long)
# [1] "GRanges"
# attr(,"package")
# [1] "GenomicRanges"
length(atac_reads_chr1X_long)
# [1] 4348926
open_region_bw <- gsub("\\.bam","_open_region\\.bw", sorted_bam)
coverage_bw <- coverage(atac_reads_chr1X_long,
                        weight = (10^6/length(atac_reads_chr1X_long))) # Normalize to total mapped reads
export.bw(coverage_bw, open_region_bw)
```