Skip to content

Commit

Permalink
Fixed bugs in samValidator.cpp
Browse files Browse the repository at this point in the history
  • Loading branch information
Bo Li committed Oct 2, 2016
1 parent d7d3273 commit 591ea2e
Show file tree
Hide file tree
Showing 2 changed files with 14 additions and 8 deletions.
8 changes: 6 additions & 2 deletions rsem-calculate-expression
Expand Up @@ -217,7 +217,6 @@ else {
if ($faiF ne "") { print "Warning: There is no need to set --fai if you ask RSEM to align reads for you.\n" }
}

pod2usage(-msg => "Forward probability should be in [0, 1]!", -exitval => 2, -verbose => 2) if ($probF < 0 || $probF > 1);
pod2usage(-msg => "Min fragment length should be at least 1!", -exitval => 2, -verbose => 2) if ($minL < 1);
pod2usage(-msg => "Min fragment length should be smaller or equal to max fragment length!", -exitval => 2, -verbose => 2) if ($minL > $maxL);
pod2usage(-msg => "The memory allocated for calculating credibility intervals should be at least 1 MB!\n", -exitval => 2, -verbose => 2) if ($NMB < 1);
Expand Down Expand Up @@ -313,6 +312,7 @@ if (!defined($probF)) {
}
}

pod2usage(-msg => "Forward probability should be in [0, 1]!", -exitval => 2, -verbose => 2) if ($probF < 0 || $probF > 1);

if ($paired_end) {
if ($no_qual) { $read_type = 2; }
Expand Down Expand Up @@ -1135,6 +1135,10 @@ Parameters for all the above models are learned from a training set. For detaile
The options in this section are deprecated. They are here only for compatibility reasons and may be removed in future releases.
=back
=over
=item B<--sam>
Inputs are alignments in SAM format. (Default: off)
Expand All @@ -1143,7 +1147,7 @@ Inputs are alignments in SAM format. (Default: off)
Inputs are alignments in BAM format. (Default: off)
=item strand-specific
=item B<--strand-specific>
Equivalent to '--strandedness forward'. (Default: off)
Expand Down
14 changes: 8 additions & 6 deletions samValidator.cpp
Expand Up @@ -80,7 +80,9 @@ int main(int argc, char* argv[]) {

printf("."); fflush(stdout);
do {
if (sam_read1(in, header, b) < 0) break;
int ret = sam_read1(in, header, b);
if (ret == -1) break;
else if (ret < 0) { isValid = false; continue; }
assert(b->core.l_qseq > 0);

qname = bam_get_canonical_name(b);
Expand All @@ -98,15 +100,15 @@ int main(int argc, char* argv[]) {
if (ispaired) {
isValid = (sam_read1(in, header, b2) >= 0) && (qname == bam_get_canonical_name(b2)) && bam_is_paired(b2);
if (!isValid) {
printf("\nOnly find one mate for paired-end read %s!\n", bam_get_qname(b));
printf("\nOnly find one mate for paired-end read %s!\nPlease make sure that the two mates of a paired-end read are adjacent to each other.\n", bam_get_qname(b));
continue;
}

assert(b2->core.l_qseq > 0);

isValid = (bam_is_read1(b) && bam_is_read2(b2)) || (bam_is_read1(b2) && bam_is_read2(b));
if (!isValid) {
printf("\nThe two mates of paired-end read %s are not adjacent!\n", bam_get_qname(b));
printf("\nThe two mates of paired-end read %s are marked as both mate1 or both mate2!\n", bam_get_qname(b));
continue;
}

Expand Down Expand Up @@ -138,7 +140,7 @@ int main(int argc, char* argv[]) {
}

bam1_t *tb = (b->core.pos < b2->core.pos ? b : b2);
isValid = tb->core.pos >= 0 || tb->core.pos + abs(tb->core.isize) <= header->target_len[tb->core.tid];
isValid = tb->core.pos >= 0 && tb->core.pos + abs(tb->core.isize) <= header->target_len[tb->core.tid];
if (!isValid) {
printf("\nPaired-end read %s aligns to [%d, %d) of transcript %s, which exceeds the transcript's boundary [0, %d)!\n",
bam_get_qname(b), tb->core.pos, tb->core.pos + abs(tb->core.isize), header->target_name[tb->core.tid], header->target_len[tb->core.tid]);
Expand All @@ -163,7 +165,7 @@ int main(int argc, char* argv[]) {

if (cqname != qname) {
isValid = used.find(qname) == used.end();
if (!isValid) { printf("\nThe alignments of read %s are not grouped together!", qname.c_str()); continue; }
if (!isValid) { printf("\nThe alignments of read %s are not grouped together!\n", qname.c_str()); continue; }
used.insert(cqname);
cqname = qname;
creadlen = readlen;
Expand All @@ -172,7 +174,7 @@ int main(int argc, char* argv[]) {
else {
assert(cqname != "");
isValid = (creadlen == readlen && (!ispaired || creadlen2 == readlen2));
if (!isValid) { printf("\nRead %s have different read/mate lengths!\n", qname.c_str()); continue; }
if (!isValid) { printf("\nRead %s have alignments showing different read/mate lengths!\n", qname.c_str()); continue; }
}

++cnt;
Expand Down

0 comments on commit 591ea2e

Please sign in to comment.