GSE64570
This report outlines the analysis of a subset of the GSE64570 data set, consisting of 3 replicates of wildtype zebrafish injected with water, and 3 replicates of tlr5a morphants injected with Flagellin.
Preparation
Reference directories and packages
basedir <- "/home/charlotte/gene_vs_tx_quantification"
refdir <- "/home/charlotte/gene_vs_tx_quantification/annotation"
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(tximport))
suppressPackageStartupMessages(library(iCOBRA))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(BiocParallel))
suppressPackageStartupMessages(library(DESeq2,
lib.loc = "/home/charlotte/R/x86_64-pc-linux-gnu-library/3.2"))
suppressPackageStartupMessages(library(DEXSeq))Metadata definition
A metadata table was generated using the information provided in Gene Expression Omnibus.
meta <- read.delim(paste0(basedir, "/data/gse64570/gse64570_phenodata.txt"),
header = TRUE, as.is = TRUE)
rownames(meta) <- meta$srr.id
meta## gsm.id srr.id condition
## SRR1736697 GSM1574496 SRR1736697 wildtype_water
## SRR1736698 GSM1574497 SRR1736698 wildtype_water
## SRR1736699 GSM1574498 SRR1736699 wildtype_water
## SRR1736715 GSM1574514 SRR1736715 tlr5a_flagellin
## SRR1736716 GSM1574515 SRR1736716 tlr5a_flagellin
## SRR1736717 GSM1574516 SRR1736717 tlr5a_flagellinFASTQ file download
The code below downloads the FASTQ files for the six samples from the SRA.
fastq_files <- c(paste0("ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR173/00",
c(7:9, 5:7), "/",
meta$srr.id, "/", meta$srr.id, ".fastq.gz"))
for (fq in fastq_files) {
if (!file.exists(paste0(basedir, "/data/gse64570/fastq/", basename(fq)))) {
cmd <- paste0("wget -P ", basedir, "/data/gse64570/fastq ", fq)
message(cmd)
system(cmd)
} else {
message(paste0(fq, " is already downloaded."))
}
}## ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR173/007/SRR1736697/SRR1736697.fastq.gz is already downloaded.## ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR173/008/SRR1736698/SRR1736698.fastq.gz is already downloaded.## ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR173/009/SRR1736699/SRR1736699.fastq.gz is already downloaded.## ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR173/005/SRR1736715/SRR1736715.fastq.gz is already downloaded.## ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR173/006/SRR1736716/SRR1736716.fastq.gz is already downloaded.## ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR173/007/SRR1736717/SRR1736717.fastq.gz is already downloaded.Reference file preparation and index building
We use the GRCz10 Danio rerio Ensembl annotation for the analysis. Here, we download the fasta file with the cDNA sequences and build an index for quantification with Salmon.
cdna_fasta <- paste0(refdir, "/Danio_Rerio_GRCz10/Danio_rerio.GRCz10.cdna.all.fa.gz")
salmon_index <- paste0(refdir, "/Danio_Rerio_GRCz10/salmon_index/Danio_Rerio_GRCz10.sidx")
############################### cDNA FASTA ################################
if (!file.exists(cdna_fasta)) {
cmd <- paste0("wget -P ", refdir, "/Danio_Rerio_GRCz10 ",
"ftp://ftp.ensembl.org/pub/release-82/fasta/danio_rerio",
"/cdna/Danio_rerio.GRCz10.cdna.all.fa.gz")
message(cmd)
system(cmd)
cmd <- paste0("gunzip -c ", cdna_fasta, "> ", gsub("\\.gz$", "", cdna_fasta))
message(cmd)
system(cmd)
} else {
message(paste0("Reference cDNA fasta is already downloaded."))
}## Reference cDNA fasta is already downloaded.## Build Salmon index
cmd <- paste("salmon index -i", salmon_index, "-t", gsub("\\.gz$", "", cdna_fasta), "-p 5 --type quasi")
message(cmd)## salmon index -i /home/charlotte/gene_vs_tx_quantification/annotation/Danio_Rerio_GRCz10/salmon_index/Danio_Rerio_GRCz10.sidx -t /home/charlotte/gene_vs_tx_quantification/annotation/Danio_Rerio_GRCz10/Danio_rerio.GRCz10.cdna.all.fa -p 5 --type quasiif (!file.exists(paste0(salmon_index, "/hash.bin"))) {
system(cmd)
} else {
message("Salmon index already exists.")
}## Salmon index already exists.Definition of gene-to-transcript mapping
Next, we derive a mapping between transcript and gene identifiers from the cDNA file downloaded above.
feature_lengths_file <- paste0(refdir, "/Danio_Rerio_GRCz10/feature_lengths.Rdata")
tx_gene_file <- paste0(refdir, "/Danio_Rerio_GRCz10/tx_gene_map.Rdata")## Calculate gene and transcript lengths, get gene-transcript mapping
if (!file.exists(feature_lengths_file)) {
calc_lengths_mapping(gtf = NULL, cdna_fasta = cdna_fasta,
feature_lengths_file = feature_lengths_file,
tx_gene_file = tx_gene_file)
} else {
message("feature lengths and tx-to-gene map already calculated.")
}## feature lengths and tx-to-gene map already calculated.load(tx_gene_file)Salmon abundance quantification
Next, we use Salmon to calculate transcript abundance estimates for each of the six samples.
salmon_basedir <- paste0(basedir, "/quantifications/gse64570/salmon")
fqs <- list.files(paste0(basedir, "/data/gse64570/fastq"), full.names = TRUE)
names(fqs) <- gsub("\\.fastq.gz", "", basename(fqs))
for (i in 1:length(fqs)) {
if (!file.exists(paste0(salmon_basedir, "/", names(fqs)[i], "/quant.sf"))) {
cmd <- sprintf("bash -c 'salmon quant -i %s -l U -r %s -p 5 -o %s'",
salmon_index,
paste0("<(gunzip -c ", fqs[i], ")"),
paste0(salmon_basedir, "/", names(fqs)[i]))
cat(cmd, "\n")
system(cmd)
} else {
cat("Salmon results for", names(fqs)[i], "already exist.\n")
}
}## Salmon results for SRR1736697 already exist.
## Salmon results for SRR1736698 already exist.
## Salmon results for SRR1736699 already exist.
## Salmon results for SRR1736715 already exist.
## Salmon results for SRR1736716 already exist.
## Salmon results for SRR1736717 already exist.Summarization of Salmon results and offset estimation
We use the tximport package (https://github.com/mikelove/tximport) to generate count matrices and offset matrices (average transcript lengths) from the Salmon transcript-level estimates. We generate two different count matrices (simplesum and scaledTPM), and additionally create offsets to be used with the simplesum matrix.
salmon_files <- list.files(salmon_basedir, pattern = "SRR", full.names = TRUE)
salmon_files <- salmon_files[file.info(salmon_files)$isdir]
salmon_files <- paste0(salmon_files, "/quant.sf")
salmon_files <- salmon_files[file.exists(salmon_files)]
names(salmon_files) <- basename(gsub("/quant.sf", "", salmon_files))
txi_salmonsimplesum <- tximport(files = salmon_files, type = "salmon", txIn = TRUE,
txOut = FALSE, countsFromAbundance = "no",
gene2tx = gene2tx)## reading in files## 1## 2## 3## 4## 5## 6## ## summarizing abundance## summarizing counts## summarizing lengthtxi_salmonscaledtpm <- tximport(files = salmon_files, type = "salmon", txIn = TRUE,
txOut = FALSE, countsFromAbundance = "scaledTPM",
gene2tx = gene2tx)## reading in files## 1## 2## 3## 4## 5## 6## ## summarizing abundance## summarizing counts## summarizing lengthtxi_salmontx <- tximport(files = salmon_files, type = "salmon", txIn = TRUE,
txOut = TRUE, countsFromAbundance = "no", gene2tx = gene2tx)## reading in files## 1## 2## 3## 4## 5## 6## salmon_quant <- list(geneCOUNT_sal_simplesum = txi_salmonsimplesum$counts,
geneCOUNT_sal_scaledTPM = txi_salmonscaledtpm$counts,
avetxlength = txi_salmonsimplesum$length,
geneTPM_sal = txi_salmonsimplesum$abundance,
txTPM_sal = txi_salmontx$abundance,
txCOUNT_sal = txi_salmontx$counts,
txi_salmonsimplesum = txi_salmonsimplesum,
txi_salmonscaledtpm = txi_salmonscaledtpm,
txi_salmontx = txi_salmontx)Differential expression analysis
Given the gene count matrices defined above we apply edgeR and DESeq2 to perform differential gene expression. For the simplesum matrix, we also apply edgeR and DESeq2 using the offsets derived from the average transcript lengths (simplesum_avetxl).
edgeR
res_sal_simplesum_edgeR <- diff_expression_edgeR(counts = salmon_quant$geneCOUNT_sal_simplesum,
meta = meta, cond_name = "condition",
sample_name = "srr.id",
gene_length_matrix = NULL)
res_sal_simplesum_avetxl_edgeR <- diff_expression_edgeR(counts = salmon_quant$geneCOUNT_sal_simplesum,
meta = meta, cond_name = "condition",
sample_name = "srr.id",
gene_length_matrix = salmon_quant$avetxlength)
res_sal_scaledTPM_edgeR <- diff_expression_edgeR(counts = salmon_quant$geneCOUNT_sal_scaledTPM,
meta = meta, cond_name = "condition",
sample_name = "srr.id",
gene_length_matrix = NULL)Diagnostics
dfh <- data.frame(pvalue = c(res_sal_scaledTPM_edgeR$tt$PValue,
res_sal_simplesum_avetxl_edgeR$tt$PValue,
res_sal_simplesum_edgeR$tt$PValue),
mth = c(rep("scaledTPM_salmon, edgeR", nrow(res_sal_scaledTPM_edgeR$tt)),
rep("simplesum_salmon_avetxl, edgeR", nrow(res_sal_simplesum_avetxl_edgeR$tt)),
rep("simplesum_salmon, edgeR", nrow(res_sal_simplesum_edgeR$tt))))
ggplot(dfh, aes(x = pvalue)) + geom_histogram() + facet_wrap(~mth) +
plot_theme() +
xlab("p-value") + ylab("count")
par(mfrow = c(1, 3))
plotBCV(res_sal_scaledTPM_edgeR$dge, main = "scaledTPM_salmon, edgeR")
plotBCV(res_sal_simplesum_avetxl_edgeR$dge, main = "simplesum_salmon_avetxl, edgeR")
plotBCV(res_sal_simplesum_edgeR$dge, main = "simplesum_salmon, edgeR")
par(mfrow = c(1, 3))
plotSmear(res_sal_scaledTPM_edgeR$dge, main = "scaledTPM_salmon, edgeR", ylim = c(-10, 10))
plotSmear(res_sal_simplesum_avetxl_edgeR$dge, main = "simplesum_salmon_avetxl, edgeR", ylim = c(-10, 10))
plotSmear(res_sal_simplesum_edgeR$dge, main = "simplesum_salmon, edgeR", ylim = c(-10, 10))