Mon. Sep. 05, 2022

Project name: Berdahl-sockeye-salmon
Funding source: unknown
Species: Oncorhynchus nerka
variable: behavior: territorial, social

« previous notebook entry « | » next notebook entry »

Background

1. 30 sockeye salmon sampled; 1-15: territorial, 16-30: social
2. brain, liver, and gonad saved in RNAlater - frozen at -80C
3. RNA extracted - see github issues: 1 and 2
4. RNA submitted to UT Austin GSAF on 5/24 for Tag-seq. Received gonad samples on 7/25. See github issue.
5. Gave go ahead to GSAF to complete rest of tagseq on 8/12.
6. Completed tagseq processing in last notebook entry.

Tag-seq analysis - Gonad Samples

Raw sequences were processed using the this R script. The resulting Gene count matrix was the used to run DEG analysis using DESEQ2, as outlined in this R script.

Using this treatment conditions datasheet, and this genome feature table, I ran DESeq2 using this gene count matrix.

# Filter data
coldata_trt <- coldata # %>% filter(group == "A" | group == "B")
cts_trt     <- subset(cts, select=row.names(coldata_trt))
# Calculate DESeq object
dds <- DESeqDataSetFromMatrix(countData = cts_trt,
                              colData = coldata_trt,
                              design = ~ trt)
dds <- DESeq(dds)
resultsNames(dds) # lists the coefficients

I filtered out genes with cumulative counts less than 10 across all samples:

# https://support.bioconductor.org/p/110307/
dds <- dds[rowSums(DESeq2::counts(dds)) >= 10,]

I then looked at the pheatmaps and multiqc report, I removed two low quality samples (one from each treatment, major outliers in PCA and ~88% overlap with rest of samples):

# Remove bad samples: C05, C17
coldata <- coldata[!(row.names(coldata) %in% c('C05', 'C17')),]
cts <- as.matrix(subset(cts, select=-c(C05, C17)))
coldata %>% dplyr::count(trt)
all(colnames(cts) %in% rownames(coldata))

Using the good samples, I compared territorial vs. social salmon and generated the following DEG list. From this list I made the following volcano plot, PCA (+ pairs plot), and pheat map comparing the impact of treatment: </br>

I also tested the impact of different shrinkage estimators (normal, apeglm, or ashr). The makers of DESeq2 suggest that the apeglm is the best.

After running all estimators on the DEG list and filtering by a log2fold change cutoff of 1.5 and a p value cutoff of 0.05 I got the following results:

"genes_before_filtering" 37942
"genes_after_filtering" 33153
"genes_dropped" 4789
"DEGs_all-genes" 33153
"DEGs_all-genes-normal" 33153
"DEGs_all-genes-apeglm" 33153
"DEGs_all-genes-ashr" 33153
"DEG_unshrunken-p0.05_lfc1.5" 146
"DEG_normal-lfc1.5" 13
"DEG_apeglm-s0.005_lfc1.5" 106
"DEG_ashr-p0.05_lfc1.5" 12

Using the apeglm shrinkage estimator and significance cutoffs, I generated the following volcano plot:

Here is the full significant apeglm-DEG list.

I threw this list into DAVID and used the O.nerka genome as a background. This resulted in 53 GO terms, 13 of which were associated with KEGG pathways.

GO terms:

KEGG Pathways: