RNA-seq多Run合并及VST标准化

仅限用于SLURM调配系统的服务器递交任务

#!/usr/bin/env Rscript
combine_columns <- function(df, info) {
  info$collapse <- factor(info$collapse)
  for (f in levels(info$collapse)) {
    samples <- rownames(subset(info, info$collapse == f))
    if (length(samples)==1){
      df[[f]] <- df[,samples]
      df[,samples] <- list(NULL)
    }else{
    df[[f]] <- rowSums(df[,samples])
    df[,samples] <- list(NULL)}
  }
  return(df)
}
#Combind counts files and caculate the cpm
library("argparse")
parser <- ArgumentParser(description="")
parser$add_argument('file', metavar="*.tsv", help="input gene counts")
parser$add_argument('info', metavar="*.tsv", help="sample information file; requires sample names to be in first column and experimental condition to be in 'condition' column")
parser$add_argument('-o', '--out', metavar="*", help="output file name prefix")
args <- parser$parse_args()
if (is.null(args$out)) { args$out <- sub('.counts.tsv', '', args$file) }
prefix <- args$out

countdata <- read.table(args$file, sep="\t", header=TRUE, row.names=1)
info<- read.table(args$info, sep="\t", header=TRUE,row.names =1 )
countdata_com<- combine_columns(countdata, info)
write.csv(countdata_com,paste0(prefix,"_combind_counts.csv"))
cpm <- t(t(countdata_com)/colSums(countdata_com) * 1000000)
row.names(cpm) <- row.names(countdata_com)
write.csv(cpm,paste0(prefix,"_combind_cpm.csv"))
pass <- rowSums(cpm > 1)
pass <- pass[pass > (ncol(countdata_com) / 2)]
countdata_pass <- subset(countdata_com, rownames(countdata_com) %in% names(pass))
write.csv(countdata_pass,paste0(prefix,"_combind_filter_counts.csv"))
cpm_pass <- subset(cpm, rownames(cpm) %in% names(pass))
write.csv(cpm_pass,paste0(prefix,"_combind_filter_cpm.csv"))
###############VST
suppressMessages(library(DESeq2, quietly=T))
tmpfile <- tempfile(fileext=".csv")
tmpdir <- tempdir()
vsd <- varianceStabilizingTransformation(as.matrix(countdata_pass))
write.table(vsd, file=paste0(args$out, ".vst.tsv"), sep="\t",
            quote=F, col.names=NA)
write.table(t(vsd), file=tmpfile, sep=",",
            quote=F, row.names=F, col.names=F)
command <- paste("peertool -f", tmpfile, "-n 10 -o", tmpdir)
message("Running command: ", command)
system(command)
peerCov <- read.table(file.path(tmpdir, "X.csv"), sep=",", header=F)
unlink(tmpfile)
unlink(tmpdir)
row.names(peerCov) <- colnames(vsd)
colnames(peerCov)

normalizeLM <- function(datarow, covariates) {
  datarow <- unlist(datarow)
  x <- mean(datarow)
  covariates$gene <- datarow
  fit <- lm(gene ~ ., data=covariates)
  return(residuals(fit) + x)
}

suppressMessages(library(plyr))
normalized <- adply(vsd, 1, normalizeLM, covariates=peerCov)
write.csv(normalized, file=paste0(args$out, ".vst.normalized.csv"))