From 3ec78aa9af79921c44d62b65f88865a4b65880be Mon Sep 17 00:00:00 2001 From: Bo Li Date: Sun, 17 Apr 2011 23:13:54 -0500 Subject: [PATCH] RSEM v1.1.8 --- README.md | 20 ++++++++++++------- RefSeq.h | 4 ++-- Refs.h | 2 +- SamParser.h | 14 ++++++------- makefile | 4 ++-- parseIt.cpp | 22 +++++++++++++------- rsem-calculate-expression | 2 +- rsem-plot-model | 42 ++++++++++++++++++++++----------------- 8 files changed, 65 insertions(+), 45 deletions(-) diff --git a/README.md b/README.md index 5acd593..2bb6015 100644 --- a/README.md +++ b/README.md @@ -27,7 +27,7 @@ the EM algorithm, single-end and paired-end read data, quality scores, variable-length reads and RSPD estimation. It can also generate genomic-coordinate BAM files and UCSC wiggle files for visualization. In addition, it provides posterior mean and 95% credibility interval -estimates for expression levels. +estimates for expression levels. ## Compilation & Installation @@ -100,6 +100,12 @@ and provide the SAM or BAM file as an argument. When using an alternative aligner, you may also want to provide the --no-bowtie option to rsem-prepare-reference so that the Bowtie indices are not built. +However, please note that RSEM does ** not ** support gapped +alignments. So make sure that your aligner does not produce alignments +with intersions/deletions. Also, please make sure that you use +'reference_name.idx.fa' , which is generated by RSEM, to build your +aligner's indices. + ### III. Visualization RSEM contains a version of samtools in the 'sam' subdirectory. When @@ -133,16 +139,16 @@ RSEM provides an R script, 'rsem-plot-model', for visulazing the model learned. Usage: - rsem-plot-model modelF outF + rsem-plot-model sample_name outF -modelF: the sample_name.model file generated by RSEM +sample_name: the name of the sample analyzed outF: the file name for plots generated from the model. It is a pdf file The plots generated depends on read type and user configuration. It may include fragment length distribution, mate length distribution, read start position distribution (RSPD), quality score vs observed quality given a reference base, position vs percentage of sequencing -error given a reference base. +error given a reference base and histogram of read alignments. fragment length distribution and mate length distribution: x-axis is fragment/mate length, y axis is the probability of generating a fragment/mate with the associated length @@ -151,6 +157,8 @@ RSPD: Read Start Position Distribution. x-axis is bin number, y-axis is the prob Quality score vs. observed quality given a reference base: x-axis is Phred quality scores associated with data, y-axis is the "observed quality", Phred quality scores learned by RSEM from the data. Q = -10log_10(P), where Q is Phred quality score and P is the probability of sequencing error for a particular base Position vs. percentage sequencing error given a reference base: x-axis is position and y-axis is percentage sequencing error + +Histogram of read alignments: x-axis is the number of alignments a read has and y-axis is the number of such reads ## Example @@ -199,9 +207,7 @@ output_name.sim.isoforms.results, output_name.sim.genes.results : Results estima ## Acknowledgements -RSEM uses randomc.h and mersenne.cpp from - for random number generation. RSEM -also uses the [Boost C++](http://www.boost.org) and +RSEM uses the [Boost C++](http://www.boost.org) and [samtools](http://samtools.sourceforge.net) libraries. ## License diff --git a/RefSeq.h b/RefSeq.h index 1b2e58e..0e7f2ca 100644 --- a/RefSeq.h +++ b/RefSeq.h @@ -105,7 +105,7 @@ class RefSeq { std::vector fmasks; // record masks for forward strand, each position occupies 1 bit }; -//internal read; option 0 : read all 1 : do not read seqence and name +//internal read; option 0 : read all 1 : do not read seqences bool RefSeq::read(std::ifstream& fin, int option) { std::string line; @@ -122,7 +122,7 @@ bool RefSeq::read(std::ifstream& fin, int option) { getline(fin, line); assert(option == 0 || option == 1); - if (option == 1) { name = seq = ""; } + if (option == 1) { seq = ""; } return true; } diff --git a/Refs.h b/Refs.h index 41fcf06..f411622 100644 --- a/Refs.h +++ b/Refs.h @@ -107,7 +107,7 @@ void Refs::makeRefs(char *inpF, RefSeqPolicy& policy, PolyARules& rules) { } //inpF in fasta format, with sequence all in one line together -//option 0 read all, 1 do not read sequences and names +//option 0 read all, 1 do not read sequences void Refs::loadRefs(char *inpF, int option) { std::ifstream fin; RefSeq seq; diff --git a/SamParser.h b/SamParser.h index 335984e..425593a 100644 --- a/SamParser.h +++ b/SamParser.h @@ -13,8 +13,8 @@ #include "utils.h" -#include "Transcript.h" -#include "Transcripts.h" +#include "RefSeq.h" +#include "Refs.h" #include "SingleRead.h" #include "SingleReadQ.h" @@ -25,7 +25,7 @@ class SamParser { public: - SamParser(char, const char*, Transcripts&, const char* = 0); + SamParser(char, const char*, Refs&, const char* = 0); ~SamParser(); /** @@ -79,7 +79,7 @@ private: char SamParser::rtTag[STRLEN] = ""; // default : no tag, thus no Type 2 reads // aux, if not 0, points to the file name of fn_list -SamParser::SamParser(char inpType, const char* inpF, Transcripts& transcripts, const char* aux) { +SamParser::SamParser(char inpType, const char* inpF, Refs& refs, const char* aux) { switch(inpType) { case 'b': sam_in = samopen(inpF, "rb", aux); break; case 's': sam_in = samopen(inpF, "r", aux); break; @@ -91,14 +91,14 @@ SamParser::SamParser(char inpType, const char* inpF, Transcripts& transcripts, c if (header == 0) { fprintf(stderr, "Fail to parse sam header!\n"); exit(-1); } // Check if the reference used for aligner is the transcript set RSEM generated - if (transcripts.getM() != header->n_targets) { + if (refs.getM() != header->n_targets) { fprintf(stderr, "Number of transcripts does not match! Please align reads against the transcript set and use RSEM generated reference for your aligner!\n"); exit(-1); } for (int i = 0; i < header->n_targets; i++) { - const Transcript& transcript = transcripts.getTranscriptAt(i + 1); + const RefSeq& refseq = refs.getRef(i + 1); // If update int to long, chance the (int) conversion - if (transcript.getTranscriptID().compare(header->target_name[i]) != 0 || transcript.getLength() != (int)header->target_len[i]) { + if (refseq.getName().compare(header->target_name[i]) != 0 || refseq.getTotLen() != (int)header->target_len[i]) { fprintf(stderr, "Transcript information does not match! Please align reads against the transcript set and use RSEM generated reference for your aligner!\n"); exit(-1); } diff --git a/makefile b/makefile index df8d2a1..12fc9d6 100644 --- a/makefile +++ b/makefile @@ -48,13 +48,13 @@ PairedEndHit.h : SingleHit.h HitContainer.h : GroupInfo.h -SamParser.h : sam/sam.h sam/bam.h utils.h SingleRead.h SingleReadQ.h PairedEndRead.h PairedEndReadQ.h SingleHit.h PairedEndHit.h Transcript.h Transcripts.h +SamParser.h : sam/sam.h sam/bam.h utils.h SingleRead.h SingleReadQ.h PairedEndRead.h PairedEndReadQ.h SingleHit.h PairedEndHit.h RefSeq.h Refs.h rsem-parse-alignments : parseIt.o sam/libbam.a $(CC) -o rsem-parse-alignments parseIt.o sam/libbam.a -lz -parseIt.o : utils.h GroupInfo.h Read.h SingleRead.h SingleReadQ.h PairedEndRead.h PairedEndReadQ.h SingleHit.h PairedEndHit.h HitContainer.h SamParser.h Transcript.h Transcripts.h sam/sam.h sam/bam.h parseIt.cpp +parseIt.o : utils.h GroupInfo.h Read.h SingleRead.h SingleReadQ.h PairedEndRead.h PairedEndReadQ.h SingleHit.h PairedEndHit.h HitContainer.h SamParser.h RefSeq.h Refs.h sam/sam.h sam/bam.h parseIt.cpp $(CC) $(COFLAGS) parseIt.cpp diff --git a/parseIt.cpp b/parseIt.cpp index 55dd954..e05585d 100644 --- a/parseIt.cpp +++ b/parseIt.cpp @@ -14,8 +14,8 @@ #include "GroupInfo.h" -#include "Transcript.h" -#include "Transcripts.h" +#include "RefSeq.h" +#include "Refs.h" #include "SingleRead.h" #include "SingleReadQ.h" @@ -34,12 +34,12 @@ int N[3]; // note, N = N0 + N1 + N2 , but may not be equal to the total number o int nHits; // # of hits int nUnique, nMulti, nIsoMulti; char fn_list[STRLEN]; -char groupF[STRLEN], tiF[STRLEN]; +char refF[STRLEN], groupF[STRLEN]; char imdName[STRLEN]; char datF[STRLEN], cntF[STRLEN]; +Refs refs; GroupInfo gi; -Transcripts transcripts; SamParser *parser; ofstream hit_out; @@ -55,7 +55,7 @@ void init(const char* imdName, char alignFType, const char* alignF) { char* aux = 0; if (strcmp(fn_list, "")) aux = fn_list; - parser = new SamParser(alignFType, alignF, transcripts, aux); + parser = new SamParser(alignFType, alignF, refs, aux); memset(cat, 0, sizeof(cat)); memset(readOutFs, 0, sizeof(readOutFs)); @@ -96,6 +96,9 @@ void parseIt(SamParser *parser) { if (record_val >= 0) { record_read.write(n_os, cat[record_val]); ++N[record_val]; + } + // flush out previous read's hits if the read is alignable reads + if (record_val == 1) { hits.updateRI(); nHits += hits.getNHits(); nMulti += hits.calcNumGeneMultiReads(gi); @@ -127,6 +130,9 @@ void parseIt(SamParser *parser) { if (record_val >= 0) { record_read.write(n_os, cat[record_val]); ++N[record_val]; + } + + if (record_val == 1) { hits.updateRI(); nHits += hits.getNHits(); nMulti += hits.calcNumGeneMultiReads(gi); @@ -186,10 +192,10 @@ int main(int argc, char* argv[]) { verbose = !quiet; + sprintf(refF, "%s.seq", argv[1]); + refs.loadRefs(refF, 1); sprintf(groupF, "%s.grp", argv[1]); gi.load(groupF); - sprintf(tiF, "%s.ti", argv[1]); - transcripts.readFrom(tiF); sprintf(imdName, "%s.temp/%s", argv[2], argv[3]); sprintf(datF, "%s.dat", imdName); @@ -220,9 +226,11 @@ int main(int argc, char* argv[]) { fout<first<<'\t'<second<