From 6eec553ba4ab1cb1c6f48cd14d9e5c92c7d85798 Mon Sep 17 00:00:00 2001 From: Bo Li Date: Tue, 5 Jun 2012 21:02:46 -0500 Subject: [PATCH 01/16] Tested the codes --- BamConverter.h | 2 +- WHAT_IS_NEW | 8 ++++++++ rsem-calculate-expression | 17 ++++++++--------- 3 files changed, 17 insertions(+), 10 deletions(-) diff --git a/BamConverter.h b/BamConverter.h index 13b8557..eaf2f66 100644 --- a/BamConverter.h +++ b/BamConverter.h @@ -101,7 +101,7 @@ void BamConverter::process() { if (!notgood) { - assert((b->core.tid == b2->core.tid) && (b->core.flag & 0x0040) && (b2->core.flag & 0x0080)); // for collapsing + if (isPaired) assert((b->core.tid == b2->core.tid) && (b->core.flag & 0x0040) && (b2->core.flag & 0x0080)); // for collapsing const Transcript& transcript = transcripts.getTranscriptViaEid(b->core.tid + 1); diff --git a/WHAT_IS_NEW b/WHAT_IS_NEW index ddcfd95..3c7f654 100644 --- a/WHAT_IS_NEW +++ b/WHAT_IS_NEW @@ -1,3 +1,11 @@ +RSEM v1.1.20 + +- Added an option to set the temporary folder name +- Removed sample_name.sam.gz. Instead, RSEM uses samtools to convert bowtie outputed SAM file into a BAM file under the temporary folder +- RSEM generated BAM files now contains all alignment lines produced by bowtie or user-specified aligners, including unalignable reads. Please note that for paired-end reads, if one mate has alignments but the other does not, RSEM will mark the alignable mate as "unmappable" (flag bit 0x4) and append an optional field "Z0:A:!" + +-------------------------------------------------------------------------------------------- + RSEM v1.1.19 - Allowed > 2^31 hits diff --git a/rsem-calculate-expression b/rsem-calculate-expression index 457f177..9e6ef2e 100755 --- a/rsem-calculate-expression +++ b/rsem-calculate-expression @@ -310,7 +310,7 @@ if ($calcCI) { &collectResults("$imdName.iso_res", "$sampleName.isoforms.results"); # isoform level &collectResults("$imdName.gene_res", "$sampleName.genes.results"); # gene level - $command = $dir."rsem-calculate-credibility-intervals $refName $sampleName $sampleToken $CONFIDENCE $NCV $NSPC $NMB"; + $command = $dir."rsem-calculate-credibility-intervals $refName $imdName $statName $CONFIDENCE $NCV $NSPC $NMB"; $command .= " -p $nThreads"; if ($quiet) { $command .= " -q"; } &runCommand($command); @@ -482,7 +482,7 @@ Generate a BAM file, 'sample_name.genome.bam', with alignments mapped to genomic =item B<--sampling-for-bam> -When RSEM generates a BAM file, instead of outputing all alignments a read has with their posterior probabilities, one alignment is sampled and outputed according to the posterior probabilities. If the sampling result is that the read comes from the "noise" transcript, nothing is outputed. (Default: off) +When RSEM generates a BAM file, instead of outputing all alignments a read has with their posterior probabilities, one alignment is sampled according to the posterior probabilities. The sampling procedure includes the alignment to the "noise" transcript, which does not appear in the BAM file. Only the sampled alignment has a weight of 1. All other alignments have weight 0. If the "noise" transcript is sampled, all alignments appeared in the BAM file should have weight 0. (Default: off) =item B<--calc-ci> @@ -646,7 +646,12 @@ is set to min(100, floor(-10 * log10(1.0 - w) + 0.5)), where w is the posterior probability of that alignment being the true mapping of a read. In addition, RSEM pads a new tag ZW:f:value, where value is a single precision floating number representing the posterior -probability. +probability. Because this file contains all alignment lines produced +by bowtie or user-specified aligners, it can also be used as a +replacement of the aligner generated BAM/SAM file. For paired-end +reads, if one mate has alignments but the other does not, this file +marks the alignable mate as "unmappable" (flag bit 0x4) and appends an +optional field "Z0:A:!". 'sample_name.transcript.sorted.bam' and 'sample_name.transcript.sorted.bam.bai' are the sorted BAM file and @@ -671,12 +676,6 @@ indicating the strand of the transcript it aligns to. 'sample_name.genome.sorted.bam' and 'sample_name.genome.sorted.bam.bai' are the sorted BAM file and indices generated by samtools (included in RSEM package). -=item B - -Only generated when the input files are raw reads instead of SAM/BAM format files - -It is the gzipped SAM output produced by bowtie aligner. - =item B Only generated when --time is specified. -- 2.39.2 From 1b2999c4407ef8419fb89b66b843b7141caff4da Mon Sep 17 00:00:00 2001 From: Bo Li Date: Mon, 2 Jul 2012 12:18:03 -0500 Subject: [PATCH 02/16] rsem v1.1.21 --- BamWriter.h | 4 +++- README.md | 9 +++++---- WHAT_IS_NEW | 10 +++++++++- bam2readdepth.cpp | 19 +++++++++++++------ bam2wig.cpp | 16 ++++++++++++---- makefile | 4 ++-- my_assert.h | 12 ++++++------ rsem-plot-transcript-wiggles | 8 ++++++-- utils.h | 6 +++--- wiggle.cpp | 29 +++++++++++++++++++---------- wiggle.h | 7 +++++++ 11 files changed, 85 insertions(+), 39 deletions(-) diff --git a/BamWriter.h b/BamWriter.h index 782f5bd..f664710 100644 --- a/BamWriter.h +++ b/BamWriter.h @@ -144,6 +144,8 @@ void BamWriter::work(HitWrapper wrapper) { b->core.mpos = b2->core.pos; b2->core.mpos = b->core.pos; } + + /* else { // if only one mate can be aligned, mask it as unaligned and put an additional tag Z0:A:! char exclamation = '!'; @@ -156,7 +158,7 @@ void BamWriter::work(HitWrapper wrapper) { bam_aux_append(b2, "Z0", 'A', bam_aux_type2size('A'), (uint8_t*)&exclamation); } } - + */ samwrite(out, b); samwrite(out, b2); diff --git a/README.md b/README.md index ed5acc0..8d3c15f 100644 --- a/README.md +++ b/README.md @@ -158,11 +158,12 @@ plot, run the 'rsem-bam2wig' program on the Usage: - rsem-bam2wig sorted_bam_input wig_output wiggle_name + rsem-bam2wig sorted_bam_input wig_output wiggle_name [--no-fractional-weight] -sorted_bam_input: sorted bam file -wig_output: output file name, e.g. output.wig -wiggle_name: the name the user wants to use for this wiggle plot +sorted_bam_input : Input BAM format file, must be sorted +wig_output : Output wiggle file's name, e.g. output.wig +wiggle_name : the name of this wiggle plot +--no-fractional-weight : If this is set, RSEM will not look for "ZW" tag and each alignment appeared in the BAM file has weight 1. Set this if your BAM file is not generated by RSEM. Please note that this option must be at the end of the command line. #### b) Loading a BAM and/or Wiggle file into the UCSC Genome Browser or Integrative Genomics Viewer(IGV) diff --git a/WHAT_IS_NEW b/WHAT_IS_NEW index 3c7f654..244066c 100644 --- a/WHAT_IS_NEW +++ b/WHAT_IS_NEW @@ -1,7 +1,15 @@ +RSEM v1.1.21 + +- Removed optional field "Z0:A:!" in the BAM outputs +- Added --no-fractional-weight option to rsem-bam2wig, if the BAM file is not generated by RSEM, this option is recommended to be set +- Fixed a bug for generating transcript level wiggle files using 'rsem-plot-transcript-wiggles' + +-------------------------------------------------------------------------------------------- + RSEM v1.1.20 - Added an option to set the temporary folder name -- Removed sample_name.sam.gz. Instead, RSEM uses samtools to convert bowtie outputed SAM file into a BAM file under the temporary folder +- Removed sample_name.sam.gz. Instead, RSEM uses samtools to convert bowtie outputted SAM file into a BAM file under the temporary folder - RSEM generated BAM files now contains all alignment lines produced by bowtie or user-specified aligners, including unalignable reads. Please note that for paired-end reads, if one mate has alignments but the other does not, RSEM will mark the alignable mate as "unmappable" (flag bit 0x4) and append an optional field "Z0:A:!" -------------------------------------------------------------------------------------------- diff --git a/bam2readdepth.cpp b/bam2readdepth.cpp index 87fb178..d252e2e 100644 --- a/bam2readdepth.cpp +++ b/bam2readdepth.cpp @@ -1,20 +1,27 @@ #include #include #include -#include +#include +#include "my_assert.h" #include "wiggle.h" using namespace std; int main(int argc, char* argv[]) { - if (argc != 2) { - printf("Usage: rsem-bam2readdepth sorted_bam_input\n"); + if (argc != 3) { + printf("Usage: rsem-bam2readdepth sorted_bam_input readdepth_output\n"); exit(-1); } - ReadDepthWriter depth_writer(std::cout); - build_wiggles(argv[1], depth_writer); + ofstream fout(argv[2]); + general_assert(fout.is_open(), "Cannot write to " + cstrtos(argv[2]) + "!"); - return 0; + ReadDepthWriter depth_writer(fout); + + build_wiggles(argv[1], depth_writer); + + fout.close(); + + return 0; } diff --git a/bam2wig.cpp b/bam2wig.cpp index 24a53c6..d205e8f 100644 --- a/bam2wig.cpp +++ b/bam2wig.cpp @@ -6,12 +6,20 @@ using namespace std; +void printUsage() { + printf("Usage: rsem-bam2wig sorted_bam_input wig_output wiggle_name [--no-fractional-weight]\n"); + printf("sorted_bam_input\t: Input BAM format file, must be sorted\n"); + printf("wig_output\t\t: Output wiggle file's name, e.g. output.wig\n"); + printf("wiggle_name\t\t: the name of this wiggle plot\n"); + printf("--no-fractional-weight\t: If this is set, RSEM will not look for \"ZW\" tag and each alignment appeared in the BAM file has weight 1. Set this if your BAM file is not generated by RSEM. Please note that this option must be at the end of the command line.\n"); + exit(-1); +} + int main(int argc, char* argv[]) { - if (argc != 4) { - printf("Usage: rsem-bam2wig sorted_bam_input wig_output wiggle_name\n"); - exit(-1); - } + if (argc < 4 || argc > 5) { printf("Number of arguments is not correct!\n"); printUsage(); } + if (argc == 5 && strcmp(argv[4], "--no-fractional-weight")) { printf("Cannot recognize option %s!\n", argv[4]); printUsage(); } + no_fractional_weight = (argc == 5 && !strcmp(argv[4], "--no-fractional-weight")); UCSCWiggleTrackWriter track_writer(argv[2], argv[3]); build_wiggles(argv[1], track_writer); diff --git a/makefile b/makefile index c9f37a8..e1a47f2 100644 --- a/makefile +++ b/makefile @@ -96,10 +96,10 @@ BamConverter.h : utils.h my_assert.h sam/sam.h sam/bam.h sam_rsem_aux.h sam_rsem rsem-tbam2gbam : utils.h Transcripts.h Transcript.h bc_aux.h BamConverter.h sam/sam.h sam/bam.h sam/libbam.a sam_rsem_aux.h sam_rsem_cvt.h tbam2gbam.cpp sam/libbam.a $(CC) -O3 -Wall tbam2gbam.cpp sam/libbam.a -lz -o $@ -rsem-bam2wig : wiggle.h wiggle.o sam/libbam.a bam2wig.cpp +rsem-bam2wig : utils.h my_assert.h wiggle.h wiggle.o sam/libbam.a bam2wig.cpp $(CC) -O3 -Wall bam2wig.cpp wiggle.o sam/libbam.a -lz -o $@ -rsem-bam2readdepth : wiggle.h wiggle.o sam/libbam.a bam2readdepth.cpp +rsem-bam2readdepth : utils.h my_assert.h wiggle.h wiggle.o sam/libbam.a bam2readdepth.cpp $(CC) -O3 -Wall bam2readdepth.cpp wiggle.o sam/libbam.a -lz -o $@ wiggle.o: sam/bam.h sam/sam.h wiggle.cpp wiggle.h diff --git a/my_assert.h b/my_assert.h index 02844e1..d63c3a2 100644 --- a/my_assert.h +++ b/my_assert.h @@ -9,29 +9,29 @@ #include #include -std::string itos(int i) { +inline std::string itos(int i) { std::ostringstream strout; strout< "Invalid number of arguments!", -exitval => 2, -verbose => 2) my ($fn, $dir, $suf) = fileparse($0); my $command = ""; +unless (-e "$ARGV[0].transcript.sorted.bam") { + $command = $dir."sam/samtools sort $ARGV[0].transcript.bam $ARGV[0].transcript.sorted"; + &runCommand($command); +} unless (-e "$ARGV[0].transcript.readdepth") { - $command = $dir."rsem-bam2readdepth $ARGV[0].transcript.sorted.bam > $ARGV[0].transcript.readdepth"; + $command = $dir."rsem-bam2readdepth $ARGV[0].transcript.sorted.bam $ARGV[0].transcript.readdepth"; &runCommand($command); } @@ -34,7 +38,7 @@ if ($show_unique) { &runCommand($command); } unless (-e "$ARGV[0].uniq.transcript.readdepth") { - $command = $dir."rsem-bam2readdepth $ARGV[0].uniq.transcript.sorted.bam > $ARGV[0].uniq.transcript.readdepth"; + $command = $dir."rsem-bam2readdepth $ARGV[0].uniq.transcript.sorted.bam $ARGV[0].uniq.transcript.readdepth"; &runCommand($command); } } diff --git a/utils.h b/utils.h index b15f154..137440e 100644 --- a/utils.h +++ b/utils.h @@ -23,7 +23,7 @@ const int RANGE = 201; const int OLEN = 25; // overlap length, number of bases must not be in poly(A) tails const int NBITS = 32; // use unsigned int, 32 bits per variable -bool verbose = true; // show detail intermediate outputs +static bool verbose = true; // show detail intermediate outputs inline bool isZero(double a) { return fabs(a) < 1e-8; } inline bool isLongZero(double a) { return fabs(a) < 1e-30; } @@ -124,7 +124,7 @@ inline std::string cleanStr(const std::string& str) { return (fr <= to ? str.substr(fr, to - fr + 1) : ""); } -void genReadFileNames(const char* readFN, int tagType, int read_type, int& s, char readFs[][STRLEN]){ +inline void genReadFileNames(const char* readFN, int tagType, int read_type, int& s, char readFs[][STRLEN]){ const char tags[3][STRLEN] = {"un", "alignable", "max"}; char suffix[STRLEN]; @@ -146,7 +146,7 @@ void genReadFileNames(const char* readFN, int tagType, int read_type, int& s, ch } } -void printTimeUsed(const time_t& a, const time_t& b, const char* program_name) { +inline void printTimeUsed(const time_t& a, const time_t& b, const char* program_name) { int hh = (b - a) / 3600; int mm = (b - a) % 3600 / 60; int ss = (b - a) % 60; diff --git a/wiggle.cpp b/wiggle.cpp index b0cb7d3..4e68b44 100644 --- a/wiggle.cpp +++ b/wiggle.cpp @@ -10,9 +10,18 @@ #include "utils.h" #include "wiggle.h" +bool no_fractional_weight = false; + void add_bam_record_to_wiggle(const bam1_t *b, Wiggle& wiggle) { - uint8_t *p_tag = bam_aux_get(b, "ZW"); - float w = (p_tag != NULL ? bam_aux2f(p_tag) : 1.0); + float w; + + if (no_fractional_weight) w = 1.0; + else { + uint8_t *p_tag = bam_aux_get(b, "ZW"); + if (p_tag == NULL) return; + w = bam_aux2f(p_tag); + } + int pos = b->core.pos; uint32_t *p = bam1_cigar(b); @@ -43,21 +52,21 @@ void build_wiggles(const std::string& bam_filename, bool *used = new bool[header->n_targets]; memset(used, 0, sizeof(bool) * header->n_targets); - int cur_tid = -1; //current tid; - HIT_INT_TYPE cnt = 0; - bam1_t *b = bam_init1(); - Wiggle wiggle; + int cur_tid = -1; //current tid; + HIT_INT_TYPE cnt = 0; + bam1_t *b = bam_init1(); + Wiggle wiggle; while (samread(bam_in, b) >= 0) { if (b->core.flag & 0x0004) continue; if (b->core.tid != cur_tid) { if (cur_tid >= 0) { used[cur_tid] = true; processor.process(wiggle); } cur_tid = b->core.tid; - wiggle.name = header->target_name[cur_tid]; - wiggle.length = header->target_len[cur_tid]; - wiggle.read_depth.assign(wiggle.length, 0.0); + wiggle.name = header->target_name[cur_tid]; + wiggle.length = header->target_len[cur_tid]; + wiggle.read_depth.assign(wiggle.length, 0.0); } - add_bam_record_to_wiggle(b, wiggle); + add_bam_record_to_wiggle(b, wiggle); ++cnt; if (cnt % 1000000 == 0) std::cout<< cnt<< std::endl; } diff --git a/wiggle.h b/wiggle.h index 1f37592..ecce0fc 100644 --- a/wiggle.h +++ b/wiggle.h @@ -1,8 +1,13 @@ +#ifndef WIGGLE_H_ +#define WIGGLE_H_ + #include #include #include #include +extern bool no_fractional_weight; // if no_frac_weight == true, each alignment counts as weight 1 + struct Wiggle { std::string name; std::vector read_depth; @@ -40,3 +45,5 @@ private: void build_wiggles(const std::string& bam_filename, WiggleProcessor& processor); + +#endif -- 2.39.2 From c678676b7da64266b4b98258a34d1b32372da6dd Mon Sep 17 00:00:00 2001 From: Bo Li Date: Sun, 15 Jul 2012 15:46:46 -0500 Subject: [PATCH 03/16] Fix a bug on convert-sam-for-rsem for the case that a paired-end read can align to the same position with different directions --- scanForPairedEndReads.cpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/scanForPairedEndReads.cpp b/scanForPairedEndReads.cpp index 751fa55..b220fa8 100644 --- a/scanForPairedEndReads.cpp +++ b/scanForPairedEndReads.cpp @@ -29,20 +29,28 @@ inline void add_to_appropriate_arr(bam1_t *b) { else arr_partial_unknown.push_back(bam_dup1(b)); } +char get_pattern_code(uint32_t flag) { + if (flag & 0x0040) return (flag & 0x0010 ? 1 : 0); + else return (flag & 0x0010 ? 0 : 1); +} + bool less_than(bam1_t *a, bam1_t *b) { int32_t ap1 = min(a->core.pos, a->core.mpos); int32_t ap2 = max(a->core.pos, a->core.mpos); int32_t bp1 = min(b->core.pos, b->core.mpos); int32_t bp2 = max(b->core.pos, b->core.mpos); + char apat = get_pattern_code(a->core.flag); // apt: a's pattern of strand and mate information + char bpat = get_pattern_code(b->core.flag); if (a->core.tid != b->core.tid) return a->core.tid < b->core.tid; if (ap1 != bp1) return ap1 < bp1; - return ap2 < bp2; + if (ap2 != bp2) return ap2 < bp2; + return apat < bpat; } int main(int argc, char* argv[]) { if (argc != 3) { - printf("Usage: rsem-scan-for-paired-end-reads input.sam output.bam\n"); + printf("UsaOAge: rsem-scan-for-paired-end-reads input.sam output.bam\n"); exit(-1); } -- 2.39.2 From 683863b75f8d8bef2461039a6911b0e9619cc113 Mon Sep 17 00:00:00 2001 From: Bo Li Date: Sat, 21 Jul 2012 21:11:07 -0500 Subject: [PATCH 04/16] Fixed a bug in perl scripts for printing error messages --- convert-sam-for-rsem | 8 ++++---- rsem-calculate-expression | 12 ++++++------ rsem-generate-ngvector | 8 ++++---- rsem-plot-transcript-wiggles | 1 + rsem-prepare-reference | 12 ++++++------ 5 files changed, 21 insertions(+), 20 deletions(-) diff --git a/convert-sam-for-rsem b/convert-sam-for-rsem index 13b9b13..decf861 100755 --- a/convert-sam-for-rsem +++ b/convert-sam-for-rsem @@ -82,16 +82,16 @@ $command = $dir."rsem-sam-validator $out_file"; # command, {err_msg} sub runCommand { - print STDERR $_[0]."\n"; + print $_[0]."\n"; my $status = system($_[0]); if ($status != 0) { my $errmsg; if (scalar(@_) > 1) { $errmsg = $_[1]; } - else { $errmsg = "\"$command\" failed! Plase check if you provide correct parameters/options for the script!"; } - print STDERR $errmsg."\n"; + else { $errmsg = "\"$command\" failed! Plase check if you provide correct parameters/options for the pipeline!"; } + print $errmsg."\n"; exit(-1); } - print STDERR "\n"; + print "\n"; } __END__ diff --git a/rsem-calculate-expression b/rsem-calculate-expression index 9e6ef2e..684e33c 100755 --- a/rsem-calculate-expression +++ b/rsem-calculate-expression @@ -345,12 +345,12 @@ if ($mTime) { sub runCommand { print $_[0]."\n"; my $status = system($_[0]); - if ($status != 0) { - my $errmsg; - if (scalar(@_) > 1) { $errmsg = $_[1]; } - else { $errmsg = "\"$command\" failed! Plase check if you provide correct parameters/options for the pipeline!"; } - print $errmsg."\n"; - exit(-1); + if ($status != 0) { + my $errmsg = ""; + if (scalar(@_) > 1) { $errmsg .= $_[1]."\n"; } + $errmsg .= "\"$_[0]\" failed! Plase check if you provide correct parameters/options for the pipeline!\n"; + print $errmsg; + exit(-1); } print "\n"; } diff --git a/rsem-generate-ngvector b/rsem-generate-ngvector index c4d26a8..cd184b6 100755 --- a/rsem-generate-ngvector +++ b/rsem-generate-ngvector @@ -28,10 +28,10 @@ sub runCommand { print $_[0]."\n"; my $status = system($_[0]); if ($status != 0) { - my $errmsg; - if (scalar(@_) > 1) { $errmsg = $_[1]; } - else { $errmsg = "\"$command\" failed! Plase check if you provide correct parameters/options for the pipeline!"; } - print $errmsg."\n"; + my $errmsg = ""; + if (scalar(@_) > 1) { $errmsg .= $_[1]."\n"; } + $errmsg .= "\"$_[0]\" failed! Plase check if you provide correct parameters/options for the pipeline!\n"; + print $errmsg; exit(-1); } print "\n"; diff --git a/rsem-plot-transcript-wiggles b/rsem-plot-transcript-wiggles index 25e8b35..37c7e27 100755 --- a/rsem-plot-transcript-wiggles +++ b/rsem-plot-transcript-wiggles @@ -46,6 +46,7 @@ if ($show_unique) { $command = $dir."rsem-gen-transcript-plots $ARGV[0] $ARGV[1] $gene_list $show_unique $ARGV[2]"; &runCommand($command); + # command, {err_msg} sub runCommand { print $_[0]."\n"; diff --git a/rsem-prepare-reference b/rsem-prepare-reference index 80a9549..56d6972 100755 --- a/rsem-prepare-reference +++ b/rsem-prepare-reference @@ -97,12 +97,12 @@ if (!$no_bowtie) { sub runCommand { print $_[0]."\n"; my $status = system($_[0]); - if ($status != 0) { - my $errmsg; - if (scalar(@_) > 1) { $errmsg = $_[1]; } - else { $errmsg = "\"$command\" failed! Plase check if you provide correct parameters/options for the pipeline!"; } - print $errmsg."\n"; - exit(-1); + if ($status != 0) { + my $errmsg = ""; + if (scalar(@_) > 1) { $errmsg .= $_[1]."\n"; } + $errmsg .= "\"$_[0]\" failed! Plase check if you provide correct parameters/options for the pipeline!\n"; + print $errmsg; + exit(-1); } print "\n"; } -- 2.39.2 From 7f503f1b9e22d7e1b44472add2263f263675b7c7 Mon Sep 17 00:00:00 2001 From: Bo Li Date: Sun, 22 Jul 2012 19:01:16 -0500 Subject: [PATCH 05/16] Fixed a bug for calling pthread_join --- EM.cpp | 7 +++---- Gibbs.cpp | 3 +-- calcCI.cpp | 5 ++--- 3 files changed, 6 insertions(+), 9 deletions(-) diff --git a/EM.cpp b/EM.cpp index 0d38728..684d639 100644 --- a/EM.cpp +++ b/EM.cpp @@ -423,7 +423,6 @@ void EM() { Params fparams[nThreads]; pthread_t threads[nThreads]; pthread_attr_t attr; - void *status; int rc; @@ -471,7 +470,7 @@ void EM() { } for (int i = 0; i < nThreads; i++) { - rc = pthread_join(threads[i], &status); + rc = pthread_join(threads[i], NULL); pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0) at ROUND " + itos(ROUND) + "!"); } @@ -521,7 +520,7 @@ void EM() { pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0) when generating files for Gibbs sampler!"); } for (int i = 0; i < nThreads; i++) { - rc = pthread_join(threads[i], &status); + rc = pthread_join(threads[i], NULL); pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0) when generating files for Gibbs sampler!"); } } @@ -581,7 +580,7 @@ void EM() { pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0) when calculating expected weights!"); } for (int i = 0; i < nThreads; i++) { - rc = pthread_join(threads[i], &status); + rc = pthread_join(threads[i], NULL); pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0) when calculating expected weights!"); } model.setNeedCalcConPrb(false); diff --git a/Gibbs.cpp b/Gibbs.cpp index 63d35a8..6bbfeb1 100644 --- a/Gibbs.cpp +++ b/Gibbs.cpp @@ -70,7 +70,6 @@ bool quiet; Params *paramsArray; pthread_t *threads; pthread_attr_t attr; -void *status; int rc; void load_data(char* reference_name, char* statName, char* imdName) { @@ -443,7 +442,7 @@ int main(int argc, char* argv[]) { pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0)!"); } for (int i = 0; i < nThreads; i++) { - rc = pthread_join(threads[i], &status); + rc = pthread_join(threads[i], NULL); pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0)!"); } release(); diff --git a/calcCI.cpp b/calcCI.cpp index 99ce148..c71d9c1 100644 --- a/calcCI.cpp +++ b/calcCI.cpp @@ -69,7 +69,6 @@ bool quiet; Params *paramsArray; pthread_t *threads; pthread_attr_t attr; -void *status; int rc; CIParams *ciParamsArray; @@ -219,7 +218,7 @@ void sample_theta_vectors_from_count_vectors() { pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0) in sample_theta_vectors_from_count_vectors!"); } for (int i = 0; i < num_threads; i++) { - rc = pthread_join(threads[i], &status); + rc = pthread_join(threads[i], NULL); pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0) in sample_theta_vectors_from_count_vectors!"); } @@ -359,7 +358,7 @@ void calculate_credibility_intervals(char* imdName) { pthread_assert(rc, "pthread_create", "Cannot create thread " + itos(i) + " (numbered from 0) in calculate_credibility_intervals!"); } for (int i = 0; i < num_threads; i++) { - rc = pthread_join(threads[i], &status); + rc = pthread_join(threads[i], NULL); pthread_assert(rc, "pthread_join", "Cannot join thread " + itos(i) + " (numbered from 0) in calculate_credibility_intervals!"); } -- 2.39.2 From 7777bfc71742d4e239d1297d2e3de058895f4ce1 Mon Sep 17 00:00:00 2001 From: Bo Li Date: Fri, 27 Jul 2012 11:46:25 -0500 Subject: [PATCH 06/16] Added a user-friendly error message for rsem-sam-validator --- samValidator.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/samValidator.cpp b/samValidator.cpp index 5b02ba4..6e60210 100644 --- a/samValidator.cpp +++ b/samValidator.cpp @@ -74,7 +74,7 @@ int main(int argc, char* argv[]) { // both mates are mapped if (!(b->core.flag & 0x0004) && !(b2->core.flag & 0x0004)) { isValid = (b->core.tid == b2->core.tid) && (b->core.pos == b2->core.mpos) && (b2->core.pos == b->core.mpos); - if (!isValid) { printf("\nOne of paired-end read %s's alignment does not have two mates adjacent to each other!\n", qname.c_str()); continue; } + if (!isValid) { printf("\nOne of paired-end read %s's alignment does not have two mates adjacent to each other! If you're running covert-sam-for-rsem now, this might mean the read contains duplicate alignments.\n", qname.c_str()); continue; } } readlen = ((b->core.flag & 0x0040) ? (uint64_t(b->core.l_qseq) << 32) + b2->core.l_qseq : (uint64_t(b2->core.l_qseq) << 32) + b->core.l_qseq); -- 2.39.2 From 9c2e46183a19d661f0a618a8eabe8ce1f6a8e2d6 Mon Sep 17 00:00:00 2001 From: Bo Li Date: Mon, 10 Sep 2012 16:22:52 -0500 Subject: [PATCH 07/16] changed output format to contain FPKM etc. ; fixed a bug for paired-end reads --- BamWriter.h | 28 +- Buffer.h | 30 +- EBSeq/DESCRIPTION | 12 - EBSeq/EBSeq_1.1.3.tar.gz | Bin 0 -> 7359273 bytes EBSeq/NAMESPACE | 38 -- EBSeq/R/CheckNg.R | 21 -- EBSeq/R/DenNHist.R | 20 -- EBSeq/R/DenNHistTable.R | 38 -- EBSeq/R/EBMultiTest.R | 336 ----------------- EBSeq/R/EBTest.R | 339 ------------------ EBSeq/R/GeneMultiSimu.R | 111 ------ EBSeq/R/GeneSimu.R | 241 ------------- EBSeq/R/GeneSimuAt.R | 291 --------------- EBSeq/R/GetData.R | 35 -- EBSeq/R/GetMultiPP.R | 6 - EBSeq/R/GetNg.R | 10 - EBSeq/R/GetPP.R | 4 - EBSeq/R/GetPatterns.R | 12 - EBSeq/R/IsoSimu.R | 122 ------- EBSeq/R/IsoSimuAt.R | 128 ------- EBSeq/R/Likefun.R | 26 -- EBSeq/R/LikefunMulti.R | 28 -- EBSeq/R/LikefunMultiDVDP.R | 28 -- EBSeq/R/LikefunMultiEMP.R | 28 -- EBSeq/R/LogN.R | 45 --- EBSeq/R/LogNMulti.R | 54 --- EBSeq/R/LogNMultiDVDP.R | 57 --- EBSeq/R/LogNMultiEMP.R | 54 --- EBSeq/R/MedianNorm.R | 5 - EBSeq/R/MergeGene.R | 107 ------ EBSeq/R/MergeIso.R | 103 ------ EBSeq/R/PlotFDTP.R | 10 - EBSeq/R/PlotFPTP.R | 10 - EBSeq/R/PlotPattern.R | 7 - EBSeq/R/PlotTopCts.R | 8 - EBSeq/R/PolyFitPlot.R | 44 --- EBSeq/R/PoolMatrix.R | 25 -- EBSeq/R/PostFC.R | 28 -- EBSeq/R/QQP.R | 14 - EBSeq/R/QuantileNorm.R | 8 - EBSeq/R/RankNorm.R | 11 - EBSeq/R/TPFDRplot.R | 39 -- EBSeq/R/TopCts.R | 23 -- EBSeq/R/beta.mom.R | 10 - EBSeq/R/crit_fun.R | 15 - EBSeq/R/f0.R | 22 -- EBSeq/R/f1.R | 10 - EBSeq/blockmodeling_0.1.8.tar.gz | Bin 0 -> 67273 bytes .../calcClusteringInfo.cpp | 0 EBSeq/data/GeneEBresultGouldBart2.rda | Bin 2549561 -> 0 bytes EBSeq/data/GeneMat.rda | Bin 249580 -> 0 bytes EBSeq/data/IsoEBresultGouldBart2.rda | Bin 3374596 -> 0 bytes EBSeq/data/IsoList.rda | Bin 163584 -> 0 bytes EBSeq/data/MultiGeneMat.rda | Bin 17921 -> 0 bytes EBSeq/data/datalist | 5 - EBSeq/demo/EBSeq.R | 110 ------ EBSeq/inst/doc/EBSeq_Vignette.pdf | Bin 5805156 -> 0 bytes EBSeq/makefile | 16 + EBSeq/man/CheckNg.Rd | 66 ---- EBSeq/man/DenNHist.Rd | 107 ------ EBSeq/man/DenNHistTable.Rd | 85 ----- EBSeq/man/EBMultiTest.Rd | 130 ------- EBSeq/man/EBSeq_NingLeng-package.Rd | 34 -- EBSeq/man/EBTest.Rd | 139 ------- EBSeq/man/GeneEBresultGouldBart2.Rd | 83 ----- EBSeq/man/GeneMultiSimu.Rd | 112 ------ EBSeq/man/GeneSimu.Rd | 100 ------ EBSeq/man/GeneSimuAt.Rd | 100 ------ EBSeq/man/GetData.Rd | 94 ----- EBSeq/man/GetMultiPP.Rd | 47 --- EBSeq/man/GetNg.Rd | 66 ---- EBSeq/man/GetPP.Rd | 45 --- EBSeq/man/GetPatterns.Rd | 46 --- EBSeq/man/IsoEBresultGouldBart2.Rd | 275 -------------- EBSeq/man/IsoSimu.Rd | 101 ------ EBSeq/man/IsoSimuAt.Rd | 107 ------ EBSeq/man/Likefun.Rd | 50 --- EBSeq/man/LikefunMulti.Rd | 50 --- EBSeq/man/LogN.Rd | 73 ---- EBSeq/man/LogNMulti.Rd | 73 ---- EBSeq/man/MedianNorm.Rd | 53 --- EBSeq/man/MergeGene.Rd | 61 ---- EBSeq/man/MergeIso.Rd | 63 ---- EBSeq/man/PlotFDTP.Rd | 71 ---- EBSeq/man/PlotFPTP.Rd | 71 ---- EBSeq/man/PlotPattern.Rd | 43 --- EBSeq/man/PlotTopCts.Rd | 66 ---- EBSeq/man/PolyFitPlot.Rd | 84 ----- EBSeq/man/PoolMatrix.Rd | 82 ----- EBSeq/man/PostFC.Rd | 47 --- EBSeq/man/QQP.Rd | 81 ----- EBSeq/man/QuantileNorm.Rd | 54 --- EBSeq/man/RankNorm.Rd | 50 --- EBSeq/man/TPFDRplot.Rd | 103 ------ EBSeq/man/TopCts.Rd | 85 ----- EBSeq/man/beta.mom.Rd | 48 --- EBSeq/man/crit_fun.Rd | 48 --- EBSeq/man/f0.Rd | 69 ---- EBSeq/man/f1.Rd | 71 ---- ...seq-generate-ngvector-from-clustering-info | 0 EM.cpp | 227 +++++++----- Gibbs.cpp | 285 ++++++++------- PairedEndModel.h | 7 +- PairedEndQModel.h | 7 +- SamParser.h | 11 +- SingleModel.h | 7 +- SingleQModel.h | 7 +- calcCI.cpp | 124 ++++--- makefile | 10 +- rsem-calculate-expression | 139 +++++-- rsem-form-counts-matrix | 4 +- rsem-generate-ngvector | 4 +- simulation.cpp | 167 +++++---- 113 files changed, 650 insertions(+), 6354 deletions(-) delete mode 100644 EBSeq/DESCRIPTION create mode 100644 EBSeq/EBSeq_1.1.3.tar.gz delete mode 100644 EBSeq/NAMESPACE delete mode 100644 EBSeq/R/CheckNg.R delete mode 100644 EBSeq/R/DenNHist.R delete mode 100644 EBSeq/R/DenNHistTable.R delete mode 100644 EBSeq/R/EBMultiTest.R delete mode 100644 EBSeq/R/EBTest.R delete mode 100644 EBSeq/R/GeneMultiSimu.R delete mode 100644 EBSeq/R/GeneSimu.R delete mode 100644 EBSeq/R/GeneSimuAt.R delete mode 100644 EBSeq/R/GetData.R delete mode 100644 EBSeq/R/GetMultiPP.R delete mode 100644 EBSeq/R/GetNg.R delete mode 100644 EBSeq/R/GetPP.R delete mode 100644 EBSeq/R/GetPatterns.R delete mode 100644 EBSeq/R/IsoSimu.R delete mode 100644 EBSeq/R/IsoSimuAt.R delete mode 100644 EBSeq/R/Likefun.R delete mode 100644 EBSeq/R/LikefunMulti.R delete mode 100644 EBSeq/R/LikefunMultiDVDP.R delete mode 100644 EBSeq/R/LikefunMultiEMP.R delete mode 100644 EBSeq/R/LogN.R delete mode 100644 EBSeq/R/LogNMulti.R delete mode 100644 EBSeq/R/LogNMultiDVDP.R delete mode 100644 EBSeq/R/LogNMultiEMP.R delete mode 100644 EBSeq/R/MedianNorm.R delete mode 100644 EBSeq/R/MergeGene.R delete mode 100644 EBSeq/R/MergeIso.R delete mode 100644 EBSeq/R/PlotFDTP.R delete mode 100644 EBSeq/R/PlotFPTP.R delete mode 100644 EBSeq/R/PlotPattern.R delete mode 100644 EBSeq/R/PlotTopCts.R delete mode 100644 EBSeq/R/PolyFitPlot.R delete mode 100644 EBSeq/R/PoolMatrix.R delete mode 100644 EBSeq/R/PostFC.R delete mode 100644 EBSeq/R/QQP.R delete mode 100644 EBSeq/R/QuantileNorm.R delete mode 100644 EBSeq/R/RankNorm.R delete mode 100644 EBSeq/R/TPFDRplot.R delete mode 100644 EBSeq/R/TopCts.R delete mode 100644 EBSeq/R/beta.mom.R delete mode 100644 EBSeq/R/crit_fun.R delete mode 100644 EBSeq/R/f0.R delete mode 100644 EBSeq/R/f1.R create mode 100644 EBSeq/blockmodeling_0.1.8.tar.gz rename calcClusteringInfo.cpp => EBSeq/calcClusteringInfo.cpp (100%) delete mode 100644 EBSeq/data/GeneEBresultGouldBart2.rda delete mode 100644 EBSeq/data/GeneMat.rda delete mode 100644 EBSeq/data/IsoEBresultGouldBart2.rda delete mode 100644 EBSeq/data/IsoList.rda delete mode 100644 EBSeq/data/MultiGeneMat.rda delete mode 100644 EBSeq/data/datalist delete mode 100644 EBSeq/demo/EBSeq.R delete mode 100644 EBSeq/inst/doc/EBSeq_Vignette.pdf create mode 100644 EBSeq/makefile delete mode 100644 EBSeq/man/CheckNg.Rd delete mode 100644 EBSeq/man/DenNHist.Rd delete mode 100644 EBSeq/man/DenNHistTable.Rd delete mode 100644 EBSeq/man/EBMultiTest.Rd delete mode 100644 EBSeq/man/EBSeq_NingLeng-package.Rd delete mode 100644 EBSeq/man/EBTest.Rd delete mode 100644 EBSeq/man/GeneEBresultGouldBart2.Rd delete mode 100644 EBSeq/man/GeneMultiSimu.Rd delete mode 100644 EBSeq/man/GeneSimu.Rd delete mode 100644 EBSeq/man/GeneSimuAt.Rd delete mode 100644 EBSeq/man/GetData.Rd delete mode 100644 EBSeq/man/GetMultiPP.Rd delete mode 100644 EBSeq/man/GetNg.Rd delete mode 100644 EBSeq/man/GetPP.Rd delete mode 100644 EBSeq/man/GetPatterns.Rd delete mode 100644 EBSeq/man/IsoEBresultGouldBart2.Rd delete mode 100644 EBSeq/man/IsoSimu.Rd delete mode 100644 EBSeq/man/IsoSimuAt.Rd delete mode 100644 EBSeq/man/Likefun.Rd delete mode 100644 EBSeq/man/LikefunMulti.Rd delete mode 100644 EBSeq/man/LogN.Rd delete mode 100644 EBSeq/man/LogNMulti.Rd delete mode 100644 EBSeq/man/MedianNorm.Rd delete mode 100644 EBSeq/man/MergeGene.Rd delete mode 100644 EBSeq/man/MergeIso.Rd delete mode 100644 EBSeq/man/PlotFDTP.Rd delete mode 100644 EBSeq/man/PlotFPTP.Rd delete mode 100644 EBSeq/man/PlotPattern.Rd delete mode 100644 EBSeq/man/PlotTopCts.Rd delete mode 100644 EBSeq/man/PolyFitPlot.Rd delete mode 100644 EBSeq/man/PoolMatrix.Rd delete mode 100644 EBSeq/man/PostFC.Rd delete mode 100644 EBSeq/man/QQP.Rd delete mode 100644 EBSeq/man/QuantileNorm.Rd delete mode 100644 EBSeq/man/RankNorm.Rd delete mode 100644 EBSeq/man/TPFDRplot.Rd delete mode 100644 EBSeq/man/TopCts.Rd delete mode 100644 EBSeq/man/beta.mom.Rd delete mode 100644 EBSeq/man/crit_fun.Rd delete mode 100644 EBSeq/man/f0.Rd delete mode 100644 EBSeq/man/f1.Rd rename rsem-for-ebseq-generate-ngvector-from-clustering-info => EBSeq/rsem-for-ebseq-generate-ngvector-from-clustering-info (100%) diff --git a/BamWriter.h b/BamWriter.h index f664710..bbdd298 100644 --- a/BamWriter.h +++ b/BamWriter.h @@ -125,24 +125,24 @@ void BamWriter::work(HitWrapper wrapper) { bool notgood = (b->core.flag & 0x0004) || (b2->core.flag & 0x0004); if (!notgood) { - //swap if b is mate 2 - if (b->core.flag & 0x0080) { - assert(b2->core.flag & 0x0040); - bam1_t *tmp = b; - b = b2; b2 = tmp; - } + //swap if b is mate 2 + if (b->core.flag & 0x0080) { + assert(b2->core.flag & 0x0040); + bam1_t *tmp = b; + b = b2; b2 = tmp; + } - hit = wrapper.getNextHit(); - assert(hit != NULL); + hit = wrapper.getNextHit(); + assert(hit != NULL); - assert(transcripts.getInternalSid(b->core.tid + 1) == hit->getSid()); - assert(transcripts.getInternalSid(b2->core.tid + 1) == hit->getSid()); + assert(transcripts.getInternalSid(b->core.tid + 1) == hit->getSid()); + assert(transcripts.getInternalSid(b2->core.tid + 1) == hit->getSid()); - convert(b, hit->getConPrb()); - convert(b2, hit->getConPrb()); + convert(b, hit->getConPrb()); + convert(b2, hit->getConPrb()); - b->core.mpos = b2->core.pos; - b2->core.mpos = b->core.pos; + b->core.mpos = b2->core.pos; + b2->core.mpos = b->core.pos; } /* diff --git a/Buffer.h b/Buffer.h index 5017796..3e45094 100644 --- a/Buffer.h +++ b/Buffer.h @@ -12,12 +12,13 @@ const int FLOATSIZE = sizeof(float); class Buffer { public: - Buffer(int nMB, int nSamples, int cvlen, const char* tmpF) { + // in_mem_arr must be allocated memory before the Buffer is constructed + Buffer(int nMB, int nSamples, int vlen, float* in_mem_arr, const char* tmpF) { cpos = 0; - size = bufsize_type(nMB) * 1024 * 1024 / FLOATSIZE / cvlen; + size = bufsize_type(nMB) * 1024 * 1024 / FLOATSIZE / vlen; if (size > (bufsize_type)nSamples) size = nSamples; general_assert(size > 0, "Memory allocated for credibility intervals is not enough!"); - size *= cvlen; + size *= vlen; buffer = new float[size]; ftmpOut.open(tmpF, std::ios::binary); @@ -25,7 +26,8 @@ public: fr = to = 0; this->nSamples = nSamples; - this->cvlen = cvlen; + this->vlen = vlen; + this->in_mem_arr = in_mem_arr; } ~Buffer() { @@ -36,14 +38,13 @@ public: ftmpOut.close(); } - void write(int n, float **vecs) { + void write(float value, float *vec) { pthread_assert(pthread_mutex_lock(&lock), "pthread_mutex_lock", "Error occurred while acquiring the lock!"); - for (int i = 0; i < n; i++) { - if (size - cpos < bufsize_type(cvlen)) flushToTempFile(); - memcpy(buffer + cpos, vecs[i], FLOATSIZE * cvlen); - cpos += cvlen; - ++to; - } + if (size - cpos < bufsize_type(vlen)) flushToTempFile(); + in_mem_arr[to] = value; + memcpy(buffer + cpos, vec, FLOATSIZE * vlen); + cpos += vlen; + ++to; pthread_assert(pthread_mutex_unlock(&lock), "pthread_mutex_unlock", "Error occurred while releasing the lock!"); } @@ -51,11 +52,12 @@ private: bufsize_type size, cpos; // cpos : current position float *buffer; + float *in_mem_arr; std::ofstream ftmpOut; pthread_mutex_t lock; int fr, to; // each flush, sample fr .. to - 1 - int nSamples, cvlen; + int nSamples, vlen; // vlen : vector length void flushToTempFile() { std::streampos gap1 = std::streampos(fr) * FLOATSIZE; @@ -63,12 +65,12 @@ private: float *p = NULL; ftmpOut.seekp(0, std::ios::beg); - for (int i = 0; i < cvlen; i++) { + for (int i = 0; i < vlen; i++) { p = buffer + i; ftmpOut.seekp(gap1, std::ios::cur); for (int j = fr; j < to; j++) { ftmpOut.write((char*)p, FLOATSIZE); - p += cvlen; + p += vlen; } ftmpOut.seekp(gap2, std::ios::cur); } diff --git a/EBSeq/DESCRIPTION b/EBSeq/DESCRIPTION deleted file mode 100644 index 5f61713..0000000 --- a/EBSeq/DESCRIPTION +++ /dev/null @@ -1,12 +0,0 @@ -Package: EBSeq -Type: Package -Title: A R package for Gene and Isoform Differential Expression Analysis On RNA-Seq Data -Version: 1.1 -Date: 2012-4-18 -Author: Ning Leng -Maintainer: Ning Leng -Depends:blockmodeling -Description: RNA-Seq Differential Expression Analysis on both gene and isoform level -License: -LazyLoad: yes -Packaged: 2012-04-25 05:25:10 UTC; ningleng diff --git a/EBSeq/EBSeq_1.1.3.tar.gz b/EBSeq/EBSeq_1.1.3.tar.gz new file mode 100644 index 0000000000000000000000000000000000000000..f0fd1214bc2af4cb0d124da8fdea8bf9378abd49 GIT binary patch literal 7359273 zcmV)2K+L}%iwFP)fIUzE1MFOTciT3y&%frUK$YY~+)J|(X~mnK%RMKy%6E-^FFLKi#@V>*s;3`gCh+N&7dU zZ){Qfd!Fa}o(Ju})$siXtg(Jk=CA+q_K%}p7G{@qYJA&z$b7eH4ZeYhz=B?WV~fDxx&WgD`i8 zVZTVTdd4U9e!aHatU0w^uVyYu~isZ;BA+nw6k-vakp%RZyxvIMBX zcRitn{XBgU<&wO9fE(Z($md7mQDq@*!0#b4HSj*oZn#|NE2uK89Kayaa8 z^b5!K2{bKg=#8&vXjScMTSKbZ+)60jB)ZY%J^ZvgWvh!++3@tRuf>SddyBy z6^Wiid(ka5nwiDXjdu9+ix-Yjh4NOk;s!IW=W#6M!oLnP+f0j9Z8S<2%_eW_N;=yT z4I|Ta-bc_f*a@M=D&^q0NGem5_N<3nV5o9lxH;#8degTptRw!*!h)0SIgZBDdi$4l zA8J*ENfB&03?APVO;LUE;E4xSYeOf?IxY15={Ox8hvP}i0}xm|BTIOKOA~8%FtMAZ zlfZ2{tP6cm$s#0x6c^=@Z(S6}eWuGF@GzMv@iK`oX&j!28f#0ds{ROjFTH{>$=e-a zXDT_w3vJkuL}})YQcKLe$iLX$=+P4C=mpJD&a_QJp6;O9tyNOXnmKC1wlX$ z2(M*~i0ziy@tHRLs?ORY4^p=u#&JFCMajT1a?q%-7rSWFRcK?c-ilEsGETz0;B{(| znnRA0YL@d!PeEO?U+d8e(yM0+EmU)P(j%)}zI`(g;4pN` zC+8W5BGTrHa_fT!OHs<72ReKeTtK@)w`euY?k229sRwC0^dp%M{j|t}bkq29FDiZs z<0%i`lfLS;7}~1lKyyN4)gu-q3|hms3~y}~J3Wog&H|azjnh&Ck9M2rTPTAZc^h+X zA1|XoCkLvTDAE%00WheLpC%~aEEgF2^DhK78!5$60*a5*QL`~GVDs$kv|1-fd6cJF zafYL1w!-*w7!}KkIQ=vZYL`5Y(~FuDUkqT(t>*8|@%L8o_ZIT^s{AKW^7dxj{WvZaj4Zz?MW{NEJ$D zd;qVMZe=AUH`|Ge%b-b2eW#axa3%X{zxSo*RC)mV!U>krF!pyC-V?NN${v$6(jZG~ zT(cLE*8bI?zrnto1NdS8w&DO7PBKSbCE6>Zg9;iTuv~`V0-*yUSd8~)(7iQl2uw%+ z)Tu&U$3VFKy1VUwYp4J1(yEx|L?e*DIvRjHb4;gi_}8KmU6dSK0o%=e+;+aO55h z*3Dz}{BOf=`P;Mezy7v&KmU6d*K0IFSF-J|t=C~3g?ah>coVE%p0dMEr_I8IJrNn6 zFz}VE^L+P5mKNtcV}Kr_J{zR{>6j-4A23MduJyWzikM$nlHPHVP5Z?(;}&ol@VuWz z6EwD07C4ai(SS1`#(H3|*#rH z1skE~b{W#2=HaMhT|eX(2&o|W;1>Wq9gs2`JM2Ro!QV^xYg<&^ad_DS`dH~FL_*Y6 zR5R`U4t^tFP|+|OQA&PgG4?tt`1ne4yuPxYGep0DlaeBglANJ@qyc!0p~U0IjD7sb z_J=IteV*rGb_s4I%{Y@;Krc{FRF0uF7WUHj910?>_Xzc_tSvSH;E=P62p}BPBjcd% zgb!SH0)>ZB24eUu%!H=V7m_Nq+ELdCNEzU!F@rD@FdOJK@F*MH_)^6d(%UM|i64t{ zB9ZqY^eHlt22BQ$xaQ+htzTz+2nvHHglc1Vn-=G3c4ZyH=wmNn_(3A7K196`d|=ZYR0|qO#|~xLfF3eP z8FYqtQpYSolnnR>n^E5()fA9D%mzZ|fYx#s-15VRYG#Kkzzte>&trxj{Oe|QR!zV) zA;MffmBaKbj%2WU;-;L$+%l$Yo*FYYwO?zDLSdS;n4OhDcmnu{2}qVxPdl0>j9inI zjNS(%u;ZKmY^&#s7Y+EY>F8g0(JEi4m-GTqUmyd~Je*A8OMvswoZ&Q~K2(ovR1WZi z2-g`zW3}8v5w8s~qz#aN2dtle12Q1pr^vK!LZ-En3{763hUezJ^fXUrC*yf=xPuc8 zogKFajNNGxy_@pxn{JxMU<;2Y5IFoid2`U->!#70!xV8fr0Yo=n$#_4ZQ#bC)@pSC z<%+6^1YafizkEjjGjTM&hsMJBLK8=0AHk4Fn{!n&87Vu=NDxyCPRn8a`exH<*$(|~ zqr%Y0v*&{WOJOj>*d!~f2j^RUd68y=9I!n{a|xE<$v-;3zIu7oK7N9J28J?}Khh3> zrexV5%1ASp;0Pp%EA&R!3nBM6b(Q$Ev)kze9M#%k zn_x<)Z5$R^^r6IU=*TCG-scRbWdK7)z-5x9y|5R>VC{Hr6Jr9T1z-eGL2_!8Hs*kM z2Cj9KUgQMCF!vxD5-cj@YvBa@X*?Y#i*l8O0Y+SF(Ysv`8-W>iWPnGYa8R9?J7&2C zJcLI;hb$CTN~Yr;YH6rBf~6B`kSY*xPo%&*Aj%`sYy)j*)3pG)A|=FRB9~wdxDen( zR1i+jbrw-cR#(qidVvK5OB>aR!N7n`HifQN(QkeY{lzUv@q^N-?r>7IG--p>khswf z?BYBHbdDeq@N`PRvST)m=&p*2Rym9vP&vw1P-Y*^XhKMh4RV=IQKv$ZBrUjl zIhy3*Y2=&*8rqxm9IeRu?>5fhey5io zk&S(de@~s?&;B=v&P2Vn`fPgx$U{9eWx>>-8QHhM6_kfPq@IhKR@4J6xlxQKfoPc= z8ev1wg+#ktT9s|EfooBi9e*hcq@fa7A8C79iyNfHK7jGF0B-92%_am3Xa57y1lUh7 z-2M7@g0iG&pZ4bjw>GY^*`Z%R!{fw6UFvJE?G%?29@OHRgY!QSGeWFrJ0Ie(7t|1_ zv^&uKV;H`;u*EJZ9Sl|tXg*ZSW?3{kFM4TFq~jVbua)gVD_;63(bC*m1oH9i#h?j6 zNHxV{*N@qoH!lytj=VhVym_N|vI$r{oc#Irqn7(?T~`_Z(O#0<2>(76|9RWn{@niW z?bdz#cNdq4|ETSC_$Qy%7X692J^~kqBN>j&MSqL4X;c!)<#{7W)SI7Z;}TTEVQTNw z+Ji!a3q^)fx(b$+A__}y0DmQvg#L)Lt5Jli8=7~65yT<&68j0$FXA>Z);zO;!L@Zu zXuTr)=o%BMVT+cTf()vYQU=LYC4gJJG?4-xaG?O#6Lgq@6rjFJLAjzP#3!-+k|~Bs zgv5b@@7j{%Ok0&=8+=qFsUovww4}>Ll_Xq2!}_CEH7cRPAHGuX`TITD4YzvUbZh7G z?VQ%(NQ~D{9hxZT%Wt@Aam{nwf56;y>HiOSHsVmV^+SM9`TwU+=jMN(Htx@V-N_~V zziN980m_6mp#ZhL4y`D8i(Kz$XE5NVaWRJoN>8t(!I68se~{v?dJ3w#A(+d6)!grE zx8d>yFX^}84(FsmeraHZu`U-m_iE;JpP1|WF}(nhl1o-Z=zp{)Egk6>R4N2Cjynnwga0fH#5|K&e6}x z`R2fFw9x*R6Z~C0J$y-fz}Ky7@HDFoM2ZA zOoO7T=rkQ|NB(=}XG;G+Y~D)n|KJAk|5I@c+Z6_WXFeCZ9iS*7uVM zI%_)ovJZF}ccI{KXFQ3r2=M$74fmHr+Hw!0Do#iC91;LR|MzAT*}3rb2A ztSG8-%$F5fTTxVIm@iq;^9lmhe5>iKs(v-!Je{&-R}e8W%~OE>u%alwEz8^vKqE?; zZ=sju)h%IBbTdc}6l2taQGTINV(hI#)3lC3C)Z|BRLtSob~ zuuLL%J1r|E-rxJW`+Z9P0sEBZhP@p1I17gpaca|ue*Zyj6~QyE%-A@>%Izd{wi(YB z(TDqCH^yxORpIMe&YLo z9wg4-zn16ShSpUWP*?QNrpASVU0Gi%Cx11 z-1zK?3;tmOs>`vJ$wBNWSGO>rPtzLUpAZgN6|6bM3NvzbaMyQ!%bz-$rj zpAl5b;YnH%ss(m&nv4*6(zy#uOFiJ zbO%(5?qRzEi{2$MH9p*;Q6{N>W^kgO!6M`lVqa~n)T$)f4%D7jpD~uz2@&W7U}OI! zFqv3BLn)WOD`gIq`l(CBA+EMhcpsF|mn{kX`DKQolTk+2yevZp!^=TNB!(u5vh|hq z0of1BEn*Y`q2fhn%crA}XfAn>C}%CVY;mw#vg&Xrv(PE@qhm0=V)b2*ATk!`&^wb9 zqOVCTyM{Xh*RbgLB5vkG%o)J?0T?9ZRHu-@Q6$CJ`->Hn7KruX_xNQO%fwF7u$N&o2 z2)_7LDDoeMj%h)iIEZb0)SdEZ9L0dL=#-XKq~J2qQGh~IFk*Q4s28DAe1(~3oO6M* zIMB-FKEPJv^6v{ay2Vy~$sX|g$|G1{eS zOcL!aaN(n9!Qds8Ap#Q}g6BV2xagPvdpw07_7Pwkt7mZAV1Hrp8|drraGB7;zyA(# zQOg0%7vVn--;A7|C_*r=E?*|*OWV>gfu?ZYPc9k6=bu3^;l-JCw3t(~&Y~XLekUT<*0#=Nh`thlF@tei2*wljBu_zj?UARFcu5|2rP=|J+BzJ&tZH$7q=ME z#K2s}P{F|U4G#{;!ADt9c zJ)%9%OhHeSpE%4`tq&) z=i3$W-*c7o|8j@#-GBda+uv@?p8wctZQsxT-^C^7|7F|jxqtn<-qRm)X0qJmF*if6 zzlIl8N*EZGa0ShzgzS)yLh{;R)x@_V%ZUTP0gp9p5^9N(FOb2^*4ruzoE;U3!Q0S zCiC8GN$o2%TPe0_nuZn{LYX9!WZGnAIx}gyP>_8ULD^J3h=QQ3Dgw%;fGmQb>>xV= zf*`xd{?9r0-uLc%^Cp?JNt(huXy(0p&pmfP_uO;OJ;w|Z3hfZkI_jD!Qd6Krsh0_9 z2C>NbGzk~XDs-c8sJe%gn5Vd80b46`=Pw~&iF&!LF#-SFZbfeUWRA96T*#Q{ zVw-=+G|Od_n138q1bP_bY97=87TBY7dvHk1LC~q(-X!sPxTF!sIe@E~(t%5<;xk1) zeoF#sX3dNpxkxL4_o#`gN}wj@4`8+yC&&f(=9QDKbbsF(XrOdDm*}!&H3?S!;pOIN zim5k_>l+MPx&zdk9+ai#Q>EQpz~jSl!R-1BrYkKry%W#7px`;Ed^(GDNu>wlRYV5!eFeI+GY0su%E=FS0%D8wXmAdY6ObgZ`AOLW$e0rdWS)^_4p1MKN+3K16hfX~<~o3b$3Av!VhYjN1tUj7CRPv& zmD#3rrpTsHuRPkNgIOAb@*j+3#dQ-ZbqAhy2W4Ty>1tOs{$0{Km zNC6=5s{FJ;RT5o2i53-Qhld0exoV6 z=50-n+7l=gYIliK z{yr`qDwH}F-a1v-jZR7;EdX(AItvON@X)0VZHCg+qp^1#-*xsL!uJr(S~~tEsf$q6 zw2Vsdu_>EfV~Uu z2RKm`S3f;hZOLRuS>Cg`JeK^Yb|Bf;)tXHtMqL6l-2GoYT;TsB2H<7+&p13x{$uvN zDFKo~Mo9*=ZxYB`XD{*Rb$tVL&;#VUBaR$!nw$H_H7}lu_pQk#b3+-39&OPFV^4;> zTI+IC#(pSkhY|Hm;Hb)MmTbro7M|rva3_mvepW;AxT%|j+@c1>d+PD9n%4p^2Elx=cwwyC_OFF&O`^8P=SAmit- zZz45`+ITjvj|N9LjQt1rFY1v0$MmxOukm@r{@<4GJ4KJNt)$KKJV26%jlWkLQKk z2=)a_*$j`k)uP2b$`r-~!g@Rd59O6MZAJt&zLa1frdQM7F#gqQ=x>bu1)PCEr2@O_ ztZ$mCc)ODQBnR;N4S`%I82En*VD9YA66SWr%h?j|FmJdRN(k`877I%^{0a-J1;Szc zQ)B;Ay%FVPJdk&&mm+0s=C>9K8dQr^V_ODeQEbLwBt&1~DE$rT8&TM1Dp$eZVg+RP z7F`P%VFiq^0@6=7QZouAq&&Sty%Z_!QnNXmNI`{{S&M`;9SaS6Zfe1x&RP;lFZ@Qe zNQ8ZdYQiRNsqmuYqLdt=x1niKOVffmoSsLtYR&5s7<2~eF=#!vz!uUV9@L{Tn4&cSO;@6NAf!}=0>K(36s``$suev5eH4;Ae~3zP zTpQCz1@-40$7>}nvX)pYy}Bk4jfRv^PzQ2@<1JKEt<=N_bbE4FZ-?1&V z$0D&B7;D|rwCiL0(ALZ?G~6TF41p@sn%>=Ua> z`Ut8N3{@k6hJMnpDq)yZS|An;<02E*U?rR=HDec=ORgCt4gi?|`T#XiYJeOnnjQ*h z6Q@4=p(91tR@(Negd)(iPQv$r-%Xv5s3w&Kr3`D3?zQS)Q#VTFLelpqp?^3`6Mj{ zlS?P2J_Lk9(<)D*Lz0W0y=Ut1+~x9Ym?3X*rA| z7NR}jW6rkuRFw^~n{Su-4Yj$&!%KD`k%JyA0;GJKJU|tV=4~z{eSC+V8--)Y^0SE^ z$^TDqlllJ%Y6bp(qGkU-(`{>-;KxWSHyCd4`J zruN{op(nsngO=oUCs+2mQ`U^0NE%AnP5w^`8Ep@cmNOAZfre$!sLQIH0vS#Smyl-FePnP5 z`N7(ys={okpg~r&`;UcXR{2KX4j<(@#EQ-VT)3&)r9&}M5bNS0lADMr?j(`{6sE+6 z7$O}-yLhpzAJU6ueM^}`DRU?jP!u#`&yqtE4>G9ZU+1X+AS0Jp=_55fH5iDf9cv8w z@Ovui8(N!N$QX1`Wv2Li0PlyGhT38=^$m#LJx#F~%m?;HZgDtGF|#3}F~=cO7mcpK zCif+g%xg2`7-*OZYQP0R#tk}Y<_+p^W4|*E_}i9jZwJU^x8mcep@x3L2>KycDJk@U zzYYG(GGGqo1JFmXNKnTxB!0w_=tFCnfC^j5@b2p^iUbl%=~%)dpn!`CN%f&tgiix4 zW}DSJXwPKopgq|LfhMWotu5dfVkZ!9a{)AI2yAvhcXYgXSXRR-XBBT$d)WxcDr6xL zjX;J;;U5HQ&{J+ z=2ckbS@wJ^s(iSqd_<{y(D=(va5wU?-~U0obF)W%131k7FBXe9?|(;vvGV!3^UfdwOek??(3UxP;l*Ub z?7^evY(XN)O1!aI+K5-s>NA=%md#y;)dE0f8MQaQDq+J6TdA-TXCCkJgK5*)3r9R6 zd)d6Fxyf5MYuPs2;bn{V?)Iuj!=$}3hBocFeRvlm_BcOnY1qgsV)x#f&#u^u0*X~739IZY!sFsBLL$JdoUR7+vC+oLdVZnteA#xsFr5!Wd^AnxkB;U1 zqK5Oi(Y8`e=bN_f`7q_lZI7$F9+lfY$@O)uJ;u!StDR;##*75I4Eob@!hU3yjF zFSGM@VLkSeX~f=rQ-zE@H)1dErVQs^$|y8mU&MC3=_E2JZ?HyezFuCU45#}oUYbCS zss*4xP!?+m8Bb&P(n;v3$)IfZU$IK#@Bc5X$0APs&u~yL%YVk>G3~$BGZQYR?N|0a zN|K+JM4Zm@TNX9USD5J6^u{z!r~14~zDzvdtI|zRjs57r1+K{?MqN~b8HtTy9JgZK za`kia-V+PML3!iYmLa4Xoc)(PMLQ7aBSoM~2(SSd_sEVG=D;$Z<$Wd=Lb1RRehYN2 z1Jv%gkwcHwQNC2+By&aPO$~6a$jd5>Z4Ub#$thf1=AamxBW=Y}3AV%|%Z`C~19Lv( zg~agsRNcP#OrOycF<`WSkpe~u)_G`gXL}o;7adS})oQ<|5#|w1EMgHntynebn+~&S zZ%F6aVZWzh9Y39>xq_1gk3K6bJe@xCU1f)NQjsuN^)BNxYFTO^7l79Es<}kQ@Ou5w z8g}e+y#9)jqsMts1h-VdPM`I`QAt=tDm=0s-B!?tZH0n17&&hDCMc-^^uh-4c(}d- zj3K{arD82ZEFiT^HE$Upz!?KcXn^L5Ro*s)ZC86CDv?E2lfmlfW9^Iu&H6`ub_v@-uc z9*@cY^S+n20A|o=dw*bj6?uRH1h8T6h<7A`Ra#?;IrG2`FwU>=B3Nl?*zIyc#$6h? zAh*0v<;i9x)wng0%X6+WVBbLS2Jto2LiDI_7znGKi-!`+g_0T2G2*ECVMa5PhGnRA zQ6SBmHpw;>g^g>9)8^0SY5O59aizd}i5{?I{2htDG#<3%(;}hiY=lxWO z3ilaUwh6O&Ja+vr*yxVG{4W@d72N+0Yi0T0xIC8r7xcYN`HOuQxA1J8jxr@IBwAbn zhY6U7&JX3!ada&@LwRKr&YPueFn4oaShKORyJgLPdRvBTDlVy&nRw?4Xy387{|@tB z3p(w;Ytgd*pK*EE{CDqrlN^jFS|S=2Ag3DyIYVh_nq5iGUX$H4h4Ofw-(U%4XLE{y zjClv}bBvVWLA71FTtHD~yXXm20>6d|An7hhF&NJ&9RtZeQhp%9w3h87ph#FE(yQ?- zAz5&j!)Sd+d<}5F{9+;&reQwt{XJRWDbxLSV}u6*M+n3*yRNSf82V<`3^NoAR_v{H zP-oVIZN(!wD26UhWU(=Lgcb;5U_B6A9s(0as|-dfwGc+kOLitwIT5!6Pel4y5ORpg zyfOgHWp-}@jjvrZFCFjVuuv%sUGL6Wd*%b7B4U-Dp_svnGBbsfWuL)ht}~G68kpt6 zrkflh9|np{>6|ywH}yRDtUea~uTP|!_D$yUB(NM+E^xT}-+C