StringTie v1.3.3b usage:
 stringtie <input.bam ..> [-G <guide_gff>] [-l <label>] [-o <out_gtf>] [-p <cpus>]
  [-v] [-a <min_anchor_len>] [-m <min_tlen>] [-j <min_anchor_cov>] [-f <min_iso>]
  [-C <coverage_file_name>] [-c <min_bundle_cov>] [-g <bdist>] [-u]
  [-e] [-x <seqid,..>] [-A <gene_abund.out>] [-h] {-B | -b <dir_path>} 
Assemble RNA-Seq alignments into potential transcripts.
 Options:
 --version : print just the version at stdout and exit
 -G reference annotation to use for guiding the assembly process (GTF/GFF3)
 --rf assume stranded library fr-firststrand
 --fr assume stranded library fr-secondstrand
 -l name prefix for output transcripts (default: STRG)
 -f minimum isoform fraction (default: 0.1)
 -m minimum assembled transcript length (default: 200)
 -o output path/file name for the assembled transcripts GTF (default: stdout)
 -a minimum anchor length for junctions (default: 10)
 -j minimum junction coverage (default: 1)
 -t disable trimming of predicted transcripts based on coverage
    (default: coverage trimming is enabled)
 -c minimum reads per bp coverage to consider for transcript assembly
    (default: 2.5)
 -v verbose (log bundle processing details)
 -g gap between read mappings triggering a new bundle (default: 50)
 -C output a file with reference transcripts that are covered by reads
 -M fraction of bundle allowed to be covered by multi-hit reads (default:0.95)
 -p number of threads (CPUs) to use (default: 1)
 -A gene abundance estimation output file
 -B enable output of Ballgown table files which will be created in the
    same directory as the output GTF (requires -G, -o recommended)
 -b enable output of Ballgown table files but these files will be 
    created under the directory path given as <dir_path>
 -e only estimate the abundance of given reference transcripts (requires -G)
 -x do not assemble any transcripts on the given reference sequence(s)
 -u no multi-mapping correction (default: correction enabled)
 -h print this usage message and exit

Transcript merge usage mode: 
  stringtie --merge [Options] { gtf_list | strg1.gtf ...}
With this option StringTie will assemble transcripts from multiple
input files generating a unified non-redundant set of isoforms. In this mode
the following options are available:
  -G <guide_gff>   reference annotation to include in the merging (GTF/GFF3)
  -o <out_gtf>     output file name for the merged transcripts GTF
                    (default: stdout)
  -m <min_len>     minimum input transcript length to include in the merge
                    (default: 50)
  -c <min_cov>     minimum input transcript coverage to include in the merge
                    (default: 0)
  -F <min_fpkm>    minimum input transcript FPKM to include in the merge
                    (default: 1.0)
  -T <min_tpm>     minimum input transcript TPM to include in the merge
                    (default: 1.0)
  -f <min_iso>     minimum isoform fraction (default: 0.01)
  -g <gap_len>     gap between transcripts to merge together (default: 250)
  -i               keep merged transcripts with retained introns; by default
                   these are not kept unless there is strong evidence for them
  -l <label>       name prefix for output transcripts (default: MSTRG)

  
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt01pos.sorted.bam -p 8 -o wt01.gtf2 -A wt01.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> wt01.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt02pos.sorted.bam -p 8 -o wt02.gtf2 -A wt02.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> wt02.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt21pos.bam -p 8 -o wt21.gtf2 -A wt21.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> wt21.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt22pos.bam -p 8 -o wt22.gtf2 -A wt22.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> wt22.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt61pos.bam -p 8 -o wt61.gtf2 -A wt61.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> wt61.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt62pos.bam -p 8 -o wt62.gtf2 -A wt62.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> wt62.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwtIFN1pos.bam -p 8 -o wtIFN1.gtf2 -A wtIFN1.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> wtIFN1.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwtIFN2pos.bam -p 8 -o wtIFN2.gtf2 -A wtIFN2.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> wtIFN2.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt01IL4pos.bam -p 8 -o wt01IL4.gtf2 -A wt01IL4.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> wt01IL4.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt02IL4pos.bam -p 8 -o wt02IL4.gtf2 -A wt02IL4.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> wt02IL4.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKo01pos.sorted.bam -p 8 -o Ko01.gtf2 -A Ko01.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> Ko01.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKo02pos.sorted.bam -p 8 -o Ko02.gtf2 -A Ko02.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> Ko02.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKO21pos.bam -p 8 -o KO21.gtf2 -A KO21.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> KO21.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKO22pos.bam -p 8 -o KO22.gtf2 -A KO22.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> KO22.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKo61pos.bam -p 8 -o Ko61.gtf2 -A Ko61.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> Ko61.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKo62pos.bam -p 8 -o Ko62.gtf2 -A Ko62.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> Ko62.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKOIFN1pos.bam -p 8 -o KOIFN1.gtf2 -A KOIFN1.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> KOIFN1.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKOIFN2pos.bam -p 8 -o KOIFN2.gtf2 -A KOIFN2.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> KOIFN2.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKo1IL4pos.bam -p 8 -o Ko1IL4.gtf2 -A Ko1IL4.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> Ko1IL4.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKo2IL4pos.bam -p 8 -o Ko2IL4.gtf2 -A Ko2IL4.tab -B -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf &> Ko2IL4.log


/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie --merge -o merged.gtf2 -G /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf Ko01.gtf2  Ko1IL4.gtf2  KO22.gtf2    Ko61.gtf2  KOIFN1.gtf2  wt01.gtf2     wt02.gtf2     wt21.gtf2  wt61.gtf2  wtIFN1.gtf2 Ko02.gtf2  KO21.gtf2    Ko2IL4.gtf2  Ko62.gtf2  KOIFN2.gtf2  wt01IL4.gtf2  wt02IL4.gtf2  wt22.gtf2  wt62.gtf2  wtIFN2.gtf2
/data/results/tools/rnaseq/stringtie/gffcompare-0.10.1.Linux_x86_64/gffcompare -r /data/images/proton/DKlab/orsalia2/stringtie/Mus_musculus.NCBIM37.64-toMM9-yellow.gtf -o comare1 merged.gtf2

/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt01pos.sorted.bam -p 8 -e -B -o wt01/wt01.merged.gtf2  -G merged.gtf2 &> wt01.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt02pos.sorted.bam -p 8 -e -B -o wt02/wt02.merged.gtf2  -G merged.gtf2 &> wt02.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt21pos.bam -p 8 -e -B -o wt21/wt21.merged.gtf2  -G merged.gtf2 &> wt21.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt22pos.bam -p 8 -e -B -o wt22/wt22.merged.gtf2  -G merged.gtf2 &> wt22.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt61pos.bam -p 8 -e -B -o wt61/wt61.merged.gtf2  -G merged.gtf2 &> wt61.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt62pos.bam -p 8 -e -B -o wt62/wt62.merged.gtf2  -G merged.gtf2 &> wt62.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwtIFN1pos.bam -p 8 -e -B -o wtIFN1/wtIFN1.merged.gtf2  -G merged.gtf2 &> wtIFN1.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwtIFN2pos.bam -p 8 -e -B -o wtIFN2/wtIFN2.merged.gtf2  -G merged.gtf2 &> wtIFN2.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt01IL4pos.bam -p 8 -e -B -o wt01IL4/wt01IL4.merged.gtf2  -G merged.gtf2 &> wt01IL4.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedwt02IL4pos.bam -p 8 -e -B -o wt02IL4/wt02IL4.merged.gtf2  -G merged.gtf2 &> wt02IL4.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKo01pos.sorted.bam -p 8 -e -B -o Ko01/Ko01.merged.gtf2  -G merged.gtf2 &> Ko01.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKo02pos.sorted.bam -p 8 -e -B -o Ko02/Ko02.merged.gtf2  -G merged.gtf2 &> Ko02.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKO21pos.bam -p 8 -e -B -o KO21/KO21.merged.gtf2  -G merged.gtf2 &> KO21.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKO22pos.bam -p 8 -e -B -o KO22/KO22.merged.gtf2  -G merged.gtf2 &> KO22.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKo61pos.bam -p 8 -e -B -o Ko61/Ko61.merged.gtf2  -G merged.gtf2 &> Ko61.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKo62pos.bam -p 8 -e -B -o Ko62/Ko62.merged.gtf2  -G merged.gtf2 &> Ko62.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKOIFN1pos.bam -p 8 -e -B -o KOIFN1/KOIFN1.merged.gtf2  -G merged.gtf2 &> KOIFN1.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKOIFN2pos.bam -p 8 -e -B -o KOIFN2/KOIFN2.merged.gtf2  -G merged.gtf2 &> KOIFN2.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKo1IL4pos.bam -p 8 -e -B -o Ko1IL4/Ko1IL4.merged.gtf2  -G merged.gtf2 &> Ko1IL4.log
/data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie /data/images/proton/DKlab/orsalia/sortedKo2IL4pos.bam -p 8 -e -B -o Ko2IL4/Ko2IL4.merged.gtf2  -G merged.gtf2 &> Ko2IL4.log

mkdir ballgown2

mv Ko01 ballgown2
      mv KO21 ballgown2
      mv Ko61 ballgown2
mv KOIFN2 ballgown2
mv wt02IL4 ballgown2
mv wt22 ballgown2
mv wtIFN1 ballgown2

mv Ko02 ballgown2
mv KO22 ballgown2
mv Ko62 ballgown2
mv wt01IL4 ballgown2
mv wt02 ballgown2
mv wt61 ballgown2
mv wtIFN2 ballgown2

mv Ko1IL4 ballgown2
mv Ko2IL4 ballgown2
mv KOIFN1 ballgown2
mv wt01 ballgown2
mv wt21 ballgown2
mv wt62 ballgown2

(change .gtf2 to .gtf)
reczko@max:/data/images/proton/DKlab/orsalia2/stringtie$ /data/results/tools/rnaseq/stringtie/stringtie-master/prepDE.py -i ballgown2


library(ballgown)
library(RSkittleBrewer)
library(genefilter)
library(dplyr)
library(devtools)
pheno_data = read.csv("DKphenodata.csv") #pheno_data = read.csv("geuvadis_phenodata.csv")
samples=c("./ballgown2/Ko01","./ballgown2/Ko02","./ballgown2/KO21","./ballgown2/KO22","./ballgown2/Ko61","./ballgown2/Ko62","./ballgown2/KOIFN1","./ballgown2/KOIFN2","./ballgown2/Ko1IL4","./ballgown2/Ko2IL4","./ballgown2/wt01","./ballgown2/wt02","./ballgown2/wt21","./ballgown2/wt22","./ballgown2/wt61","./ballgown2/wt62","./ballgown2/wtIFN1","./ballgown2/wtIFN2","./ballgown2/wt01IL4","./ballgown2/wt02IL4"	 )

#@ save all:
#bg_chrX = ballgown(samples=samples,  pData=pheno_data) #bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
bg_chrX = ballgown(samples=samples,  pData=pheno_data) #bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
sampleNames(bg_chrX)
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)
full_table <- texpr(bg_chrX , 'all')
write.csv(full_table, "stringtie_transcript_results2.csv",row.names=FALSE)

#@ all wt vs all ko
pheno_data = read.csv("DKphenodata-wtVSko.csv") #pheno_data = read.csv("geuvadis_phenodata.csv")
bg_chrX = ballgown(samples=samples,  pData=pheno_data) #bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)
results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="cond", getFC=TRUE, meas="FPKM")
#results_genes = stattest(bg_chrX_filt, feature="gene",covariate="sex", adjustvars = c("population"), getFC=TRUE,meas="FPKM")
results_genes = stattest(bg_chrX_filt, feature="gene",covariate="cond", getFC=TRUE,meas="FPKM")
results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)
results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)
write.csv(results_transcripts, "wtVSko_transcript_results.csv",row.names=FALSE)
write.csv(results_genes, "wtVSko_gene_results.csv", row.names=FALSE)

require(NBPseq)

pheno_data = read.csv("DKphenodata-ko0VSwt0.csv") #pheno_data = read.csv("geuvadis_phenodata.csv")=

pheno_data = read.table("DKphenodata-ko0VSwt0.csv",header=T,sep=",")
str(pheno_data)
#'data.frame':	4 obs. of  2 variables:
# $ ids : Factor w/ 4 levels "Ko01","Ko02",..: 1 2 3 4
# $ cond: Factor w/ 2 levels "ko0h","wt0h": 1 1 2 2
samples=c("./ballgown/Ko01","./ballgown/Ko02","./ballgown/wt01","./ballgown/wt02")

bg_chrX = ballgown(samples=samples,  pData=pheno_data) #bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
sampleNames(bg_chrX)
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)
#results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="sex",adjustvars =c("population"), getFC=TRUE, meas="FPKM")
results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="cond", getFC=TRUE, meas="FPKM")
#results_genes = stattest(bg_chrX_filt, feature="gene",covariate="sex", adjustvars = c("population"), getFC=TRUE,meas="FPKM")
results_genes = stattest(bg_chrX_filt, feature="gene",covariate="cond", getFC=TRUE,meas="FPKM")
results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)
results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)
write.csv(results_transcripts, "wt0VSko0_transcript_results.csv",row.names=FALSE)
write.csv(results_genes, "wt0VSko0_gene_results.csv", row.names=FALSE)

pheno_data = read.table("DKphenodata-ko2VSwt2.csv",header=T,sep=",")
samples=c("./ballgown/KO21","./ballgown/KO22","./ballgown/wt21","./ballgown/wt22")

bg_chrX = ballgown(samples=samples,  pData=pheno_data) #bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)
#results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="sex",adjustvars =c("population"), getFC=TRUE, meas="FPKM")
results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="cond", getFC=TRUE, meas="FPKM")
#results_genes = stattest(bg_chrX_filt, feature="gene",covariate="sex", adjustvars = c("population"), getFC=TRUE,meas="FPKM")
results_genes = stattest(bg_chrX_filt, feature="gene",covariate="cond", getFC=TRUE,meas="FPKM")
results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)
results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)
write.csv(results_transcripts, "wt2VSko2_transcript_results.csv",row.names=FALSE)
write.csv(results_genes, "wt2VSko2_gene_results.csv", row.names=FALSE)

pheno_data = read.table("DKphenodata-ko6VSwt6.csv",header=T,sep=",")
samples=c("./ballgown/Ko61","./ballgown/Ko62","./ballgown/wt61","./ballgown/wt62")
bg_chrX = ballgown(samples=samples,  pData=pheno_data) #bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)
#results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="sex",adjustvars =c("population"), getFC=TRUE, meas="FPKM")
results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="cond", getFC=TRUE, meas="FPKM")
#results_genes = stattest(bg_chrX_filt, feature="gene",covariate="sex", adjustvars = c("population"), getFC=TRUE,meas="FPKM")
results_genes = stattest(bg_chrX_filt, feature="gene",covariate="cond", getFC=TRUE,meas="FPKM")
results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)
results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)
write.csv(results_transcripts, "wt6VSko6_transcript_results.csv",row.names=FALSE)
write.csv(results_genes, "wt6VSko6_gene_results.csv", row.names=FALSE)

pheno_data = read.table("DKphenodata-koIFN-VS-wtIFN.csv",header=T,sep=",")
samples=c("./ballgown/KOIFN1","./ballgown/KOIFN2","./ballgown/wtIFN1","./ballgown/wtIFN2")
bg_chrX = ballgown(samples=samples,  pData=pheno_data) #bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)
#results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="sex",adjustvars =c("population"), getFC=TRUE, meas="FPKM")
results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="cond", getFC=TRUE, meas="FPKM")
#results_genes = stattest(bg_chrX_filt, feature="gene",covariate="sex", adjustvars = c("population"), getFC=TRUE,meas="FPKM")
results_genes = stattest(bg_chrX_filt, feature="gene",covariate="cond", getFC=TRUE,meas="FPKM")
results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)
results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)
write.csv(results_transcripts, "wtIFN-VS-koIFN_transcript_results.csv",row.names=FALSE)
write.csv(results_genes, "wtIFN-VS-koIFN_gene_results.csv", row.names=FALSE)

pheno_data = read.table("DKphenodata-koIL4-VS-wtIL4.csv",header=T,sep=",")
samples=c("./ballgown/Ko1IL4","./ballgown/Ko2IL4","./ballgown/wt01IL4","./ballgown/wt02IL4")
bg_chrX = ballgown(samples=samples,  pData=pheno_data) #bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)
#results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="sex",adjustvars =c("population"), getFC=TRUE, meas="FPKM")
results_transcripts = stattest(bg_chrX_filt,feature="transcript",covariate="cond", getFC=TRUE, meas="FPKM")
#results_genes = stattest(bg_chrX_filt, feature="gene",covariate="sex", adjustvars = c("population"), getFC=TRUE,meas="FPKM")
results_genes = stattest(bg_chrX_filt, feature="gene",covariate="cond", getFC=TRUE,meas="FPKM")
results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)
results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)
write.csv(results_transcripts, "wtIL4-VS-koIL4_transcript_results.csv",row.names=FALSE)
write.csv(results_genes, "wtIL4-VS-koIL4_gene_results.csv", row.names=FALSE)

#see getAltspliced.sh


awk -f mergeDEresults1.awk all_pairs.txt > wtVSko_transcript_results_all_pairs.csv
awk -f mergeDEresultsAndGetAltSpliced1.awk all_pairs.txt > wtVSko_transcript_results_all_pairs_alt_spliced.csv


tropical= c('darkorange', 'dodgerblue','hotpink', 'limegreen', 'yellow')
palette(tropical)

/data/images/proton/DKlab/orsalia/stringtie/merged.gtf2
chr8	StringTie	transcript	4284782	4325100	1000	-	.	gene_id "MSTRG.14897"; transcript_id "ENSMUST00000098950"; gene_name "Elavl1"; ref_gene_id "ENSMUSG00000040028"; 

match( "MSTRG.14897" ,ballgown::geneIDs(bg_chrX))
[1] 93200

png("Isoform expression in Elavl1-wt0h-vs-ko0h_2.png",width=1024,height=1200,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[93200], bg_chrX, main=c('Isoform expression in Elavl1'),samples=c( "Ko01","Ko02","wt01","wt02"))
dev.off()

png("Isoform expression in Dlst-wt0h-vs-ko0h_2.png",width=1024,height=1200,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.3721" ,ballgown::geneIDs(bg_chrX))], bg_chrX, main=c('Isoform expression in Dlst'),samples=c( "Ko01","Ko02","wt01","wt02"))
dev.off()

 png("Isoform expression in Elavl1_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[93200], bg_chrX, main=c('Isoform expression in Elavl1'),samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()

png("Isoform expression in Dlst_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.3721" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()

png("Isoform expression in Rnf6_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.12773" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Csde1_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.10262" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Serinc3_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.9564" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Ddx58_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.10670" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Ddx23_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.5697" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Lgals8_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.3950" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Pgap2_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.14430" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Srp9_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.933" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Fnbp4_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.8943" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Slc39a1_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.10110" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Srek1_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.4490" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Bmp2k_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.12191" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Smndc1_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.8260" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Cast_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.4357" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Fxyd2_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.16091" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Mbnl1_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.9921" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Dctn2_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.1740" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("Isoform expression in Atpaf2_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.2268" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()



fpkm = texpr(bg_chrX,meas="FPKM")
fpkm = log2(fpkm+1)
boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')
plot(fpkm[12,] ~ pheno_data$sex, border=c(1,2),main=paste(ballgown::geneNames(bg_chrX)[12],' : ',ballgown::transcriptNames(bg_chrX)[12]),pch=19, xlab="Sex",ylab='log2(FPKM+1)')
points(fpkm[12,] ~ jitter(as.numeric(pheno_data$sex)),col=as.numeric(pheno_data$sex))
plotTranscripts(ballgown::geneIDs(bg_chrX)[1729], bg_chrX, main=c('Gene XIST in sample ERR188234'), sample=c('ERR188234'))
plotMeans('MSTRG.56', bg_chrX_filt,groupvar="sex",legend=FALSE)

library(samr)
colnames(rnaseqMatrix)
rnaseqMatrix= read.table("transcript_count_matrix.csv",header=T,sep=",")
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
x=rnaseqMatrix[,c(5,6,11,13)]
y=c(2,2,1,1)
#     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)

samfit <- SAMseq(x, y, resp.type = "Two class unpaired", geneid=row.names(x), fdr=0.05,nperms=200,nresamp=40)
samfit$siggenes.table$ngenes.up
[1] 7667
Genes up
        Gene ID Gene Name          Score(d) Fold Change      q-value(%)
   [1,] g12     ENSMUST00000135559 2        1.533            84.022    
   [2,] g26     MSTRG.13892.3      2        1.539            84.022    
samfit <- SAMseq(x, y, resp.type = "Two class unpaired", geneid=row.names(x), fdr=0.05,nperms=200,nresamp=1000)
Genes up
        Gene ID Gene Name          Score(d) Fold Change      q-value(%)
   [1,] g12     ENSMUST00000135559 2        1.533            84.171    
   [2,] g26     MSTRG.13892.3      2        1.539            84.171    

samfit <- SAMseq(x, y, resp.type = "Two class unpaired", geneid=row.names(x), fdr=0.05,nperms=200,nresamp=5000)
samfit$siggenes.table$ngenes.up
[1] 8241
Number of thresholds chosen (all possible thresholds) = 2874
Genes up
        Gene ID Gene Name          Score(d) Fold Change      q-value(%)
   [1,] g44     ENSMUST00000019026 2        1.506            84.415    
   [2,] g55     MSTRG.10069.6      2        4.432            84.415    

Number of thresholds chosen (all possible thresholds) = 3932
Genes up
        Gene ID Gene Name          Score(d) Fold Change      q-value(%)
   [1,] g44     ENSMUST00000019026 2        1.506            84.617    
   [2,] g55     MSTRG.10069.6      2        4.432            84.617    
SAMseq(x = x, y = y, resp.type = "Two class unpaired", geneid = row.names(x), 
    nperms = 200, nresamp = 20000, fdr.output = 0.05)
   
samfit <- SAMseq(x, y, resp.type = "Two class unpaired", geneid=row.names(x), fdr=0.05)
 print(samfit)
#Genes up
#        Gene ID Gene Name          Score(d) Fold Change      q-value(%)
#   [1,] g9      ENSMUST00000135550 2        7.98             82.244    
#   [2,] g12     ENSMUST00000135559 2        1.533            82.244    
#   [3,] g25     MSTRG.8388.6       2        2.671            82.244    
plotTranscripts(ballgown::geneIDs(bg_chrX)[1729], bg_chrX, main=c('Gene XIST in sample ERR188234'), sample=c('ERR188234'))
pheno_data = read.csv("geuvadis_phenodata.csv")
bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)


require(edgeR)
rnaseqMatrix0= read.table("transcript_count_matrix.csv",header=T,sep=",")
#data = read.table("/data/images/proton/run125/www/genes.counts.matrix", header=T, row.names=1, com='')
#col_ordering = c(1,2,6,7)

rnaseqMatrix =rnaseqMatrix0;
#rnaseqMatrix = round(rnaseqMatrix)
rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,]
#conditions = factor(c(rep("FAGsA", 2), rep("FAGsB", 2)))
#                       1    2     3       4       5     6     7       8       9     10    11    12       13    14       15    16    17    18    19      20 
#conditions = factor(c("KO2","KO2","KOIFN","KOIFN","Ko0","Ko0","KoIL4","KoIL4","Ko6","Ko6","wt0","wt0IL4","wt0","wt0IL4","wt2","wt2","wt6","wt6","wtIFN","wtIFN"))
conditions = factor(c("KO","KO","KO","KO","KO","KO","KO","KO","KO","KO","wt","wt","wt","wt","wt","wt","wt","wt","wt","wt"))

exp_study = DGEList(counts=rnaseqMatrix, group=conditions)
exp_study = calcNormFactors(exp_study)
exp_study = estimateCommonDisp(exp_study)
exp_study = estimateTagwiseDisp(exp_study)
et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
r = rownames(tTags)
c=data[r,col_ordering]
#tab <- topTags(lrt)
#write.table(tTags, file="edgeRtranscripts.csv")

col_ordering = c(11,13,5,6)#0h
rnaseqMatrix =rnaseqMatrix0[,col_ordering];rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];conditions = factor(c("WT","WT","KO","KO"));
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
r <- rownames(topTags(et)$table)
c=cpm(exp_study)[r, order(exp_study$samples$group)]
x=cbind(as.data.frame(tTags),c)
summary(de <- decideTestsDGE(et, p=0.05))
#detags <- rownames(exp_study)[as.logical(de)]
#plotSmear(et, de.tags=detags)
#abline(h = c(-1, 1), col = "blue")
write.table(tTags, file="edgeRtranscripts-0h.csv")
write.table(x, file="edgeRtranscriptsExpr-0h.csv")

col_ordering = c(15,16,1,2)#2h
rnaseqMatrix =rnaseqMatrix0[,col_ordering];rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];conditions = factor(c("WT","WT","KO","KO"));
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
write.table(tTags, file="edgeRtranscripts-2h.csv")
col_ordering = c(17,18,9,10)#6h
rnaseqMatrix =rnaseqMatrix0[,col_ordering];rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];conditions = factor(c("WT","WT","KO","KO"));
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
write.table(tTags, file="edgeRtranscripts-6h.csv")
col_ordering = c(19,20,3,4)#IFN
rnaseqMatrix =rnaseqMatrix0[,col_ordering];rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];conditions = factor(c("WT","WT","KO","KO"));
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
write.table(tTags, file="edgeRtranscripts-IFN.csv")
col_ordering = c(12,14,7,8)#IL4
rnaseqMatrix =rnaseqMatrix0[,col_ordering];rnaseqMatrix = rnaseqMatrix[rowSums(rnaseqMatrix)>=2,];conditions = factor(c("WT","WT","KO","KO"));
exp_study = DGEList(counts=rnaseqMatrix, group=conditions);exp_study = calcNormFactors(exp_study);exp_study = estimateCommonDisp(exp_study);exp_study = estimateTagwiseDisp(exp_study);et = exactTest(exp_study)
tTags = topTags(et,n=NULL,adjust.method = "fdr")
write.table(tTags, file="edgeRtranscripts-IL4.csv")

awk -f mergeEdgerDEresults1.awk all_pairs_edger.txt > edgeRtranscripts_all_pairs.csv
ref transcripts 74118
transcripts in edgeRtranscripts-0h : 65358
transcripts in edgeRtranscripts-2h : 63764
transcripts in edgeRtranscripts-6h : 62878
transcripts in edgeRtranscripts-IFN : 60486
transcripts in edgeRtranscripts-IL4 : 63540

# wt vs ko alt splice
reczko@max:/data/images/proton/DKlab/orsalia/stringtie$ awk -f mergeEdgerDEresultsAndGetAltSpliced1.awk all_pairs_edger.txt > foo
ref transcripts 74118
transcripts in edgeRtranscripts-0h : 65358
transcripts in edgeRtranscripts-2h : 63764
transcripts in edgeRtranscripts-6h : 62878
transcripts in edgeRtranscripts-IFN : 60486
transcripts in edgeRtranscripts-IL4 : 63540

# pairs alt splice
for i in edgeRtr*.withGID.csv
 do
 echo $i
 awk -f getEdgerDEresultsAndGetAltSpliced1.awk $i > $i".alt_spliced.csv"
 wc $i".alt_spliced.csv"
 done

for i in edgeRtr*.withGID.csv
 do
 echo $i
 awk -f getEdgerDEresultsAndGetAltSplicedBest1.awk $i > $i".alt_spliced_best.csv"
 wc $i".alt_spliced_best.csv"
 done

edgeRtranscripts-0h.withGID.csv
ref transcripts 
  141  1967 31334 edgeRtranscripts-0h.withGID.csv.alt_spliced_best.csv
edgeRtranscripts-2h.withGID.csv
ref transcripts 
  161  2247 35961 edgeRtranscripts-2h.withGID.csv.alt_spliced_best.csv
edgeRtranscripts-6h.withGID.csv
ref transcripts 
   93  1295 20486 edgeRtranscripts-6h.withGID.csv.alt_spliced_best.csv
edgeRtranscripts-IFN.withGID.csv
ref transcripts 
  108  1505 23872 edgeRtranscripts-IFN.withGID.csv.alt_spliced_best.csv
edgeRtranscripts-IL4.withGID.csv
ref transcripts 
  102  1421 22645 edgeRtranscripts-IL4.withGID.csv.alt_spliced_best.csv
edgeRtranscripts_wtVSko.withGID.csv
ref transcripts 
   65   903 14167 edgeRtranscripts_wtVSko.withGID.csv.alt_spliced_best.csv


# Do a volcano plot
pdf("genes.counts.matrix.FAGsA_vs_FAGsB.edgeR.DE_results.MA_n_Volcano2.pdf", w=11, h=8.5)
plot(tTags$table$logFC, -log(tTags$table$PValue,10), main="Volcano plot of raw p-values", 
  pch=20, xlab="log2 ( Group B / Group A )", ylab="-log10 (raw p-value)")
#plot(tTags$table$logFC, -log(tTags$table$FDR,10), main="Volcano plot of FDR-adjusted p-values", 
#  pch=20, xlab="log2 ( Group B / Group A )", ylab="-log10 (FDR p-value)")
#plot(tTags$table$logFC, -log(tTags$table$PValue,10), main="Volcano plot showing genes with a raw p-value < 1e-05", 
#  pch=20, xlab="log2 ( Group B / Group A )", ylab="-log10 (raw p-value)", ylim=c(0,5))
#plot(tTags$table$logFC, -log(tTags$table$FDR,10), main="Volcano plot showing genes with a FDR p-value < 1e-05", 
#  pch=20, xlab="log2 ( Group B / Group A )", ylab="-log10 (FDR p-value)", ylim=c(0,5))
dev.off()
library(rtracklayer)
gff <- import.gff("Brapa_reference/Brapa_gene_v1.5.gff")
gff #take a look

#create a column "gene_id" that contains the gene name for every entry
gff$gene_id <- ifelse(is.na(gff$ID),gff$Parent,gff$ID)

export(gff,"Brapa_reference/Brapa_gene_v1.5.gtf",format="gtf")

#@
https://www.biostars.org/p/218136/
1
Hi,

you can get a table with all the information in your ballgown object (bg) , doing:

full_table <- texpr(bg , 'all')

The table will include: transcript id, chromosome, strand, start positon, end position, number of exons, length, gene id, gene name, and the FPKM values for each sample.
Hi!

I want to contribute the answer that AlyssaFrazee gave on Ballgown's github to retrieve GenesName (it's not intuitive, but rather a good workaround) Add geneNames to result of stattest on feature gene:

indices <- match(results_genes$id, texpr(bg, 'all')$gene_id)
gene_names_for_result <- texpr(bg, 'all')$gene_name[indices]
results_genes <- data.frame(geneNames=gene_names_for_result, results_genes)

#@
png("EdgeR Isoform expression in Ncoa6_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.9456" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("EdgeR Isoform expression in Bcl2l14_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("ENSMUSG00000030200" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("EdgeR Isoform expression in Fcgr2b_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.799" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("EdgeR Isoform expression in Mcart6_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("ENSMUSG00000044348" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("EdgeR Isoform expression in Ppp6r1_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.13695",ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()
png("EdgeR Isoform expression in Tmem176a_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.13043" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()



Fam49b
png("EdgeR Isoform expression in Fam49b_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.5375" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()

png("EdgeR top1-0h-Isoform expression in Zfp143_2.png",width=4000,height=2000,pointsize = 22)
plotTranscripts(ballgown::geneIDs(bg_chrX)[match("MSTRG.14498" ,ballgown::geneIDs(bg_chrX))], bg_chrX,		samples=c("wt01" ,   "wt02" ,   "wt21" ,   "wt22" ,   "wt61" ,   "wt62" ,   "wtIFN1" , "wtIFN2" , "wt01IL4" ,"wt02IL4", "Ko01" ,   "Ko02" ,   "KO21" ,   "KO22" ,   "Ko61" ,   "Ko62" ,   "KOIFN1" , "KOIFN2" , "Ko1IL4" , "Ko2IL4" ))
dev.off()