StringTie v1.3.3b usage: stringtie [-G ] [-l ] [-o ] [-p ] [-v] [-a ] [-m ] [-j ] [-f ] [-C ] [-c ] [-g ] [-u] [-e] [-x ] [-A ] [-h] {-B | -b } 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 -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 reference annotation to include in the merging (GTF/GFF3) -o output file name for the merged transcripts GTF (default: stdout) -m minimum input transcript length to include in the merge (default: 50) -c minimum input transcript coverage to include in the merge (default: 0) -F minimum input transcript FPKM to include in the merge (default: 1.0) -T minimum input transcript TPM to include in the merge (default: 1.0) -f minimum isoform fraction (default: 0.01) -g 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 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.gtf -A wt01.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A wt02.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A wt21.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A wt22.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A wt61.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A wt62.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A wtIFN1.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A wtIFN2.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A wt01IL4.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A wt02IL4.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A Ko01.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A Ko02.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A KO21.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A KO22.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A Ko61.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A Ko62.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A KOIFN1.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A KOIFN2.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A Ko1IL4.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.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.gtf -A Ko2IL4.tab -B -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf &> Ko2IL4.log /data/results/tools/rnaseq/stringtie/stringtie-1.3.3b.Linux_x86_64/stringtie --merge -o merged.gtf -G /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf Ko01.gtf Ko1IL4.gtf KO22.gtf Ko61.gtf KOIFN1.gtf wt01.gtf wt02.gtf wt21.gtf wt61.gtf wtIFN1.gtf Ko02.gtf KO21.gtf Ko2IL4.gtf Ko62.gtf KOIFN2.gtf wt01IL4.gtf wt02IL4.gtf wt22.gtf wt62.gtf wtIFN2.gtf /data/results/tools/rnaseq/stringtie/gffcompare-0.10.1.Linux_x86_64/gffcompare -r /data/results/reference/mmu/Mus_musculus.NCBIM37.64-toMM9.gtf -o comare1 merged.gtf /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.gtf -G merged.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 -e -B -o wt02/wt02.merged.gtf -G merged.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 -e -B -o wt21/wt21.merged.gtf -G merged.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 -e -B -o wt22/wt22.merged.gtf -G merged.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 -e -B -o wt61/wt61.merged.gtf -G merged.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 -e -B -o wt62/wt62.merged.gtf -G merged.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 -e -B -o wtIFN1/wtIFN1.merged.gtf -G merged.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 -e -B -o wtIFN2/wtIFN2.merged.gtf -G merged.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 -e -B -o wt01IL4/wt01IL4.merged.gtf -G merged.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 -e -B -o wt02IL4/wt02IL4.merged.gtf -G merged.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 -e -B -o Ko01/Ko01.merged.gtf -G merged.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 -e -B -o Ko02/Ko02.merged.gtf -G merged.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 -e -B -o KO21/KO21.merged.gtf -G merged.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 -e -B -o KO22/KO22.merged.gtf -G merged.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 -e -B -o Ko61/Ko61.merged.gtf -G merged.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 -e -B -o Ko62/Ko62.merged.gtf -G merged.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 -e -B -o KOIFN1/KOIFN1.merged.gtf -G merged.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 -e -B -o KOIFN2/KOIFN2.merged.gtf -G merged.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 -e -B -o Ko1IL4/Ko1IL4.merged.gtf -G merged.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 -e -B -o Ko2IL4/Ko2IL4.merged.gtf -G merged.gtf &> Ko2IL4.log mkdir ballgown mv Ko01 ballgown mv KO21 ballgown mv Ko61 ballgown mv KOIFN2 ballgown mv wt02IL4 ballgown mv wt22 ballgown mv wtIFN1 ballgown mv Ko02 ballgown mv KO22 ballgown mv Ko62 ballgown mv wt01IL4 ballgown mv wt02 ballgown mv wt61 ballgown mv wtIFN2 ballgown mv Ko1IL4 ballgown mv Ko2IL4 ballgown mv KOIFN1 ballgown mv wt01 ballgown mv wt21 ballgown mv wt62 ballgown /data/results/tools/rnaseq/stringtie/stringtie-master/prepDE.py 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("./ballgown/Ko01","./ballgown/Ko02","./ballgown/KO21","./ballgown/KO22","./ballgown/Ko61","./ballgown/Ko62","./ballgown/KOIFN1","./ballgown/KOIFN2","./ballgown/Ko1IL4","./ballgown/Ko2IL4","./ballgown/wt01","./ballgown/wt02","./ballgown/wt21","./ballgown/wt22","./ballgown/wt61","./ballgown/wt62","./ballgown/wtIFN1","./ballgown/wtIFN2","./ballgown/wt01IL4","./ballgown/wt02IL4" ) #@ save all: 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_results.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.gtf 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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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()