awk -f /data/images/proton/link-bam-files-to-sample-ids1.awk ../samples.txt > link-bam-files-to-sample-ids1.sh
reczko@max:/data/images/proton/run340/www$ source !$
source link-bam-files-to-sample-ids1.sh
reczko@max:/data/images/proton/run340/www$ cp reads.txt reads.txt0
GK3R96-UI1.bam 548098
set to
GK3R96-UI1.bam 1446690
reczko@max:/data/images/proton/run340/www$ awk -f /data/images/proton/make-stranded-sampled-mm10-tracks1.awk  reads.txt  > make-stranded-sampled-mm10-tracks.sh

awk -f /data/images/proton/make-metaseqr-targets1.awk ../samples.txt > targets.txt

library(metaseqR)

the.path <- "/data/images/proton/run340/www"

the.contrasts.1 <- c(
	"UI_vs_TI",
	"UI_vs_GIH",	
	"UI_vs_GIL",
	"TI_vs_GIH",
	"TI_vs_GIL",
	"GIH_vs_GIL"
)

# 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)

# Read targets files
message("Reading targets file...")
targets <- read.targets("targets.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_run340"),
    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("natural","log2","rpgm"),
    export.values="normalized",
    export.stats="mean"
)


#@ without GIH2
# Read targets files
message("Reading targets file...")
targets <- read.targets("targets.txt2")

# 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_run340_without_GIH2"),
    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("natural","log2","rpgm"),
    export.values="normalized",
    export.stats="mean"
)


metaseqr(
    sample.list=file.path(the.path,"targets.txt2"),
    contrast=the.contrasts.1,
    annotation="download",
    org="mm10",
    refdb="ensembl",
    count.type="utr",
    normalization="deseq",
    statistics="deseq",
    fig.format=c("png","pdf"),
    export.where=file.path(the.path,"metaseqr_quantseq_run340_without_GIH2"),
    restrict.cores=0.5,
    qc.plots=c(
        "mds","biodetection","countsbio","saturation","readnoise","filtered",
        "correl","pairwise","boxplot","gcbias","lengthbias","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
        ),
        biotype=get.defaults("biotype.filter","mm10")
    ),
    pcut=0.05,
    export.what=c("annotation","p.value","adj.p.value","fold.change",
		"counts","flags"),
    export.scale=c("log2","rpgm"),
    export.values="normalized",
    export.counts.table=TRUE,
	report.top=0.05
)


#@ stranded analysis (fixing metaseqR.1.16.*'s 'bug')
require(metaseqR)
# [1] metaseqR_1.17.2
the.path <- "/data/images/proton/run340/www"

the.contrasts.1 <- c(
	"UI_vs_TI",
	"UI_vs_GIH",	
	"UI_vs_GIL",
	"TI_vs_GIH",
	"TI_vs_GIL",
	"GIH_vs_GIL"
)

metaseqr(
    sample.list=file.path(the.path,"targets.txt2"),
    contrast=the.contrasts.1,
    annotation="download",
    org="mm10",
    refdb="ensembl",
    count.type="utr",
    normalization="deseq",
    statistics="deseq",
    fig.format=c("png","pdf"),
    export.where=file.path(the.path,"metaseqr_quantseq_run340_without_GIH2_stranded"),
    restrict.cores=0.5,
    qc.plots=c(
        "mds","biodetection","countsbio","saturation","readnoise","filtered",
        "correl","pairwise","boxplot","gcbias","lengthbias","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
        ),
        biotype=get.defaults("biotype.filter","mm10")
    ),
    pcut=0.05,
    export.what=c("annotation","p.value","adj.p.value","fold.change",
		"counts","flags"),
    export.scale=c("log2","rpgm"),
    export.values="normalized",
    export.counts.table=TRUE,
	report.top=0.05
)


#@ 12092017 omit UI1
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)
targets <- read.targets("targets-woUI1.txt")
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)

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",
#             meta.p = "pandora",
#statistics = c("deseq", "edger", "noiseq", "bayseq", "limma", "nbpseq"),
#             meta.p = "pandora",
#    statistics="deseq", # or whatever supported, more than one also
    fig.format=c("png","pdf"),
     export.where=file.path(the.path,"metaseqr_quantseq_run340_omitGIH2_omitUI1_stranded"),
    restrict.cores=0.5, # fraction of available cores to use
    qc.plots=c(
        "mds","biodetection","countsbio","saturation","readnoise","filtered",
        "correl","pairwise","boxplot","gcbias","lengthbias",
        "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","meta.p.value","adj.meta.p.value","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.counts.table=TRUE,
	report.top=0.05
#    export.stats="mean"
)

#@ merge run340 and run296
GIHvsG
GILvsG
TIvsT 

the.path <- "/data/images/proton/run340/www"

the.contrasts.1 <- c(
	"GIH_vs_G",
	"GIL_vs_G",	
	"TI_vs_T"
)
the.contrasts.1 <- c(
	"G_vs_GIH",
	"G_vs_GIL",	
	"T_vs_TI"
)

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)
targets <- read.targets("targets-340-296.txt")
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)
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",
#             meta.p = "pandora",
#statistics = c("deseq", "edger", "noiseq", "bayseq", "limma", "nbpseq"),
#             meta.p = "pandora",
#    statistics="deseq", # or whatever supported, more than one also
    fig.format=c("png","pdf"),
     export.where=file.path(the.path,"metaseqr_quantseq_run340_omitGIH2_omitUI1_vs_run296_stranded"),
    restrict.cores=0.5, # fraction of available cores to use
    qc.plots=c(
        "mds","biodetection","countsbio","saturation","readnoise","filtered",
        "correl","pairwise","boxplot","gcbias","lengthbias",
        "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","meta.p.value","adj.meta.p.value","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.counts.table=TRUE,
	report.top=0.05
#    export.stats="mean"
)


#@
targets <- read.targets("targets-340-296.txt")
r2c.out <- read2count(targets,transcript.data,has.all.fields=TRUE,multic=multic)
the.counts <- cbind(r2c.out$mergedann,r2c.out$counts)


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",
#             meta.p = "pandora",
#statistics = c("deseq", "edger", "noiseq", "bayseq", "limma", "nbpseq"),
#             meta.p = "pandora",
#    statistics="deseq", # or whatever supported, more than one also
    fig.format=c("png","pdf"),
     export.where=file.path(the.path,"metaseqr_quantseq_run340_omitGIH2_omitUI1_vs_run296_stranded2"),
     restrict.cores=0.5, # fraction of available cores to use
    qc.plots=c(
        "mds","biodetection","countsbio","saturation","readnoise","filtered",
        "correl","pairwise","boxplot","gcbias","lengthbias",
        "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","meta.p.value","adj.meta.p.value","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.counts.table=TRUE,
	report.top=0.05
#    export.stats="mean"
)