R_code_deseq_te_DS3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e559c5c6-3bc2-4cc3-bbf8-3c239f0f1cc0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "library(DESeq2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef6c4b90-ad9f-4ee5-88d2-ccb824177404",
   "metadata": {},
   "source": [
    "# Methods section 2.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f560609-7561-4b95-9933-f009c874237e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# DESeq analysis for mRNA+RPF of DS3\n",
    "# split by the 3 metrics used; 60 first nt, full-length CDS, full-length CDS minus TY1-GAG\n",
    "for(metric in c('60nt','mrna','Ty1')){\n",
    "    # import count table derived from Plastid\n",
    "    count.table <- read.table(paste('~/mhenault_landrylab/md_ty_expression/deseq2/', metric, '.counts.tab', sep=''), header=TRUE, row.names='region_name')\n",
    "    #import design table\n",
    "    design.table <- read.table('~/mhenault_landrylab/md_ty_expression/deseq2/samples.tab', header=TRUE, row.names='X0', sep='\\t')\n",
    "    design.table$Experiment <- factor(design.table$Experiment, levels=c('rna','rpf'))\n",
    "    design.table$Species <- factor(design.table$Species, ordered=FALSE)\n",
    "    \n",
    "    #split the analysis per subgenome (subsampling both samples and genes accordingly)\n",
    "    for (subgenome in c('S288c','CBS432')){\n",
    "        if (subgenome == 'S288c'){\n",
    "            sub.design.table <- design.table[design.table$Species %in% c('Scer','hyb'), ]\n",
    "            sub.design.table$Species <- relevel(sub.design.table$Species, ref='Scer')\n",
    "            sub.count.table <- count.table[grep(subgenome, rownames(count.table)), rownames(sub.design.table)]\n",
    "        }\n",
    "        if (subgenome == 'CBS432'){\n",
    "            sub.design.table <- design.table[design.table$Species %in% c('Spar','hyb'), ]\n",
    "            sub.design.table$Species <- relevel(sub.design.table$Species, ref='Spar')\n",
    "            sub.count.table <- count.table[grep(subgenome, rownames(count.table)), rownames(sub.design.table)]\n",
    "        }\n",
    "\n",
    "        dat <- DESeqDataSetFromMatrix(countData=sub.count.table, colData=sub.design.table, design= ~ Experiment + Species + Experiment:Species)\n",
    "        res <- DESeq(dat)\n",
    "\n",
    "        #parse relevant contrasts and export results (raw and shrinked LFC)\n",
    "        for(f in resultsNames(res)[2:4]){\n",
    "            Res <- results(res, name=f)\n",
    "            write.csv(Res, paste('~/mhenault_landrylab/md_ty_expression/deseq2/results.',subgenome,'.',metric,'.raw_lfc.', f, '.csv', sep=''))\n",
    "\n",
    "            ResLFC <- lfcShrink(res, coef=f, type='apeglm')\n",
    "            write.csv(ResLFC, paste('~/mhenault_landrylab/md_ty_expression/deseq2/results.',subgenome,'.',metric,'.shrink_lfc.', f, '.csv', sep=''))\n",
    "        }\n",
    "        \n",
    "        # export normalized counts\n",
    "        if(metric == 'mrna'){\n",
    "            write.csv(counts(res, normalize=TRUE), paste('~/mhenault_landrylab/md_ty_expression/deseq2/norm_counts.',subgenome,'.',metric,'.csv', sep=''))\n",
    "        }\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce09c3e3-7dfe-4536-8c4c-c15377ae1756",
   "metadata": {},
   "outputs": [],
   "source": [
    "# DESeq analysis for mRNA+RPF of DS3\n",
    "# select only samples for mRNA-seq, excluding RPF data\n",
    "samples.rnaseq <- rownames(subset(design.table, Experiment=='rna'))\n",
    "# select only one metric, full-length CDS\n",
    "metric <- 'mrna'\n",
    "# import count table derived from Plastid\n",
    "count.table <- read.table(paste('~/mhenault_landrylab/md_ty_expression/deseq2/', metric, '.counts.tab', sep=''), header=TRUE, row.names='region_name')\n",
    "#import design table\n",
    "design.table <- read.table('~/mhenault_landrylab/md_ty_expression/deseq2/samples.tab', header=TRUE, row.names='X0', sep='\\t')\n",
    "design.table$Species <- factor(design.table$Species, ordered=FALSE)\n",
    "#split the analysis per subgenome (subsampling both samples and genes accordingly)\n",
    "for (subgenome in c('S288c','CBS432')){\n",
    "    if (subgenome == 'S288c'){\n",
    "        sub.design.table <- design.table[design.table$Species %in% c('Scer','hyb') & design.table$Experiment=='rna', ]\n",
    "        sub.design.table$Species <- relevel(sub.design.table$Species, ref='Scer')\n",
    "        sub.count.table <- count.table[grep(subgenome, rownames(count.table)), rownames(sub.design.table)]\n",
    "    }\n",
    "    if (subgenome == 'CBS432'){\n",
    "        sub.design.table <- design.table[design.table$Species %in% c('Spar','hyb') & design.table$Experiment=='rna', ]\n",
    "        sub.design.table$Species <- relevel(sub.design.table$Species, ref='Spar')\n",
    "        sub.count.table <- count.table[grep(subgenome, rownames(count.table)), rownames(sub.design.table)]\n",
    "    }\n",
    "    \n",
    "    dat <- DESeqDataSetFromMatrix(countData=sub.count.table, colData=sub.design.table, design= ~ Species)\n",
    "    res <- DESeq(dat)\n",
    "\n",
    "    #parse relevant contrast and export results (raw and shrinked LFC)\n",
    "    f <- resultsNames(res)[2]\n",
    "    Res <- results(res, name=f)\n",
    "    write.csv(Res, paste('~/mhenault_landrylab/md_ty_expression/deseq2/results.',subgenome,'.raw_lfc.', f, '.rnaseq.csv', sep=''))\n",
    "\n",
    "    ResLFC <- lfcShrink(res, coef=f, type='apeglm')\n",
    "    write.csv(ResLFC, paste('~/mhenault_landrylab/md_ty_expression/deseq2/results.',subgenome,'.shrink_lfc.', f, '.rnaseq.csv', sep=''))\n",
    "\n",
    "    # export normalized counts\n",
    "    write.csv(counts(res, normalize=TRUE), paste('~/mhenault_landrylab/md_ty_expression/deseq2/norm_counts.',subgenome,'.rnaseq.csv', sep=''))\n",
    "    \n",
    "}"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}