annotate test-data/reference.R @ 5:423ad61697c4

Bugfix 1: [triplet] Lengths (frameCounting) if given should be a range (not zero). readingFrame function fails with subscript out of bounds. Bugfix 2: [triplet] Check if transcript name in SAM matches the name in FASTA. Produce an error message if it's not. This fixes the problem where an empty plot is produced (no bars). [ribosome_profile] - A proper error message is now produced if an invalid transcript name is provided. Updated test data

# Script used to prepare the input data for Triplet periodicity step. It is done this way
# as SAM files were not available for this dataset at the time this release was made
library(riboSeqR)
fastaCDS <- findCDS(fastaFile='rsem_chlamy236_deNovo.transcripts.fa', startCodon=c("ATG"), stopCodon=c("TAG", "TAA", "TGA"))

# RiboSeq file 1-4 are just links to chlamy236_plus_deNovo_plusOnly_Index[17,3,5,7]
ribofiles <- paste(c("RiboSeq file 1", "RiboSeq file 2", "RiboSeq file 3", "RiboSeq file 4"))
rnafiles <- paste(c("RNASeq file 1", "RNASeq file 2", "RNASeq file 3", "RNASeq file 4"))
riboDat <- readRibodata(ribofiles, rnafiles, replicates=c("WT", "WT", "M", "M"))
save('riboDat', file='Prepare riboSeqR input (R data file)', compress=FALSE)

fCs <- frameCounting(riboDat, fastaCDS)
fS <- readingFrame(rC = fCs)

#pdf(file="/tmp/Periodicity-plot.pdf")
#plotFS(fS, legend.text = c('Frame 0', 'Frame 1', 'Frame 2'))
#dev.off()

ffCs <- filterHits(fCs, lengths=c(27,28), frames=list(1,0), hitMean = 50, unqhitMean = 10, fS=fS)
save('fastaCDS', 'riboDat', 'fCs', 'fS', 'ffCs', file='Metagene analysis (R data file)', compress=FALSE)

#pdf(file='/tmp/Metagene-analysis-plot.pdf')
#plotCDS(coordinates=ffCs@CDS, riboDat=riboDat, lengths=27)
#plotCDS(coordinates=ffCs@CDS, riboDat=riboDat, lengths=28)
#dev.off()
#
#pdf(file='/tmp/Ribosome-profile-plot.pdf')
#plotTranscript("CUFF.37930.1", coordinates=ffCs@CDS, riboData=riboDat, length=27, cap=200)
#dev.off()
#
#riboCounts <- sliceCounts(ffCs, lengths = c(27,28), frames = list(0,2))
#rnaCounts <- rnaCounts(riboDat, ffCs@CDS)

#save('fastaCDS', 'riboDat', 'fCs', 'fS', 'ffCs', 'riboCounts', 'rnaCounts', file='Robjects.rda')