Introduction

This workshop will be an overview of an RNASeq workflow using widely accepted analysis techniques. We will perform a case study examining gene knockout of Atg5 in macrophages of B6 mice.

This Workflow will be posted to my github repository with the outputs from the commands

WorkFlow Overview

  1. Quality control (Skipping this step)
  2. Map reads to genome
  3. Count reads mapped to annotated genes
  4. Data Assessment/Quality Control
  5. Detect differentially expressed genes
  6. Data Visualization (Pretty Plots)
  7. Pathway analysis

Data originate from a study on norovirus infection (https://www.ncbi.nlm.nih.gov/sra/?term=PRJEB10074)

Case Study

Black 6 Mice +/- Atg5 (5 knockouts and 9 controls)   - Macrophage cells
Sequencing – Ilumina HiSeq
  - Single-end
  - 50 bp reads

READ MAPPING

STEP 1: READ MAPPING with STAR

OPTIONS (more options available in the STAR manual)

--runThreadN NumberOfThreads
--genomeDir /path/to/genomeDir
--readFilesIn /path/to/read1 [/path/to/read2] 
--readFilesCommand (uncompression command)*
--outFileNamePrefix 
--outFilterMismatchNmax N (recommended 0.06*readLength)
--outSamtype (sorted or unsorted)

STAR Command

$ STAR --runThreadN 2 --outBAMsortingThreadN 2 \
  --genomeDir <GENOME_INDEX_DIRECTORY> \
  --readFilesIn <READS_DIR><file> \
  --outFileNamePrefix <OUTPUT_DIR><prefix>. \
  --outFilterMismatchNmax 3
  --outReadsUnmapped Fastx \
  --outSAMtype BAM SortedByCoordinate

Now run the STAR command with “AACGCATT.fq” as the input file and output it to the STAR_Mapping directory. ####Note: For each exercise, “” will be replaced with “AACGCATT”

Once the command has completed Take a look at the files inside STAR_Mapping generated by STAR to make sure all files have been generated. You should have files with the same extensions as in the image below. alt text

We can use the following command to print the contents of the Log file to your screen to get some general mapping information

$ cat AACGCATT.Log.final.out

You will notice that the number of input reads is listed in this file. Was your answer about the number of reads correct?

PAUSE 🛑

STEP 1.2: VIEW ALIGNMENTS IN IGV

  1. Load up x2g0 and select the terminal icon
  2. enter the command “igv.sh” and a browser should load
  3. Select File tab -> Load from File (~/workshop_materials/Transcriptomics/BAMS/AACGCATT.bam)
  4. Select Genome tab -> Load Genome From File (~/workshop_materials/Transcriptomics/MouseGenome_GRCm38/mm_ref_GRCm38.p2_chr1.fa)
  5. In the Search box type in a chromosome location
    i.) We only have chromosome one in this case ii.) for an example, type in “chr1:9,841,441-9,841,599” iii.) double click to zoom in and let’s look at the read mapping together

return to Transcriptomics.

For convenience, lets copy the bam file generated by STAR into the BAMS directory.

STEP 2: Counting

Before we execute the counting command, we need to convert the bam files into sam format.
Move into the BAMS directory if you aren’t there already.
Run the below command to convert the bam file into sam format

$  samtools view -h <prefix>.bam > <prefix>.sam

Now that we have the alignment file in sam output, use one of your unix commands to take a looks at some of the lines in this file.

PAUSE 🛑

STEP 2: COUNTING With HTSeq

PAUSE 🛑

Considerations for real studies:

When performing an RNASeq study, we have multiple samples and multiple output files. We don’t want to perform these commands one by one, so we can use some unix tools to automate these repetitive tasks. This is beyond the scope of this workshop, but I recommend you take a look at the links below at some point, so you can learn how to write some basic scripts execute a command like STAR iteratively over multiple files

   for-do-loop - iterate over multiple files
   using variables
   grep
   “|” - redirect output of a command to another command

Below is an example of how we use for-do-loops to execute the HTSeq command across multiple files. Try not to let yourself be overwhelmed with the syntax. Slowly, over time, you will become equipped with all the tools to write scripts similary to automate tasks.

alt text

alt text

In addition, there are many variations of the above script that could accomplish the same output. As you become more familiar with unix and coding in the terminal, you will acquire skills to make your analyses much more efficient.

Up until now, we have been working on a single file and a subset of it’s reads. I have placed the HTSeq output for each file in a github repository. You can go to the Counts directory and type the following command to retrieve the files.

$ git pull

Open a browser and log into RStudio

#load libraries
library("phyloseq"); packageVersion("phyloseq")
library("DESeq2"); packageVersion("DESeq2")
library("ggplot2"); packageVersion("ggplot2")
library("ggrepel"); packageVersion("ggrepel")
library("data.table"); packageVersion("data.table")

Read HTSeq output into DESeq2 ###The following commands are adopted from the DESeq2 manual

First we need to read in the files

directory <- "/home/genomics/workshop_materials/Transcriptomics/Counts/"
sampleFiles <- list.files(path=directory,pattern=".txt")
sampleFiles

Now read in your metadata from the “Conditions.txt” file

sampleCondition<- read.table("/home/genomics/workshop_materials/Transcriptomics/Conditions.txt",head=TRUE) #file with sample data to be compared
sampleCondition

sampleCondition <- sampleCondition[order(sampleCondition$SampleID),] #We need to order this table by SampleID before merging with the sampleFiles
sampleCondition

combine sample data into a data frame with file names

sampleTable <- data.frame(sampleName = sampleFiles, fileName = sampleFiles, condition = sampleCondition)
sampleTable #verify that the "fileName" column matches the condition.SampleID column

Now we have our metadata formatted in a way that can be used by DESeq2, but we need to create a data matrix that can be used by the DESeq function

ddsCounts <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design = ~ condition.condition) #creates DESEqDataSet for DESeq function

PAUSE 🛑

Before we run the DESeq function, we need to do a final quality control step to make sure there are no batch effects or outliers

#The data first need to be transformed using the variance stabilizing transformation
vsd <- vst(ddsCounts, blind=FALSE) #variance stabilizing transformation
pcaData <- plotPCA(vsd, intgroup=c("condition.condition"),returnData=TRUE)
percentVar <- round(100 * attr(pcaData,"percentVar"))
head(pcaData)

ggplot(pcaData, aes(PC1, PC2, color=condition.condition))+
  geom_text(mapping = aes(label = name),size=2.5)+
  geom_point(size=2.5) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  xlim(-25,25)+
  coord_fixed()

We can see that the sample “CGCAATTC.txt” performs drastically differently than the rest of the samples, so we need to remove it

vsd2 <- vsd[,vsd$condition.SampleID != "CGCAATTC.txt"]

pcaData2 <- plotPCA(vsd2, intgroup=c("condition.condition"),returnData=TRUE)
percentVar <- round(100 * attr(pcaData2,"percentVar"))
head(pcaData2)

ggplot(pcaData2, aes(PC1, PC2, color=condition.condition))+
  #geom_text(mapping = aes(label = name),size=1.5)+
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed()

You can see, there is now separation between the two conditions, so we won’t do anymore sample removal based on the PCA results

Now we have to remove the outlier(s) from the dds counts matrix since this is what the DESeq function will be run on

ddsCounts2 <- ddsCounts[,ddsCounts$condition.SampleID != "CGCAATTC.txt"]

#Now we can run the differential expression analysis
dds <- DESeq(ddsCounts2) #DESeq2 command for differential expression test

Extracting DESeq2 Results

Below is a description of the column headers that will appear in the DESeq2 results table: * Basemean - mean of normalized counts for all samples * log2foldchange * lfcse - standard error * stat - wald statistic (log2foldchange / lfcse) * pvalue * padj - pvalue adjusted for false discovery ####Warning: You MUST resist the temptation to use the p-value when evaluating DESeq2 results. The adjusted p-value (padj) accounts for multiple hypothesis testing

#text file with ensembl ID's and corresponding gene names
mmGenes = read.table("/home/genomics/workshop_materials/Transcriptomics/MouseGeneTable.txt", header=TRUE, row.names = 1)
mmGenes = as.matrix(mmGenes) #I convert the data frame to a matrix to get a downstream function, cbind, to work properly
head(mmGenes)

extract results table from DESeq2

res <- results(dds)
res.dds <- res[order(res$padj),] #order table by adj p value
summary(res.dds) #lets take a look at how many genes are up-regulated/down-regulated
  
  # Quick check of factor levels
mcols(res.dds, use.names = TRUE) #this is important for results interpretation

We can add the gene names to the results by merging the mmGenes table with the res.dds table

res.dds = cbind(as(res.dds, "data.frame"), Gene_Name = as(mmGenes[rownames(res.dds), ],"matrix"))
head(res.dds)

Make a new object to add significance boolean for plotting

resdt.dds = data.table(as(res.dds, "data.frame"),
                 keep.rownames = TRUE)
setnames(resdt.dds, "rn", "Ensembl_ID") 

resdt.dds
resdt.dds[, Significant := padj < .1]
resdt.dds[!is.na(Significant)]
resdt.dds
  
write.table(res.dds,file="/home/genomics/workshop_materials/Transcriptomics/FF_vs_Cre.DESeq.out.txt",quote=F,sep="\t")

make rnk table for Gene Set Enrichment Analysis

RNK = data.table(Gene_Name = res.dds$Gene_Name, stat = res.dds$stat) #These are 2 columns from our deseq2 output that we need for GSEA
RNK = subset(RNK, stat != "NA")
head(RNK)
write.table(RNK, "/home/genomics/workshop_materials/Transcriptomics/FF_vs_Cre.DESeq.rnk",quote=F,sep="\t", row.names = F)

Volcano plots are a good way of visualizing the effects of the experimental conditions on the cells

volcano = ggplot(
    data = resdt.dds[!is.na(Significant)],
    mapping = aes(x = log2FoldChange,
                  y = -log10(padj),
                  color = Significant,
                  label = Ensembl_ID, label1 = Gene_Name)) +
    theme_bw() +
    geom_point() + 
    geom_point(data = resdt.dds[(Significant)], size = 7, alpha = 0.7) + #Larger circles for the significant values
    theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5)) +
    geom_hline(yintercept = -log10(.1)) +
    ggtitle("DESeq2 Negative Binomial Test Volcano Plot \nProximal Colon Abx vs no Abx") +
    theme(axis.title = element_text(size=12)) +
    theme(axis.text = element_text(size=12)) +
    theme(legend.text = element_text(size=12)) +
    geom_vline(xintercept = 0, lty = 2)
#volcano
volcano + xlim(-7,7) #xlim isn't necessary, but I add it here to make the plot symmetrical about the y axix
summary(resdt.dds)

save.image("./Transcriptomics.RData") #it's to to save the image so all the object that have been initiated can be accessed at any time

After running GSEA, we may want to make some heat maps. To do that we need to merge HTSeq counts files with Gene list file to make a counts table

counts = mmGenes #initialize counts table as gene list with gene names and ensemble IDs
#subset sampleFiles to remove the outlier we removed earlier
sampleFiles

#iterate over each HTSeq file
for(i in 1:length(sampleFiles)){ 
  file <- read.table(paste(directory,sampleFiles[i],sep=""),header=FALSE,sep="\t",row.names = 1)

  colnames(file) = c(sub('\\.txt$', '', sampleFiles[i])) #add the file name minus the ".txt" suffix to the header so we can identify each column
  
  dat = subset(file, !(row.names(file)%in%c("__no_feature","__ambiguous","__too_low_aQual","__not_aligned","__alignment_not_unique"))) #remove the extra fields from the counts table
  
 counts = cbind(as(counts,"matrix"),as(dat,"data.frame")) #each time the loop iterates, it will add a new column for the file's counts
  
}
head(counts)
write.table(counts,"/home/genomics/workshop_materials/Transcriptomics/CountsTable.txt",quote=F,sep="\t")

Because R…We have to open the file “CountsTable.txt” and add the column header “Ensembl_ID” followed by a tab to the file, so the columns match with their headers. You can open the file in RStudio by double clicking the file

Gene Set Enrichment Analysis (GSEA)

Open x2goclient Because R…We have to open the file “CountsTable.txt” and add the column header “Ensembl_ID” followed by a tab to the file, so the columns match with their headers. Click the GSEA tab (a java browser should open)

Enrichment Plots

1. Open GSEA

2. select the “Load data” tab and select the rnk file (blue square by the name)

alt text

alt text

3. Tools/GSEAPreranked

3.1 Select Gene Sets

   i) gseaftp.broadinstitute.org://pub/gsea/gene_sets/h.all.v6.2.symbols.gmt (Hallmark Gene Symbols)
   ii) gseaftp.broadinstitute.org://pub/gsea/gene_sets/c2.cp.v6.2.symbols.gmt (Canonical Pathways)