RunDeSeq2_with_args.R

#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)

### Author - Amarinder Singh Thind - # https://github.com/amarinderthind/RNA-seq-tutorial-for-gene-differential-expression-analysis
# modifications by Kerr Wall to get to work with Dec2 dataset
# Rscript --vanilla r_arguments.r arg1 arg2

library(edgeR)
library(DESeq2)
library("biomaRt")

###################### load the raw count matrix #######################

firstC <- args[1] # N2 Spin SVIP
SecondC <- args[2] # N2 Spin SVIP
Method <- args[3] # salmon | hisat2
print(firstC)
print(SecondC)
print(Method)

outdir <- paste('.', sep='')
infile <- paste(outdir, Method, '.gene_counts.tsv', sep='')
anno_file <-'sample_info.txt'

# read input file
rawcount<-read.table(file=infile, header=TRUE, sep="\t",  row.names=1, check.names=FALSE)

# convert to integer
rawcount[, 1:ncol(rawcount)] <- sapply(rawcount[, 1:ncol(rawcount)], as.integer)

###################### Data annotation  #################################
anno <-read.table(anno_file,header=TRUE,  sep="\t") ##In this case Two coulmns (a) sample (b) Condition
rownames(anno) <- anno$sample
p.threshold <- 0.05   ##define threshold for filtering

### subset raw and conditional data for defined pairs
anno <- anno[(anno$Condition ==firstC |anno$Condition ==SecondC),]
rawcount <- rawcount[,names(rawcount) %in% anno$sample]

############################### Create DESeq2 datasets #############################
dds <- DESeqDataSetFromMatrix(countData = rawcount, colData = anno, design = ~Condition )
## Run DESEQ2
dds <- DESeq(dds)

################# contrast based  comparison ##########################
#In case of multiple comparisons ## we need to change the contrast for every comparision
contrast<- c("Condition",firstC,SecondC)
res <- results(dds, contrast=contrast)
res$threshold <- as.logical(res$padj < p.threshold)  #Threshold defined earlier
nam <- paste('down_in',firstC, sep = '_')
res[, nam] <- as.logical(res$log2FoldChange < 0)
genes.deseq <- row.names(res)[which(res$threshold)]
genes_deseq2_sig <- res[which(res$threshold),]

file <- paste(outdir, firstC,'_', SecondC, '.', Method, '.dec2_results_padj0.05.csv', sep = '')
all_results <- paste(outdir, firstC, '_', SecondC, '.', Method, '.dec2_results_all.csv', sep = '')
write.table(genes_deseq2_sig,file,sep = ",")
write.table(res,all_results,sep = ",")