awk -f /data/images/proton/link-bam-files-to-sample-ids1.awk ../samples.txt > link-bam-files-to-sample-ids1.sh reczko@fix:/data/images/proton2/run387/www$ source !$ source link-bam-files-to-sample-ids1.sh awk -f /data/images/proton/make-metaseqr-targets1.awk ../samples.txt > targets.txt library(metaseqR) the.path <- "/data/images/proton2/run387/www" the.contrasts.1 <- c( "GFP_vs_TOM", "GFP_vs_UNS", "TOM_vs_UNS") # Read transcript data from external file transcript.data <- read.delim("/data/results/tools/rnaseq/metaseqr/mm10/transcript_data_mm10.txt.gz") rownames(transcript.data) <- as.character(transcript.data$transcript_id) # metaseqR related variables outside the pipeline multic <- check.parallel(0.5) # If wish to use multiple cores assign("VERBOSE",TRUE,envir=metaseqR:::meta.env) #@ TEST1 # Read targets files message("Reading targets file...") targets <- read.targets("targets-TEST1.txt") # Do the counting based with bam files in the targets file r2c.out <- read2count(targets,transcript.data,has.all.fields=TRUE,multic=multic) # Create a counts table to be passed as "embedded" annotation the.counts <- cbind(r2c.out$mergedann,r2c.out$counts) # Run metaseqR metaseqr( counts=the.counts, sample.list=targets$samples, contrast=the.contrasts.1, annotation="embedded", gene.file="/data/results/tools/rnaseq/metaseqr/mm10/gene_data_mm10.txt.gz", id.col=4, # required with embedded annotation gc.col=5, # required with embedded annotation bt.col=8, # required with embedded annotation name.col=7, # required with embedded annotation org="custom", count.type="utr", normalization="deseq", # or whatever supported statistics="deseq", # or whatever supported, more than one also fig.format=c("png","pdf"), export.where=file.path(the.path,"metaseqr_quantseq_run387_TEST1"), restrict.cores=0.5, # fraction of available cores to use qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise","filtered", "correl","pairwise","boxplot","gcbias","lengthbias","meandiff", "meanvar","biodist","volcano","deheatmap" ), exon.filters=NULL, gene.filters=list( length=list( length=500 ), avg.reads=list( average.per.bp=100, quantile=0.25 ), expression=list( median=TRUE, mean=FALSE, quantile=NA, known=NA, custom=NA ), # it's the default anyway biotype=get.defaults("biotype.filter","mm10") ), pcut=0.05, # only for the truncated significant list output, all results are exported anyway export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts","flags"), # if you use pandora, the fields "meta.p.value" and "adj.meta.p.value" should be added export.scale=c("log2","rpgm"), export.values="normalized", export.stats="mean" ) #TEST2 message("Reading targets file...") targets <- read.targets("targets-TEST2.txt") # Do the counting based with bam files in the targets file r2c.out <- read2count(targets,transcript.data,has.all.fields=TRUE,multic=multic) # Create a counts table to be passed as "embedded" annotation the.counts <- cbind(r2c.out$mergedann,r2c.out$counts) # Run metaseqR metaseqr( counts=the.counts, sample.list=targets$samples, contrast=the.contrasts.1, annotation="embedded", gene.file="/data/results/tools/rnaseq/metaseqr/mm10/gene_data_mm10.txt.gz", id.col=4, # required with embedded annotation gc.col=5, # required with embedded annotation bt.col=8, # required with embedded annotation name.col=7, # required with embedded annotation org="custom", count.type="utr", normalization="deseq", # or whatever supported statistics="deseq", # or whatever supported, more than one also fig.format=c("png","pdf"), export.where=file.path(the.path,"metaseqr_quantseq_run387_TEST2"), restrict.cores=0.5, # fraction of available cores to use qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise","filtered", "correl","pairwise","boxplot","gcbias","lengthbias","meandiff", "meanvar","biodist","volcano","deheatmap" ), exon.filters=NULL, gene.filters=list( length=list( length=500 ), avg.reads=list( average.per.bp=100, quantile=0.25 ), expression=list( median=TRUE, mean=FALSE, quantile=NA, known=NA, custom=NA ), # it's the default anyway biotype=get.defaults("biotype.filter","mm10") ), pcut=0.05, # only for the truncated significant list output, all results are exported anyway export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts","flags"), # if you use pandora, the fields "meta.p.value" and "adj.meta.p.value" should be added export.scale=c("log2","rpgm"), export.values="normalized", export.stats="mean" ) #@ public # Read targets files message("Reading targets file...") targets <- read.targets("targets2.txt") # Do the counting based with bam files in the targets file r2c.out <- read2count(targets,transcript.data,has.all.fields=TRUE,multic=multic) # Create a counts table to be passed as "embedded" annotation the.counts <- cbind(r2c.out$mergedann,r2c.out$counts) # Run metaseqR metaseqr( counts=the.counts, sample.list=targets$samples, contrast=the.contrasts.1, annotation="embedded", gene.file="/data/results/tools/rnaseq/metaseqr/mm10/gene_data_mm10.txt.gz", id.col=4, # required with embedded annotation gc.col=5, # required with embedded annotation bt.col=8, # required with embedded annotation name.col=7, # required with embedded annotation org="custom", count.type="utr", normalization="deseq", # or whatever supported statistics="deseq", # or whatever supported, more than one also fig.format=c("png","pdf"), export.where=file.path(the.path,"metaseqr_quantseq_run387"), restrict.cores=0.5, # fraction of available cores to use qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise","filtered", "correl","pairwise","boxplot","gcbias","lengthbias","meandiff", "meanvar","biodist","volcano","deheatmap" ), exon.filters=NULL, gene.filters=list( length=list( length=500 ), avg.reads=list( average.per.bp=100, quantile=0.25 ), expression=list( median=TRUE, mean=FALSE, quantile=NA, known=NA, custom=NA ), # it's the default anyway biotype=get.defaults("biotype.filter","mm10") ), pcut=0.05, # only for the truncated significant list output, all results are exported anyway export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts","flags"), # if you use pandora, the fields "meta.p.value" and "adj.meta.p.value" should be added export.scale=c("log2","rpgm"), export.values="normalized", export.stats="mean" )