--- title: "02-RNASeq" author: "Hannia Larino" date: "2025-04-10" output: html_document --- ```{bash} #working on RNA-Seq assignment. Skip instructions until running Kallisto. ``` ```{bash} pwd cd /home/shared/8TB_HDD_02/hannia pwd # put this to navigate kallisto /home/shared/kallisto/kallisto ``` Downloading reference ```{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 #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. ``` Line 30-indexing file rna.fna while renaming as "cgigas_roslin_rna.index" ```{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 ``` Line 39-downloading several data files at once listed on URL using wget ```{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/ ``` Line 46- making new folder in "output" folder ```{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 ``` line 57- creating subdirectory, view assignment instructions to annotate what each line is doing ```{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 ``` ```{bash} #mkdir /home/shared/8TB_HDD_02/hannia/2.RNASeq/output/kallisto_01 #mkdir=make new output directory, already did this 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 #line 69 Finds all .fastq.gz files in the ../data/ folder. #line 70: Extracts just the sample name by removing the suffix _L001_R1_001.fastq.gz. #So a file like SAMPLE1_L001_R1_001.fastq.gz becomes SAMPLE1. #line 70: xargs -I{}For each sample name (like SAMPLE1), runs the Kallisto quantification command. #line 71: ``` ```{bash} #erro says i need to use R, but it did download stuff cd /home/shared/8TB_HDD_02/hannia/2.RNASeq 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 ``` ```{bash} find /home/shared -name abundance_estimates_to_matrix.pl 2>/dev/null #seeing if this is installed ``` Running DESeq2 This code performs differential expression analysis to identify differentially expressed genes (DEGs) between a control condition and a desiccated condition. ```{r} install.packages("DESeq2") #might not have installed due to version of R install.packages("tidyverse") install.packages("pheatmap") install.packages("RColorBrewer") install.packages("data.table") ``` ```{r} library(DESeq2) #cant load cant download, line 142 fixing this library(tidyverse) library(pheatmap) library(RColorBrewer) library(data.table) ``` ```{r} #this package says if R is up2date install.packages("installr") library(installr) check.for.updates.R() ``` ```{r} if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("DESeq2") ``` Read in count matrix- what is count matrix? countmatrix is a transcript-level expression matrix (from Kallisto), where: Each row = one isoform (transcript). Each column = one sample. Each cell = how many reads were mapped to that isoform in that sample. ```{r} #setwd("/home/shared/8TB_HDD_02/hannia/2.RNASeq") #library("DESeq2") countmatrix <- read.delim("../2.RNASeq/output/kallisto_01.isoform.counts.matrix", header = TRUE, sep = '\t') rownames(countmatrix) <- countmatrix$X countmatrix <- countmatrix[,-1] head(countmatrix) ``` result: 2 lines D54_S145 D56_S136 D58_S144 M45_S140 M48_S137 M89_S138 D55_S146 D57_S143 D59_S142 M46_S141 M49_S139 M90_S147 XM_034472528.1 0.00000 0.0000 0.00000 0.00000 0.00000 1.00000 0.000 0.00000 0.00000 0.0000 0.00000 0.0000 XM_034461546.1 5.50031 13.0531 6.90814 11.05400 13.00000 5.04393 11.904 9.63401 20.88880 20.0601 10.75450 12.3664 Round integers up to whole numbers for further analysis: ```{r} countmatrix <- round(countmatrix, 0) str(countmatrix) #like summary() ``` Get DEGs based on Desication. Moisture removed from treatment (?) You're telling DESeq2: "Here's my count data. Here are the sample conditions (control/desicated). I want to analyze differences in expression based on condition." ```{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) ``` ```{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) ``` ```{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, ]) ``` You’re visualizing how gene expression changes due to desiccation, where: Gray dots are all genes, Red dots are significantly differentially expressed (FDR ≤ 0.05), Blue lines help you see genes that change more than 2-fold. ```{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") ``` ```{r} #getting table #since setwd() in 2.RNASeq folder, this is your relative path "output/DEGlist.tab write.table(tmp.sig, "output/DEGlist.tab", row.names = T) `` #completed assignment#