---
title: "Counting Reads with featureCounts"
author: "Taylor Jones (edited by Rutendo Sigauke)"
date: "`r Sys.Date()`"
output:
  pdf_document:
    toc: yes
  html_document:
    toc: yes
    toc_float:
      collapsed: yes
      smooth_scroll: yes
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
options(knitr.purl.inline = TRUE)
```


# Introduction

Here we will learn how to download a package, what metadata table is (and why it is important), and run featureCounts, which counts reads over genes.

# We will want to start fresh and clear our environment.

Start by clearing your console. To do this hit `Ctrl+l` or go to Edit --> Clear Console. Clear your environment and plots by hitting the broom icon in both those cells reset our working directory. 

## Alternatively, the counting can be submitted as an sbatch job
For this, you can open the `d6_featureCounts.R` and `d6_featureCounts.sbatch` files in vim to make edits.Then submit the job specifying resources needed. The `d6_featureCounts.R` will only contain R code shown in this notebook.

```
#!/bin/bash

#SBATCH --job-name=<NAME OF JOB>				 # Job name
#SBATCH --mail-type=ALL				 # Mail events (NONE, BEGIN, END, FAIL, ALL)
#SBATCH --mail-user=<YOUR E-MAIL ADDRESS>		 # Where to send mail
#SBATCH --nodes=1					 # Number of cores job will run on
#SBATCH --ntasks=4					 # Number of CPU (processers, tasks)
#SBATCH --time=1:00:00					 # Time limit hrs:min:sec
#SBATCH --partition short				 # Job queue
#SBATCH --mem=4gb					 # Memory limit
#SBATCH --output=/YOUR/EOFILES/PATH/%x_%j.out
#SBATCH --error=/YOUR/EOFILES/PATH/%x_%j.err

################### SET VARIABLES ######################################

FEATURECOUNTS=/PATH/TO/YOUR/d6_featureCounts.R

########################################################################
################## PRINT JOB INFO ######################################

printf "Sample ID: $ROOTNAME"
printf "\nDirectory: $PROJECT"
printf "\nRun on: $(hostname)"
printf "\nRun from: $(pwd)"
printf "\nScript: $0\n"
date

printf "\nYou've requested $SLURM_CPUS_ON_NODE core(s).\n"

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

Rscript $FEATURECOUNTS
```

# Set your working directory

```{r setwd, eval=FALSE}
workdir <- '/path/to/wd/'
setwd(workdir)
getwd()
```

# Install packages 

Like before: most packages you can install with this syntax: `install.packages("PACKAGE")`. Such as this package: `install.packages("tidyverse")` --great package!

HOWEVER, some of the bioinformatic software is located in Bioconductor, including `RSubread`(R wrapper for [**Subread**](http://subread.sourceforge.net/)) which contains `featureCounts`.


## Install `BiocManager`

Note: `eval=FALSE` means the code in the chunk is NOT evaluated. To evaluate the code chunk, change `eval=FALSE` to `eval=TRUE`.

The installations only needs to be run once. Turn back `eval=FALSE` after installation.

```{r installBiocmanager, eval=FALSE}
##This only needs to be run once. Turn back eval=FALSE after installation.

if (!requireNamespace("BiocManager", quietly = TRUE))  ## Install BiocManager
  install.packages("BiocManager")  #this will take a moment to run

```

## Install `RSubread`.

```{r installRsubread, eval=FALSE}
BiocManager::install("Rsubread")  #This actually installs RSubread
```

## Load the `Rsubread` package

This actually loads the library into the environment

```{r loadRsubread, eval=FALSE}

library("Rsubread")

```


Here, this calls for help of a package. This is useful if you have no idea where to start.

```{r rsubreadHelp, eval=FALSE}

??Rsubread # click on the Help pages search results

```


```{r featureCountsHelp, eval=FALSE}

??Rsubread::featureCounts # we can look specifically at featureCounts and the featureCount flags

```

# Counting reads with `Rsubread`

NOTE: We run `featureCounts` on a cluster because bam files are LARGE.
We want to run this on the cluster but using command line `R` is not always great.
An alternative is to generate a script like this then submit this script to the cluster in an sbatch script.

## Inputs for `featureCounts`

### Bam files

Read in bam file list
```{r bamdir, eval=FALSE}
bamdir <- '/path/to/bam/files/'

filelist <- c(paste0(bamdir, "chr21Eric_repA.RNA.sorted.bam"),
           paste0(bamdir, "chr21Ethan_repA.RNA.sorted.bam"))

outdir <- paste(workdir,'/', 'counts', '/', sep='') ##naming our outdir
dir.create(outdir) ###creating the directory


```

### Annotation file in GTF format 

We use a GTF -- this is necessary for `featureCounts` and more useful for RNA-seq data as it is cleaner/faster for counting over exons

```{r loadGTF, eval=FALSE}

hg38gtf <- "/path/to/hg38_ucsc_genes_chr21.gtf"

```


## Read counting using `featureCounts`

These parameters below are for specifying how to count over the gtf file.

```{r featurecounts, eval=FALSE}
# Sink will save out stdout -- check all of these settings (e.g. isPairedEnd)!

# Define root name for output files
fileroot <- "chr21_Ethan_Eric"

sink(paste0(outdir, fileroot, "_featureCounts_gene_rnaseq.txt")) #we want to save the output as a file
fc <- featureCounts(files=filelist,
                    annot.ext=hg38gtf,
                    isGTFAnnotationFile=TRUE,
                    GTF.featureType="exon",
                    GTF.attrType="gene_id",
                    useMetaFeatures=TRUE,
                    allowMultiOverlap=TRUE,
                    largestOverlap=TRUE,
                    countMultiMappingReads=TRUE,
                    isPairedEnd=TRUE,
                    strandSpecific=0,
                    nthreads=4)  #when you move to a bigger machine change to 8

#Flag info here: ??Rsubread::featureCounts for specifics on each flag

sink()
```

## Write results 

We'll keep our accession number (GeneID) and length for filtering genes in DESeq2. We will also be writing different files to separate files (listed below).

- `featureCounts_gene_rnaseq.txt` : GeneID, Length, Counts
- `coverage.csv` : Counts
- `.stat.csv` : Coverage Statistics
- `.annotation.csv` : GeneID, Chromosome, Start, End, Strand, Length
- `.targets.csv` : file names for input `bam` files

```{r writecounts, eval=FALSE}
# Write outout file with write.table
write.table(x=data.frame(fc$annotation[,c("GeneID","Length")],
                         fc$counts,stringsAsFactors=FALSE),
            paste0(outdir, fileroot, file="_featureCounts_gene_rnaseq.txt"),
            quote=FALSE,sep="\t",
            row.names=FALSE)

# Write out separate files for counts and other count statistics
fileroot <- "chr21_Ethan_Eric_repA"
write.csv(fc$counts, 
          paste(outdir, fileroot,".coverage.csv", sep=""))
write.csv(fc$stat,
          paste(outdir, fileroot,".stat.csv", sep=""))
write.csv(fc$annotation, 
          paste(outdir, fileroot,".annotation.csv", sep=""))
write.csv(fc$targets,
          paste(outdir, 
                fileroot,".targets.csv", sep=""))

```

