--- title: "03- RNASeq md" author: "Hannia Larino" date: "`r format(Sys.time(), '%d %B, %Y')`" output: html_document: theme: readable highlight: zenburn toc: true toc_float: true number_sections: true code_folding: show code_download: true --- **HOW TO COMPLETE RNASeq Analysis** **PART I. Downloading reference genes & installing Kallisto.** #**Introduction** Using the instructions on https://sr320.github.io/course-fish546-2025/assignments/02-DGE.html, this document will show the steps required to complete and differential gene analysis from RNASeq. ##Links used: *Pacific Oyster genes fasta file: https://gannet.fish.washington.edu/seashell/bu-github/nb-2023/Cgigas/data/rna.fna * ###Libraries used: *DESeq2 *tidyverse *pheatmap *RColorBrewer *data.table *BiocManager **1. Downloading reference** ```{r, engine='bash'} cd /home/shared/8TB_HDD_02/hannia/2.RNASeq/data curl --insecure -O https://gannet.fish.washington.edu/seashell/bu-github/nb-2023/Cgigas/data/rna.fna head -n 10 <(curl -s https://gannet.fish.washington.edu/seashell/bu-github/nb-2023/Cgigas/data/rna.fna) #grabbing fasta file of pacific oyster genes. #-- insecure ignores the fact that gannet does not have security certificate to auntheticate. Not reccomended but we know server. ``` **2. Running Kallisto index command.** *Indexing file rna.fna while renaming as "cgigas_roslin_rna.index.* Pre-processing transcript sequences into format Kallisto can easily read. ```{r, engine='bash'} /home/shared/kallisto/kallisto \ index -i \ /home/shared/8TB_HDD_02/hannia/2.RNASeq/data/cgigas_roslin_rna.index \ /home/shared/8TB_HDD_02/hannia/2.RNASeq/data/rna.fna ``` **3.Downloading several data files at once listed on URL using wget.** ```{r, engine='bash'} cd /home/shared/8TB_HDD_02/hannia/2.RNASeq/data wget --recursive --no-parent --no-directories \ --no-check-certificate \ --accept '*.fastq.gz' \ https://gannet.fish.washington.edu/seashell/bu-github/nb-2023/Cgigas/data/nopp/ ``` **4. Making new folder in "output" folder.** ```{r, engine='bash'} cd /home/shared/8TB_HDD_02/hannia/2.RNASeq pwd mkdir /home/shared/8TB_HDD_02/hannia/2.RNASeq/output/kallisto_01 #mkdir makes new folder called kallisto01 up one directory level ``` **5.Creating sub directory.** ```{r, engine='bash'} find /home/shared/8TB_HDD_02/hannia/2.RNASeq/data/*fastq.gz \ | xargs basename -s _L001_R1_001.fastq.gz | xargs -I{} /home/shared/kallisto/kallisto \ quant -i /home/shared/8TB_HDD_02/hannia/2.RNASeq/data/cgigas_roslin_rna.index \ -o /home/shared/8TB_HDD_02/hannia/2.RNASeq/output/kallisto_01/{} \ -t 40 \ --single -l 100 -s 10 /home/shared/8TB_HDD_02/hannia/2.RNASeq/data/{}_L001_R1_001.fastq.gz ``` **6.This command is a bash one-liner that runs kallisto quantification in a loop over multiple single-end RNA-seq .fastq.gz files.** ```{r, engine='bash'} find /home/shared/8TB_HDD_02/hannia/2.RNASeq/data/*fastq.gz \ | xargs basename -s _L001_R1_001.fastq.gz | xargs -I{} /home/shared/kallisto/kallisto \ quant -i /home/shared/8TB_HDD_02/hannia/2.RNASeq/data/cgigas_roslin_rna.index \ -o /home/shared/8TB_HDD_02/hannia/2.RNASeq/output/kallisto_01/{} \ -t 40 \ --single -l 100 -s 10 /home/shared/8TB_HDD_02/hannia/2.RNASeq/data/{}_L001_R1_001.fastq.gz ``` **7.Creating abundance estimates to gene expression matrix.** ```{r, engine='bash'} perl /home/shared/trinityrnaseq-v2.12.0/util/abundance_estimates_to_matrix.pl \ --est_method kallisto \ --gene_trans_map none \ --out_prefix ../2.RNASeq/output/kallisto_01 \ --name_sample_by_basedir \ ../2.RNASeq/output/kallisto_01/D54_S145/abundance.tsv \ ../2.RNASeq/output/kallisto_01/D56_S136/abundance.tsv \ ../2.RNASeq/output/kallisto_01/D58_S144/abundance.tsv \ ../2.RNASeq/output/kallisto_01/M45_S140/abundance.tsv \ ../2.RNASeq/output/kallisto_01/M48_S137/abundance.tsv \ ../2.RNASeq/output/kallisto_01/M89_S138/abundance.tsv \ ../2.RNASeq/output/kallisto_01/D55_S146/abundance.tsv \ ../2.RNASeq/output/kallisto_01/D57_S143/abundance.tsv \ ../2.RNASeq/output/kallisto_01/D59_S142/abundance.tsv \ ../2.RNASeq/output/kallisto_01/M46_S141/abundance.tsv \ ../2.RNASeq/output/kallisto_01/M49_S139/abundance.tsv \ ../2.RNASeq/output/kallisto_01/M90_S147/abundance.tsv \ ../2.RNASeq/output/kallisto_01/N48_S194/abundance.tsv \ ../2.RNASeq/output/kallisto_01/N50_S187/abundance.tsv \ ../2.RNASeq/output/kallisto_01/N52_S184/abundance.tsv \ ../2.RNASeq/output/kallisto_01/N54_S193/abundance.tsv \ ../2.RNASeq/output/kallisto_01/N56_S192/abundance.tsv \ ../2.RNASeq/output/kallisto_01/N58_S195/abundance.tsv \ ../2.RNASeq/output/kallisto_01/N49_S185/abundance.tsv \ ../2.RNASeq/output/kallisto_01/N51_S186/abundance.tsv \ ../2.RNASeq/output/kallisto_01/N53_S188/abundance.tsv \ ../2.RNASeq/output/kallisto_01/N55_S190/abundance.tsv \ ../2.RNASeq/output/kallisto_01/N57_S191/abundance.tsv \ ../2.RNASeq/output/kallisto_01/N59_S189/abundance.tsv ``` **PART II. RUNNING DESeq2** **8. Reading in count matrix** \*Switching to R language and displaying first few rows. ```{r} install.packages("knitr") library(knitr) ``` ```{r, engine='r'} countmatrix <- read.delim("../2.RNASeq/output/kallisto_01.isoform.counts.matrix", header = TRUE, sep = '\t') rownames(countmatrix) <- countmatrix$X countmatrix <- countmatrix[,-1] head(countmatrix) ``` **9. Round integers up to hole numbers for further analysis.** ```{r, engine='r'} countmatrix <- round(countmatrix, 0) str(countmatrix) #like summary() ``` **10.Getting DEG based on Desication.** ```{r, engine='r'} deseq2.colData <- data.frame(condition=factor(c(rep("control", 12), rep("desicated", 12))), type=factor(rep("single-read", 24))) rownames(deseq2.colData) <- colnames(data) deseq2.dds <- DESeqDataSetFromMatrix(countData = countmatrix, colData = deseq2.colData, design = ~ condition) kable(head(deseq2.res), digits = 3, caption = "Top 6 Differential Expression Results") ``` **11. DESeq data table.** ```{r, engine='r'} deseq2.dds <- DESeq(deseq2.dds) #taking deseq data set (deseq2.dds) and runs DESeq-> this creates the desired analysis to get the table deseq2.res <- results(deseq2.dds) deseq2.res <- deseq2.res[order(rownames(deseq2.res)), ] #re-ordering alphabetically by row name head(deseq2.res) ``` \*Making adjustments. ```{r, engine='r'} # Count number of hits with adjusted p-value less then 0.05 dim(deseq2.res[!is.na(deseq2.res$padj) & deseq2.res$padj <= 0.05, ]) ``` **DESeq RESULTS** ```{r, engine='r'} tmp <- deseq2.res # The main plot plot(tmp$baseMean, tmp$log2FoldChange, pch=20, cex=0.45, ylim=c(-3, 3), log="x", col="darkgray", main="DEG Dessication (pval <= 0.05)", xlab="mean of normalized counts", ylab="Log2 Fold Change") # Getting the significant points and plotting them again so they're a different color tmp.sig <- deseq2.res[!is.na(deseq2.res$padj) & deseq2.res$padj <= 0.05, ] points(tmp.sig$baseMean, tmp.sig$log2FoldChange, pch=20, cex=0.45, col="red") # 2 FC lines abline(h=c(-1,1), col="blue") ``` #Creating heat map summary of data. #Top 50 most deferentially expressed genes. This heat map shows the relative expressions levels of the top genes. ```{r heatmap, eval=TRUE} library(DESeq2) library(pheatmap) # Select top 50 differentially expressed genes res <- results(deseq2.dds) res_ordered <- res[order(res$padj), ] top_genes <- row.names(res_ordered)[1:50] # Extract counts and normalize counts <- counts(deseq2.dds, normalized = TRUE) counts_top <- counts[top_genes, ] # Log-transform counts log_counts_top <- log2(counts_top + 1) # Generate heatmap pheatmap(log_counts_top, scale = "row") ```