http://rseqc.sourceforge.net/#infer-experiment-py For pair-end RNA-seq, there are two different ways to strand reads (such as Illumina ScriptSeq protocol): 1++,1–,2+-,2-+ read1 mapped to ‘+’ strand indicates parental gene on ‘+’ strand read1 mapped to ‘-‘ strand indicates parental gene on ‘-‘ strand read2 mapped to ‘+’ strand indicates parental gene on ‘-‘ strand read2 mapped to ‘-‘ strand indicates parental gene on ‘+’ strand 1+-,1-+,2++,2– read1 mapped to ‘+’ strand indicates parental gene on ‘-‘ strand read1 mapped to ‘-‘ strand indicates parental gene on ‘+’ strand read2 mapped to ‘+’ strand indicates parental gene on ‘+’ strand read2 mapped to ‘-‘ strand indicates parental gene on ‘-‘ strand For single-end RNA-seq, there are also two different ways to strand reads: ++,– read mapped to ‘+’ strand indicates parental gene on ‘+’ strand read mapped to ‘-‘ strand indicates parental gene on ‘-‘ strand +-,-+ read mapped to ‘+’ strand indicates parental gene on ‘-‘ strand read mapped to ‘-‘ strand indicates parental gene on ‘+’ strand Options: --version show program’s version number and exit -h, --help show this help message and exit -i INPUT_FILE, --input-file=INPUT_FILE Input alignment file in SAM or BAM format -r REFGENE_BED, --refgene=REFGENE_BED Reference gene model in bed fomat. -s SAMPLE_SIZE, --sample-size=SAMPLE_SIZE Number of reads sampled from SAM/BAM file. default=200000 -q MAP_QUAL, --mapq=MAP_QUAL Minimum mapping quality (phred scaled) for an alignment to be considered as “uniquely mapped”. default=30 Example 1: Pair-end non strand specific: infer_experiment.py -r hg19.refseq.bed12 -i Pairend_nonStrandSpecific_36mer_Human_hg19.bam #Output:: This is PairEnd Data Fraction of reads failed to determine: 0.0172 Fraction of reads explained by "1++,1--,2+-,2-+": 0.4903 Fraction of reads explained by "1+-,1-+,2++,2--": 0.4925 Interpretation: 1.72% of total reads were mapped to genome regions that we cannot determine the “standness of transcripts” (such as regions that having both strands transcribed). For the remaining 98.28% (1 - 0.0172 = 0.9828) of reads, half can be explained by “1++,1–,2+-,2-+”, while the other half can be explained by “1+-,1-+,2++,2–”. We conclude that this is NOT a strand specific dataset because “strandness of reads” was independent of “standness of transcripts” Example 2: Pair-end strand specific: infer_experiment.py -r hg19.refseq.bed12 -i Pairend_StrandSpecific_51mer_Human_hg19.bam #Output:: This is PairEnd Data Fraction of reads failed to determine: 0.0072 Fraction of reads explained by "1++,1--,2+-,2-+": 0.9441 Fraction of reads explained by "1+-,1-+,2++,2--": 0.0487 Interpretation: 0.72% of total reads were mapped to genome regions that we cannot determine the “standness of transcripts” (such as regions that having both strands transcribed). For the remaining 99.28% (1 - 0.0072 = 0.9928) of reads, the vast majority was explained by “1++,1–,2+-,2-+”, suggesting a strand-specific dataset. Example 3: Single-end strand specific: infer_experiment.py -r hg19.refseq.bed12 -i SingleEnd_StrandSpecific_36mer_Human_hg19.bam #Output:: This is SingleEnd Data Fraction of reads failed to determine: 0.0170 Fraction of reads explained by "++,--": 0.9669 Fraction of reads explained by "+-,-+": 0.0161 Interpretation: This is single-end, strand specific RNA-seq data. Strandness of reads are concordant with strandness of reference gene. reczko@max:/data/images/proton/DKlab/orsalia$ infer_experiment.py --r /data/results/reference/mmu/mm9/mm9_Ensembl.bed -i sortedwt21pos.bam -s 999999 Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 999999 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0152 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0040 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9808 for i in *bam > do > echo $i > infer_experiment.py --r /data/results/reference/mmu/mm9/mm9_Ensembl.bed -i $i -s 1000000 > done sortedKo01pos.bam.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0559 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0072 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9369 sortedKo01pos.sorted.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0192 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0073 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9735 sortedKo02pos.bam.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0616 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0069 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9316 sortedKo02pos.sorted.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0253 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0095 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9652 sortedKo1IL4pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0235 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0076 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9689 sortedKO21pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0211 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0045 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9744 sortedKO22pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0211 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0046 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9743 sortedKo2IL4pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0233 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0069 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9698 sortedKo61pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0202 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0053 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9745 sortedKo62pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0220 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0055 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9725 sortedKOIFN1pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0178 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0050 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9772 sortedKOIFN2pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0171 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0050 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9779 sortedwt01IL4pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0216 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0077 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9707 sortedwt01pos.bam.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0581 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0064 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9355 sortedwt02IL4pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0166 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0055 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9779 sortedwt02pos.bam.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0574 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0081 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9345 sortedwt21pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0152 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0040 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9808 sortedwt22pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0162 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0043 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9795 sortedwt61pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0175 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0055 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9771 sortedwt62pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0169 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0044 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9788 sortedwtIFN1pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0144 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0048 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9808 sortedwtIFN2pos.bam Reading reference gene model /data/results/reference/mmu/mm9/mm9_Ensembl.bed ... Done Loading SAM/BAM file ... Total 1000000 usable reads were sampled This is PairEnd Data Fraction of reads failed to determine: 0.0195 Fraction of reads explained by "1++,1--,2+-,2-+": 0.0054 Fraction of reads explained by "1+-,1-+,2++,2--": 0.9751 #@ # contrast=c("WT0_vs_KO0","WT2_vs_KO2","WT6_vs_KO6","WTIFN_vs_KOIFN","WTIL4_vs_KOIL4"), reczko@max:/data/images/proton/DKlab/orsalia$ mv ~/R ~/R2 reczko@max:/data/images/proton/DKlab/orsalia$ mv ~/R.bak ~/R library(metaseqR) the.path="/data/images/proton/DKlab/orsalia" the.path="/mnt/raid/data/results/fleming/dk" metaseqr( sample.list=file.path(the.path,"targets.txt"), contrast=c("WT2_vs_KO2"), org="mm9", refdb="refseq", normalization="edger", statistics="edger", export.where=file.path(the.path,"metaseqr_IFN_IL4"), qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise", "correl","pairwise","boxplot","volcano","biodist","deheatmap" ), pcut=0.05, export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts"), export.scale=c("log2"), export.values="normalized", export.counts.table=TRUE, restrict.cores=0.3, fig.format=c("png","pdf") ) #@ on estia /home/moulos/Rbase/bin/R the.path="/mnt/raid/data/results/fleming/dk" metaseqr( sample.list=file.path(the.path,"targets.txt"), contrast=c("WT2_vs_KO2"), org="mm9", refdb="refseq", normalization="edger", statistics="edger", export.where=file.path(the.path,"metaseqr_IFN_IL4"), qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise", "correl","pairwise","boxplot","volcano","biodist","deheatmap" ), pcut=0.05, export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts"), export.scale=c("log2"), export.values="normalized", export.counts.table=TRUE, restrict.cores=0.3, fig.format=c("png","pdf") ) metaseqr( sample.list=file.path(the.path,"targets.txt"), contrast=c("WT0_vs_KO0","WT2_vs_KO2","WT6_vs_KO6","WTIFN_vs_KOIFN","WTIL4_vs_KOIL4"), org="mm9", refdb="refseq", normalization="edger", statistics="edger", export.where=file.path(the.path,"metaseqr_WT_KO_IFN_IL4"), qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise", "correl","pairwise","boxplot","volcano","biodist","deheatmap" ), pcut=0.05, export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts"), export.scale=c("log2"), export.values="normalized", export.counts.table=TRUE, restrict.cores=0.3, fig.format=c("png","pdf") ) awk -f /data/images/proton/make-stranded-sampled-mm9-tracks1.awk reads.txt > make-stranded-sampled-mm9-tracks.sh awk -f /data/images/proton/orderTracks.awk /data/images/proton/DKlab/orsalia/SStrackDb.txt > /data/images/proton/DKlab/orsalia/SStrackDb.txt2 cp SStrackDb.txt2 SStrackDb.txt rm SStrackDb.txt2 #@ merge replicates for tracks in Phub.txt awk -f /data/images/proton/make-unstranded-sampled-mm9-tracks-merging-pairs1.awk reads.txt /data/images/proton/DKlab/orsalia/make-unstranded-sampled-mm9-pair-tracks.sh normalization = c("edaseq", "deseq", "edger", "noiseq", "nbpseq", "each", "none"), meta.p="pandora", statistics=c("edaseq", "deseq", "edger", "noiseq", "nbpseq" ), #@ metaseqr( sample.list=file.path(the.path,"targets.txt-r"), contrast=c("WT0_vs_KO0","WT2_vs_KO2","WT6_vs_KO6","WTIFN_vs_KOIFN","WTIL4_vs_KOIL4"), org="mm9", refdb="refseq", normalization="edger", statistics="edger", export.where=file.path(the.path,"metaseqr_WT_KO_IFN_IL4_V2"), qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise", "correl","pairwise","boxplot","volcano","biodist","deheatmap" ), pcut=0.05, export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts"), export.scale=c("log2"), export.values="normalized", export.counts.table=TRUE, restrict.cores=0.3, fig.format=c("png","pdf") ) #@09062017 per exon the.path="/data/images/proton/DKlab/orsalia" metaseqr( sample.list=file.path(the.path,"targets.txt-r"), contrast=c("WT0_vs_KO0","WT2_vs_KO2","WT6_vs_KO6","WTIFN_vs_KOIFN","WTIL4_vs_KOIL4"), count.type="exon", org="mm9", refdb="refseq", normalization="edger", statistics="edger", export.where=file.path(the.path,"metaseqr_WT_KO_IFN_IL4_V2_exon"), qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise", "correl","pairwise","boxplot","volcano","biodist","deheatmap" ), pcut=0.05, export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts"), export.scale=c("log2"), export.values="normalized", export.counts.table=TRUE, restrict.cores=0.5, fig.format=c("png","pdf") ) # again as unstranded the.path="/data/images/proton/DKlab/orsalia" metaseqr( sample.list=file.path(the.path,"targets-unstranded.txt"), contrast=c("WT0_vs_KO0","WT2_vs_KO2","WT6_vs_KO6","WTIFN_vs_KOIFN","WTIL4_vs_KOIL4"), org="mm9", refdb="refseq", normalization="edger", statistics="edger", export.where=file.path(the.path,"metaseqr_WT_KO_IFN_IL4_V2u"), qc.plots=c( "mds","biodetection","countsbio","saturation","readnoise", "correl","pairwise","boxplot","volcano","biodist","deheatmap" ), pcut=0.05, export.what=c("annotation","p.value","adj.p.value","fold.change","stats", "counts"), export.scale=c("log2"), export.values="normalized", export.counts.table=TRUE, restrict.cores=0.5, fig.format=c("png","pdf") ) #@ /data/images/proton/DKlab/orsalia/stringtie/runStringtie1.sh SAMseq package:samr R Documentation Significance analysis of sequencing data - simple user interface Description: Correlates a large number of features (eg. genes) with an outcome variable, such as a group indicator, quantitative variable or survival time. This is a simple user interface for the samr function applied to sequencing data. For array data applications, see the function SAM. Usage: SAMseq(x, y, censoring.status = NULL, resp.type = c("Quantitative", "Two class unpaired", "Survival", "Multiclass", "Two class paired"), geneid = NULL, genenames = NULL, nperms = 100, random.seed = NULL, nresamp = 20, fdr.output = 0.20) Arguments: x: Feature matrix: p (number of features) by n (number of : samples), one observation per column (missing values allowed) y: n-vector of outcome measurements censoring.status: n-vector of censoring censoring.status (1=died or event occurred, 0=survived, or event was censored), needed for a censored survival outcome resp.type: Problem type: "Quantitative" for a continuous parameter; "Two class unpaired" for two classes with unpaired observations; "Survival" for censored survival outcome; "Multiclass": more than 2 groups; "Two class paired" for two classes with paired observations. geneid: Optional character vector of geneids for output. genenames: Optional character vector of genenames for output. nperms: Number of permutations used to estimate false discovery rates random.seed: Optional initial seed for random number generator (integer) nresamp: Number of resamples used to construct test statistic. Default : 20. fdr.output: (Approximate) False Discovery Rate cutoff for output in significant genes table Details: This is a simple, user-friendly interface to the samr package used on sequencing data. It automatically disables arguments/features that do not apply to sequencing data. It calls samr, samr.compute.delta.table and samr.compute.siggenes.table. samr detects differential expression for microarray data, and sequencing data, and other data with a large number of features. samr is the R package that is called by the "official" SAM Excel Addin. The format of the response vector y and the calling sequence is illustrated in the examples below. A more complete description is given in the SAM manual at http://www-stat.stanford.edu/~tibs/SAM Value: A list with components samr.obj: Output of samr. See documentation for samr for details : siggenes.table: Table of significant genes, output of samr.compute.siggenes.table. This has components: genes.up-matrix of significant genes having positive correlation with the outcome and genes.lo-matrix of significant genes having negative correlation with the outcome. For survival data, genes.up are those genes having positive correlation with risk- that is, increased expression corresponds to higher risk (shorter survival) genes.lo are those whose increased expression corresponds to lower risk (longer survival). delta.table: Output of samr.compute.delta.table. del: Value of delta (distance from 45 degree line in SAM plot) for used for creating delta.table and siggenes.table. Changing the input value fdr.output will change the resulting del. call: The calling sequence Author(s): Jun Li and Balasubrimanian Narasimhan and Robert Tibshirani : : References: Tusher, V., Tibshirani, R. and Chu, G. (2001): Significance analysis of microarrays applied to the ionizing radiation response PNAS 2001 98: 5116-5121, (Apr 24). http://www-stat.stanford.edu/~tibs/SAM Li, Jun and Tibshirani, R. (2011). Finding consistent patterns: a nonparametric approach for identifying differential expression in RNA-Seq data. To appear, Statistical Methods in Medical Research. Examples: ######### two class unpaired comparison set.seed(100) mu <- matrix(100, 1000, 20) mu[1:100, 11:20] <- 200 mu <- scale(mu, center=FALSE, scale=runif(20, 0.5, 1.5)) x <- matrix(rpois(length(mu), mu), 1000, 20) y <- c(rep(1, 10), rep(2, 10)) samfit <- SAMseq(x, y, resp.type = "Two class unpaired") # examine significant gene list : print(samfit) # plot results plot(samfit) ######### two class paired comparison set.seed(100) mu <- matrix(100, 1000, 20) mu[1:100, 11:20] <- 200 mu <- scale(mu, center=FALSE, scale=runif(20, 0.5, 1.5)) x <- matrix(rpois(length(mu), mu), 1000, 20) y <- c(-(1:10), 1:10) samfit <- SAMseq(x, y, resp.type = "Two class paired") # examine significant gene list print(samfit) # plot results plot(samfit) ######### Multiclass comparison set.seed(100) mu <- matrix(100, 1000, 20) : mu[1:20, 1:5] <- 120 mu[21:50, 6:10] <- 80 mu[51:70, 11:15] <- 150 mu[71:100, 16:20] <- 60 mu <- scale(mu, center=FALSE, scale=runif(20, 0.5, 1.5)) x <- matrix(rpois(length(mu), mu), 1000, 20) y <- c(rep(1:4, rep(5, 4))) samfit <- SAMseq(x, y, resp.type = "Multiclass") # examine significant gene list print(samfit) # plot results plot(samfit) ######### Quantitative comparison set.seed(100) mu <- matrix(100, 1000, 20) y <- runif(20, 1, 3) mu[1 : 100, ] <- matrix(rep(100 * y, 100), ncol=20, byrow=TRUE) mu <- scale(mu, center=FALSE, scale=runif(20, 0.5, 1.5)) x <- matrix(rpois(length(mu), mu), 1000, 20) samfit <- SAMseq(x, y, resp.type = "Quantitative") : # examine significant gene list print(samfit) # plot results plot(samfit) ######### Survival comparison set.seed(100) mu <- matrix(100, 1000, 20) y <- runif(20, 1, 3) mu[1 : 100, ] <- matrix(rep(100 * y, 100), ncol=20, byrow=TRUE) mu <- scale(mu, center=FALSE, scale=runif(20, 0.5, 1.5)) x <- matrix(rpois(length(mu), mu), 1000, 20) y <- y + runif(20, -0.5, 0.5) censoring.status <- as.numeric(y < 2.3) y[y >= 2.3] <- 2.3 samfit <- SAMseq(x, y, censoring.status = censoring.status, resp.type = "Survival") # examine significant gene list print(samfit) # plot results : plot(samfit) #@ "Elavl1","MSTRG.14897","transcript","93200",3.40392012592381,0.101864713739399,0.898647182772707 chr8 4311599 4311786 ENSMUSG00000040028_MEX_5 0 - Elavl1 protein_coding (>2fold down in KO) /data/images/proton/DKlab/orsalia/stringtie/comare1.annotated.gtf chr8 StringTie transcript 4284782 4325100 . - . transcript_id "ENSMUST00000098950"; gene_id "MSTRG.14897"; gene_name "Elavl1"; xloc "XLOC_032330"; ref_gene_id "ENSMUSG00000040028"; cmp_ref "ENSMUST00000098950"; class_code "="; tss_id "TSS61333"; chr8 StringTie exon 4284782 4289924 . - . transcript_id "ENSMUST00000098950"; gene_id "MSTRG.14897"; exon_number "1"; chr8 StringTie exon 4295336 4295561 . - . transcript_id "ENSMUST00000098950"; gene_id "MSTRG.14897"; exon_number "2"; chr8 StringTie exon 4301685 4301838 . - . transcript_id "ENSMUST00000098950"; gene_id "MSTRG.14897"; exon_number "3"; chr8 StringTie exon 4304927 4305030 . - . transcript_id "ENSMUST00000098950"; gene_id "MSTRG.14897"; exon_number "4"; chr8 StringTie exon 4311599 4311786 . - . transcript_id "ENSMUST00000098950"; gene_id "MSTRG.14897"; exon_number "5"; chr8 StringTie exon 4324886 4325100 . - . transcript_id "ENSMUST00000098950"; gene_id "MSTRG.14897"; exon_number "6"; chr8 StringTie transcript 4288342 4325100 . - . transcript_id "MSTRG.14897.2"; gene_id "MSTRG.14897"; gene_name "Elavl1"; xloc "XLOC_032330"; cmp_ref "ENSMUST00000098950"; class_code "j"; tss_id "TSS61333"; chr8 StringTie exon 4288342 4289924 . - . transcript_id "MSTRG.14897.2"; gene_id "MSTRG.14897"; exon_number "1"; chr8 StringTie exon 4295336 4295561 . - . transcript_id "MSTRG.14897.2"; gene_id "MSTRG.14897"; exon_number "2"; chr8 StringTie exon 4301685 4301838 . - . transcript_id "MSTRG.14897.2"; gene_id "MSTRG.14897"; exon_number "3"; chr8 StringTie exon 4304927 4305030 . - . transcript_id "MSTRG.14897.2"; gene_id "MSTRG.14897"; exon_number "4"; chr8 StringTie exon 4311599 4311786 . - . transcript_id "MSTRG.14897.2"; gene_id "MSTRG.14897"; exon_number "5"; chr8 StringTie exon 4312557 4312614 . - . transcript_id "MSTRG.14897.2"; gene_id "MSTRG.14897"; exon_number "6"; chr8 StringTie exon 4324886 4325100 . - . transcript_id "MSTRG.14897.2"; gene_id "MSTRG.14897"; exon_number "7"; chr8 StringTie transcript 4288345 4325100 . - . transcript_id "MSTRG.14897.3"; gene_id "MSTRG.14897"; gene_name "Elavl1"; xloc "XLOC_032330"; cmp_ref "ENSMUST00000098950"; class_code "j"; tss_id "TSS61333"; chr8 StringTie exon 4288345 4289924 . - . transcript_id "MSTRG.14897.3"; gene_id "MSTRG.14897"; exon_number "1"; chr8 StringTie exon 4295336 4295561 . - . transcript_id "MSTRG.14897.3"; gene_id "MSTRG.14897"; exon_number "2"; chr8 StringTie exon 4301685 4301838 . - . transcript_id "MSTRG.14897.3"; gene_id "MSTRG.14897"; exon_number "3"; chr8 StringTie exon 4304927 4305030 . - . transcript_id "MSTRG.14897.3"; gene_id "MSTRG.14897"; exon_number "4"; chr8 StringTie exon 4324886 4325100 . - . transcript_id "MSTRG.14897.3"; gene_id "MSTRG.14897"; exon_number "5"; chr8 StringTie transcript 4288346 4324980 . - . transcript_id "MSTRG.14897.4"; gene_id "MSTRG.14897"; gene_name "Elavl1"; xloc "XLOC_032330"; cmp_ref "ENSMUST00000098950"; class_code "j"; tss_id "TSS61334"; chr8 StringTie exon 4288346 4289924 . - . transcript_id "MSTRG.14897.4"; gene_id "MSTRG.14897"; exon_number "1"; chr8 StringTie exon 4295336 4295561 . - . transcript_id "MSTRG.14897.4"; gene_id "MSTRG.14897"; exon_number "2"; chr8 StringTie exon 4301685 4301838 . - . transcript_id "MSTRG.14897.4"; gene_id "MSTRG.14897"; exon_number "3"; chr8 StringTie exon 4304927 4305030 . - . transcript_id "MSTRG.14897.4"; gene_id "MSTRG.14897"; exon_number "4"; chr8 StringTie exon 4311599 4311786 . - . transcript_id "MSTRG.14897.4"; gene_id "MSTRG.14897"; exon_number "5"; chr8 StringTie exon 4322703 4322832 . - . transcript_id "MSTRG.14897.4"; gene_id "MSTRG.14897"; exon_number "6"; chr8 StringTie exon 4324886 4324980 . - . transcript_id "MSTRG.14897.4"; gene_id "MSTRG.14897"; exon_number "7";