03.00-D-Apul-RNAseq-gene-expression-DESeq2

1 Set Bash variables

{
echo "#### Assign Variables ####"
echo ""

echo "# Data directories"
echo 'export timeseries_dir=/home/shared/8TB_HDD_02/shedurkin/timeseries_molecular'
echo 'export output_dir_top="${timeseries_dir}/D-Apul/output/03.00-D-Apul-RNAseq-gene-expression-DESeq2"'
echo ""



echo "# Output files"
echo 'export coldata="${output_dir_top}/DESeq2-coldata.tab"'
echo ""

echo "# Paths to programs"
echo 'export programs_dir="/home/shared"'


echo "# Set number of CPUs to use"
echo 'export threads=40'
echo ""


echo "# Print formatting"
echo 'export line="--------------------------------------------------------"'
echo ""
} > .bashvars

cat .bashvars

#### Assign Variables ####

# Data directories
export timeseries_dir=/home/shared/8TB_HDD_02/shedurkin/timeseries_molecular
export output_dir_top="${timeseries_dir}/D-Apul/output/03.00-D-Apul-RNAseq-gene-expression-DESeq2"

# Output files
export coldata="${output_dir_top}/DESeq2-coldata.tab"

# Paths to programs
export programs_dir="/home/shared"
# Set number of CPUs to use
export threads=40

# Print formatting
export line="--------------------------------------------------------"

2 Set R variables

# Define the output directory path
output_dir <- "../output/03.00-D-Apul-RNAseq-gene-expression-DESeq2/"

# Set desired false discovery rate threshold (i.e. adjusted p-value, padj)
fdr <- 0.05

# Set log2 fold change threshold (a value of '1' is equal to a fold change of '2')
log2fc <- 1

3 Create coldat file

# Load bash variables into memory
source .bashvars

# Create output directory, if it doesn't exist
mkdir --parents "${output_dir_top}"

# Create associative array with sample and timepoint
metadata="../../M-multi-species/data/rna_metadata.csv"

# Create DESeq2-formatted coldata file

## Create header
printf "\t%s\t%s\n" "time.point" "colony.id"> "${coldata}"

## Read the metadata file line by line
while IFS=',' read -r sample_number sample_name plate well_number azenta_sample_name colony_id timepoint sample_type species_strain SampleBuffer; do
    # Check if the species is "Acropora pulchra"
    if [[ "${species_strain}" == "Acropora pulchra" ]]; then
      printf "%s\t%s\t%s\n" "${azenta_sample_name}" "${timepoint}" "${colony_id}"
    fi
done < <(tail -n +2 "${metadata}") \
| sort -k1,1 \
>> "${coldata}"



## Tab-delimited output of sample and timepoint
for sample in "${!sample_timepoint_map[@]}"
do
  printf "%s\t%s\n" "$sample" "${sample_timepoint_map[$sample]}"
done | sort -k1,1 \
>> "${coldata}"

# Peek at output
head "${coldata}" | column -t

echo ""
echo "${line}"
echo ""

wc -l "${coldata}"

echo ""
echo "${line}"
echo ""

echo "Colony counts:"
echo ""
awk 'NR > 1 {print $3}' "${coldata}" | sort | uniq -c

4 Read in gene counts and coldata files

gene.counts <- as.matrix(read.csv(file = "../output/02.20-D-Apul-RNAseq-alignment-HiSat2/apul-gene_count_matrix.csv", row.names="gene_id", check.names = FALSE))

coldata <- read.csv(file = "../output/03.00-D-Apul-RNAseq-gene-expression-DESeq2/DESeq2-coldata.tab", row.names=1, sep = "\t")
coldata$time.point <- factor(coldata$time.point)

head(gene.counts)

head(coldata)

4.1 Verify rownames match

all(rownames(coldata) == colnames(gene.counts))

5 Create DESeq2 data set

dds <- DESeqDataSetFromMatrix(countData = gene.counts,
                              colData = coldata,
                              design = ~ time.point + colony.id)
dds

5.1 Add gene columun feature

featureData <- data.frame(gene=rownames(gene.counts))
mcols(dds) <- DataFrame(mcols(dds), featureData)
mcols(dds)

dds$time.point <- factor(dds$time.point, levels = c("TP1","TP2", "TP3", "TP4"))

6 DESeq analysis

dds <- DESeq(dds)

6.1 Pairwise results tables

6.1.1 Full

tp1.v.tp2.results <- results(dds, contrast=c("time.point","TP1","TP2"), alpha = fdr, lfcThreshold = log2fc)
tp1.v.tp3.results <- results(dds, contrast=c("time.point","TP1","TP3"), alpha = fdr, lfcThreshold = log2fc)
tp1.v.tp4.results <- results(dds, contrast=c("time.point","TP1","TP4"), alpha = fdr, lfcThreshold = log2fc)
tp2.v.tp3.results <- results(dds, contrast=c("time.point","TP2","TP3"), alpha = fdr, lfcThreshold = log2fc)
tp2.v.tp4.results <- results(dds, contrast=c("time.point","TP2","TP4"), alpha = fdr, lfcThreshold = log2fc)
tp3.v.tp4.results <- results(dds, contrast=c("time.point","TP3","TP4"), alpha = fdr, lfcThreshold = log2fc)

tp2.v.tp4.results

summary(tp2.v.tp4.results)

table(tp2.v.tp4.results$padj < 0.05)

7 Write DDS results tables to CSVs

# Create a named list of the data frames
results_list <- list(
  tp1.v.tp2.results = tp1.v.tp2.results,
  tp1.v.tp3.results = tp1.v.tp3.results,
  tp1.v.tp4.results = tp1.v.tp4.results,
  tp2.v.tp3.results = tp2.v.tp3.results,
  tp2.v.tp4.results = tp2.v.tp4.results,
  tp3.v.tp4.results = tp3.v.tp4.results
)

# Loop through the list and write each data frame to a CSV file in the specified directory
for (df_name in names(results_list)) {
  write.csv(results_list[[df_name]], file = paste0(output_dir, df_name, ".table.csv"), row.names = FALSE, quote = FALSE)
}

8 Variance stabilizing transformations (VST)

• Here we transform counts using a variance stabilizing transformation (VST), since we have many samples.

vsd <- varianceStabilizingTransformation(dds, blind=FALSE)

NOTE: Hover over points to see the sample numbers

9 Plotting

9.1 Sample distances

sampleDists <- dist(t(assay(vsd)))

sampleDistMatrix <- as.matrix( sampleDists )

rownames(sampleDistMatrix) <- paste( vsd$colony.id, vsd$time.point, sep = " - " )

colnames(sampleDistMatrix) <- NULL

colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)

pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colors)

9.2 PCA - All time points

Visualize sample clustering via PCA (after transformation)

# PCA with points color coded by time point 
plotPCA(vsd, intgroup = c("time.point"))

# PCA with points color coded by colony ID 
plotPCA(vsd, intgroup = c("colony.id"))

Time points 1 and 4 are clustering together, and time points 2 and 3 are clustering together. It also looks like colonies cluster somewhat.

9.3 Heatmap - All time points

top_20_counts_all <- order(rowMeans(counts(dds,normalized=TRUE)),
                decreasing=TRUE)[1:200]

timepoint_annotation = colData(dds) %>% as.data.frame() %>% select(time.point)


pheatmap(assay(vsd)[top_20_counts_all,], 
         cluster_rows=FALSE, 
         show_rownames=FALSE,
         cluster_cols=TRUE, 
         annotation_col = timepoint_annotation)