[med-svn] [Git][med-team/ivar][master] 8 commits: New upstream version 1.4.2+dfsg
Étienne Mollier (@emollier)
gitlab at salsa.debian.org
Sun Jul 14 10:17:51 BST 2024
Étienne Mollier pushed to branch master at Debian Med / ivar
Commits:
c446e8d5 by Étienne Mollier at 2023-07-20T14:26:49+02:00
New upstream version 1.4.2+dfsg
- - - - -
4e892866 by Étienne Mollier at 2024-07-14T10:46:58+02:00
New upstream version 1.4.3+dfsg
- - - - -
be4f3528 by Étienne Mollier at 2024-07-14T10:46:59+02:00
Update upstream source from tag 'upstream/1.4.3+dfsg'
Update to upstream version '1.4.3+dfsg'
with Debian dir d9da123ac01f2e8351a8c711768912aacbda1b27
- - - - -
d9ca5540 by Étienne Mollier at 2024-07-14T11:03:59+02:00
forward-build-options.patch: refresh.
- - - - -
7b379c20 by Étienne Mollier at 2024-07-14T11:04:15+02:00
gcc-13-2.patch: refresh.
- - - - -
2f7f265a by Étienne Mollier at 2024-07-14T11:14:52+02:00
d/rules: move permission fixup to dh_fixperms step.
- - - - -
49489b2b by Étienne Mollier at 2024-07-14T11:16:03+02:00
d/control: declare compliance to standards version 4.7.0.
- - - - -
de6dc239 by Étienne Mollier at 2024-07-14T11:16:43+02:00
ready to upload to unstable.
- - - - -
20 changed files:
- Dockerfile
- configure.ac
- + data/test.strand.pileup
- + data/test_strand.gff
- + data/test_unpaired.sorted.bam
- + data/test_unpaired.sorted.bam.bai
- debian/changelog
- debian/control
- debian/patches/forward-build-options.patch
- debian/patches/gcc-13-2.patch
- debian/rules
- src/ivar.cpp
- src/parse_gff.cpp
- src/ref_seq.cpp
- src/ref_seq.h
- src/trim_primer_quality.cpp
- tests/Makefile.am
- tests/check_quality_trim.cpp → tests/test_quality_trim.cpp
- + tests/test_quality_trim_unpaired.cpp
- + tests/test_variants_reverse_strand.cpp
Changes:
=====================================
Dockerfile
=====================================
@@ -1,18 +1,18 @@
-FROM ubuntu:18.04
+FROM ubuntu:20.04
MAINTAINER Karthik G <gkarthik at scripps.edu>
RUN apt-get update
RUN apt-get install -y build-essential autoconf zlib1g-dev python3 wget libbz2-dev liblzma-dev libncurses-dev git bedtools python3-pip vim nano
# HTSlib
RUN cd root/ &&\
- wget https://github.com/samtools/htslib/releases/download/1.9/htslib-1.9.tar.bz2 &&\
- tar xvf htslib-1.9.tar.bz2 &&\
- cd htslib-1.9/ &&\
+ wget https://github.com/samtools/htslib/releases/download/1.10.2/htslib-1.10.2.tar.bz2 &&\
+ tar xvf htslib-1.10.2.tar.bz2 &&\
+ cd htslib-1.10.2/ &&\
./configure &&\
make &&\
make install &&\
cd ../ &&\
- rm htslib-1.9.tar.bz2
+ rm htslib-1.10.2.tar.bz2
ENV LD_LIBRARY_PATH /usr/local/lib:$LD_LIBRARY_PATH
# SAMtools
RUN cd root &&\
@@ -34,12 +34,12 @@ RUN cd root/ &&\
make install
# bwa
RUN cd root/ &&\
- wget https://github.com/lh3/bwa/archive/v0.7.17.tar.gz &&\
- tar xvf v0.7.17.tar.gz &&\
- cd bwa-0.7.17/ &&\
+ wget https://github.com/lh3/bwa/archive/v0.7.18.tar.gz &&\
+ tar xvf v0.7.18.tar.gz &&\
+ cd bwa-0.7.18/ &&\
make &&\
cd ../ &&\
- rm v0.7.17.tar.gz
-ENV PATH /root/bwa-0.7.17:$PATH
+ rm v0.7.18.tar.gz
+ENV PATH /root/bwa-0.7.18:$PATH
# Snakemake
RUN pip3 install pandas snakemake
=====================================
configure.ac
=====================================
@@ -2,7 +2,7 @@
# Process this file with autoconf to produce a configure script.
AC_PREREQ([2.63])
-AC_INIT([ivar], [1.4.2], [gkarthik at scripps.edu])
+AC_INIT([ivar], [1.4.3], [gkarthik at scripps.edu])
AM_INIT_AUTOMAKE([foreign subdir-objects])
AC_CONFIG_HEADERS([config.h])
=====================================
data/test.strand.pileup
=====================================
@@ -0,0 +1,20 @@
+test 1 G 20 GGGGGGGGGGGGGGGGGGGG ????????????????????
+test 2 A 20 AAAAAAAAAAAAAAAAAAAA ????????????????????
+test 3 G 20 GGGGGGGGGGGGGGGGGGGG ????????????????????
+test 4 G 20 GGGGGGGGGGGGGGGGGGGG ????????????????????
+test 5 C 20 CCCCCCCCCCTTTTTTTTTT ????????????????????
+test 6 T 20 TTTTTTTTTTTTTTTTTTTT ????????????????????
+test 7 G 20 GGGGGGGGGGGGGGGGGGGG ????????????????????
+test 8 C 20 CCCCCCCCCCCCCCCCCCCC ????????????????????
+test 9 C 20 CCCCCCCCCCCCCCCCCCCC ????????????????????
+test 10 A 20 AAAAAAAAAAAAAAAAAAAA ????????????????????
+test 11 G 20 GGGGGGGGGGGGGGGGGGGG ????????????????????
+test 12 C 20 CCCCCCCCCCCCCCCCCCCC ????????????????????
+test 13 C 20 CCCCCCCCCCCCCCCCCCCC ????????????????????
+test 14 G 20 GGGGGGGGGGGGGGGGGGGG ????????????????????
+test 15 G 20 GGGGGGGGGGCCCCCCCCCC ????????????????????
+test 16 A 20 AAAAAAAAAAAAAAAAAAAA ????????????????????
+test 17 C 20 CCCCCCCCCCCCCCCCCCCC ????????????????????
+test 18 T 20 TTTTTTTTTTTTTTTTTTTT ????????????????????
+test 19 T 20 TTTTTTTTTTTTTTTTTTTT ????????????????????
+test 20 C 20 CCCCCCCCCCCCCCCCCCCC ????????????????????
=====================================
data/test_strand.gff
=====================================
@@ -0,0 +1,3 @@
+## GFF3 file format
+test Genbank CDS 1 9 . + . ID=test1;Note=positiveStrand;gene=A;
+test Genbank CDS 11 20 . - . ID=test2;Note=negativeStrand;gene=B;
=====================================
data/test_unpaired.sorted.bam
=====================================
Binary files /dev/null and b/data/test_unpaired.sorted.bam differ
=====================================
data/test_unpaired.sorted.bam.bai
=====================================
Binary files /dev/null and b/data/test_unpaired.sorted.bam.bai differ
=====================================
debian/changelog
=====================================
@@ -1,3 +1,13 @@
+ivar (1.4.3+dfsg-1) unstable; urgency=medium
+
+ * New upstream version 1.4.3+dfsg
+ * forward-build-options.patch: refresh.
+ * gcc-13-2.patch: refresh.
+ * d/rules: move permission fixup to dh_fixperms step.
+ * d/control: declare compliance to standards version 4.7.0.
+
+ -- Étienne Mollier <emollier at debian.org> Sun, 14 Jul 2024 11:16:30 +0200
+
ivar (1.4.2+dfsg-3) unstable; urgency=medium
* d/clean: new: clean build artifacts. (Closes: #1047333)
=====================================
debian/control
=====================================
@@ -6,7 +6,7 @@ Section: science
Priority: optional
Build-Depends: debhelper-compat (= 13),
libhts-dev
-Standards-Version: 4.6.2
+Standards-Version: 4.7.0
Vcs-Browser: https://salsa.debian.org/med-team/ivar
Vcs-Git: https://salsa.debian.org/med-team/ivar.git
Homepage: https://github.com/andersen-lab/ivar
=====================================
debian/patches/forward-build-options.patch
=====================================
@@ -19,8 +19,8 @@ This patch header follows DEP-3: http://dep.debian.net/deps/dep3/
@@ -1,6 +1,6 @@
LIBS = -lhts -lz -lpthread
--CXXFLAGS = -g -std=c++11 -Wall -Wextra -Werror
-+CXXFLAGS += -g -std=c++11 -Wall -Wextra -Werror
+-CXXFLAGS = -g -std=c++11 -Wall -Wextra -Werror -Wno-unused-but-set-variable
++CXXFLAGS += -g -std=c++11 -Wall -Wextra -Werror -Wno-unused-but-set-variable
- TESTS = check_primer_trim check_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_consensus_seq_id check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases check_isize_trim check_interval_tree check_amplicon_search check_consensus_min_insert_threshold
- check_PROGRAMS = check_primer_trim check_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_consensus_seq_id check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases check_isize_trim check_interval_tree check_amplicon_search check_consensus_min_insert_threshold
+ TESTS = check_primer_trim check_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_consensus_seq_id check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases check_isize_trim check_interval_tree check_amplicon_search check_consensus_min_insert_threshold check_quality_trim_unpaired check_strand_variants
+ check_PROGRAMS = check_primer_trim check_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_consensus_seq_id check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases check_isize_trim check_interval_tree check_amplicon_search check_consensus_min_insert_threshold check_quality_trim_unpaired check_strand_variants
=====================================
debian/patches/gcc-13-2.patch
=====================================
@@ -3,8 +3,8 @@ Last-Update: Fri, 04 Aug 2023 17:48:43 +0300
Author: Adrian Bunk <bunk at debian.org>
Bug-Debian: https://bugs.debian.org/1043025
---- ivar-1.4.2+dfsg.orig/src/Makefile.am
-+++ ivar-1.4.2+dfsg/src/Makefile.am
+--- ivar.orig/src/Makefile.am
++++ ivar/src/Makefile.am
@@ -1,6 +1,6 @@
LIBS = -lhts -lz -lpthread
@@ -13,13 +13,13 @@ Bug-Debian: https://bugs.debian.org/1043025
# this lists the binaries to produce, the (non-PHONY, binary) targets in
# the previous manual Makefile
---- ivar-1.4.2+dfsg.orig/tests/Makefile.am
-+++ ivar-1.4.2+dfsg/tests/Makefile.am
+--- ivar.orig/tests/Makefile.am
++++ ivar/tests/Makefile.am
@@ -1,6 +1,6 @@
LIBS = -lhts -lz -lpthread
--CXXFLAGS += -g -std=c++11 -Wall -Wextra -Werror
-+CXXFLAGS += -g -Wall -Wextra -Werror
+-CXXFLAGS += -g -std=c++11 -Wall -Wextra -Werror -Wno-unused-but-set-variable
++CXXFLAGS += -g -Wall -Wextra -Werror -Wno-unused-but-set-variable
- TESTS = check_primer_trim check_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_consensus_seq_id check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases check_isize_trim check_interval_tree check_amplicon_search check_consensus_min_insert_threshold
- check_PROGRAMS = check_primer_trim check_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_consensus_seq_id check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases check_isize_trim check_interval_tree check_amplicon_search check_consensus_min_insert_threshold
+ TESTS = check_primer_trim check_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_consensus_seq_id check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases check_isize_trim check_interval_tree check_amplicon_search check_consensus_min_insert_threshold check_quality_trim_unpaired check_strand_variants
+ check_PROGRAMS = check_primer_trim check_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_consensus_seq_id check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases check_isize_trim check_interval_tree check_amplicon_search check_consensus_min_insert_threshold check_quality_trim_unpaired check_strand_variants
=====================================
debian/rules
=====================================
@@ -16,7 +16,7 @@ ifeq (,$(filter nocheck,$(DEB_BUILD_OPTIONS)))
dh_auto_test --no-parallel
endif
-# FIXME: maybe use the clearer dh_fixperm target instead (but first locate the
-# file in the installation tree; could be a target for find…)
-execute_before_dh_auto_install:
- chmod -x $(CURDIR)/data/primer_pair_test/test_primer_pair.bed
+EXAMPLES_DIR = $(CURDIR)/debian/ivar/usr/share/doc/ivar/examples
+execute_before_dh_fixperms:
+ find -name test_primer_pair.bed
+ chmod -v -x $(EXAMPLES_DIR)/data/primer_pair_test/test_primer_pair.bed
=====================================
src/ivar.cpp
=====================================
@@ -17,7 +17,7 @@
#include "suffix_tree.h"
#include "trim_primer_quality.h"
-const std::string VERSION = "1.4.2";
+const std::string VERSION = "1.4.3";
struct args_t {
std::string bam; // -i
=====================================
src/parse_gff.cpp
=====================================
@@ -37,18 +37,18 @@ gff3_feature::gff3_feature(std::string line) {
}
ctr++;
}
- if (ctr < 9) std::cout << "GFF file is not in GFF3 file format!" << std::endl;
+ if (ctr < 9) std::cerr << "GFF file is not in GFF3 file format!" << std::endl;
line_stream.clear();
}
int gff3_feature::print() {
- std::cout << seqid << "\t" << source << "\t" << type << "\t" << start << "\t"
+ std::cerr << seqid << "\t" << source << "\t" << type << "\t" << start << "\t"
<< end << "\t" << score << "\t" << strand << "\t" << phase << "\t";
std::map<std::string, std::string>::iterator it;
for (it = attributes.begin(); it != attributes.end(); it++) {
- std::cout << it->first << ": " << it->second << "; ";
+ std::cerr << it->first << ": " << it->second << "; ";
}
- std::cout << std::endl;
+ std::cerr << std::endl;
return 0;
}
@@ -85,6 +85,8 @@ int gff3_feature::get_phase() { return phase; }
std::string gff3_feature::get_type() { return type; }
+char gff3_feature::get_strand() { return strand; }
+
int64_t gff3_feature::get_edit_position() {
int64_t edit_pos = -1;
std::map<std::string, std::string>::iterator it;
@@ -129,16 +131,22 @@ gff3::gff3(std::string path) {
int gff3::read_file(std::string path) {
std::ifstream fin = std::ifstream(path);
if (!fin) {
- std::cout << "GFF file does not exist at " << path << std::endl;
+ std::cerr << "GFF file does not exist at " << path << std::endl;
return -1;
}
std::string line;
while (std::getline(fin, line)) {
if (line[0] == '#') // Avoid comments in GFF file
continue;
- features.push_back(gff3_feature(line));
+ if(!line.empty())
+ features.push_back(gff3_feature(line));
+ }
+ if(!features.empty()){
+ this->is_empty = false;
+ } else {
+ std::cerr << "GFF file is empty!" << std::endl;
}
- this->is_empty = false;
+
return 0;
}
=====================================
src/ref_seq.cpp
=====================================
@@ -1,5 +1,26 @@
#include "ref_seq.h"
+// Complement base array
+// Source: https://github.com/samtools/samtools/blob/024e3d5ef0a2b0b7049f1fde6faebe8249988a06/faidx.c#L56
+const unsigned char comp_base[256] = {
+ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
+ 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
+ 32, '!', '"', '#', '$', '%', '&', '\'', '(', ')', '*', '+', ',', '-', '.', '/',
+ '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', ':', ';', '<', '=', '>', '?',
+ '@', 'T', 'V', 'G', 'H', 'E', 'F', 'C', 'D', 'I', 'J', 'M', 'L', 'K', 'N', 'O',
+ 'P', 'Q', 'Y', 'S', 'A', 'A', 'B', 'W', 'X', 'R', 'Z', '[', '\\',']', '^', '_',
+ '`', 't', 'v', 'g', 'h', 'e', 'f', 'c', 'd', 'i', 'j', 'm', 'l', 'k', 'n', 'o',
+ 'p', 'q', 'y', 's', 'a', 'a', 'b', 'w', 'x', 'r', 'z', '{', '|', '}', '~', 127,
+ 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
+ 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159,
+ 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175,
+ 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191,
+ 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207,
+ 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223,
+ 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239,
+ 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255,
+};
+
char ref_antd::get_base(int64_t pos, std::string region) { // 1-based position
int len;
char base = 0;
@@ -11,6 +32,13 @@ char ref_antd::get_base(int64_t pos, std::string region) { // 1-based position
return base;
}
+void ref_antd::reverse_complement_codon(char* codon) {
+ char temp = comp_base[(unsigned char)codon[2]];
+ codon[2] = comp_base[(unsigned char)codon[0]];
+ codon[0] = temp;
+ codon[1] = comp_base[(unsigned char)codon[1]];
+}
+
char *ref_antd::get_codon(int64_t pos, std::string region,
gff3_feature feature) {
int len;
@@ -19,7 +47,6 @@ char *ref_antd::get_codon(int64_t pos, std::string region,
std::string edit_sequence = feature.get_edit_sequence();
int64_t edit_sequence_size = edit_sequence.size();
char *codon = new char[3];
- int i;
int64_t edit_offset = 0;
if (pos > edit_pos + edit_sequence_size && edit_pos != -1) {
edit_offset =
@@ -27,12 +54,22 @@ char *ref_antd::get_codon(int64_t pos, std::string region,
? edit_sequence_size
: (pos - edit_pos); // Account for edits in position of insertion
}
- codon_start_pos =
- (feature.get_start() - 1) + feature.get_phase() +
- (((pos + edit_offset - (feature.get_start() + feature.get_phase()))) /
- 3) *
- 3;
- for (i = 0; i < 3; ++i) {
+ // codon_start_pos is w.r.t the reference sequence
+ if (feature.get_strand() == '-'){
+ codon_start_pos =
+ (feature.get_end() - 1) - feature.get_phase() -
+ (((feature.get_end() - feature.get_phase() - (pos + edit_offset))) /
+ 3) *
+ 3;
+ codon_start_pos -= 2; // This is to get to codon start from the 3' end and then take reverse complement
+ } else {
+ codon_start_pos =
+ (feature.get_start() - 1) + feature.get_phase() +
+ (((pos + edit_offset - (feature.get_start() + feature.get_phase()))) /
+ 3) *
+ 3;
+ }
+ for (int i = 0; i < 3; i++) {
if (codon_start_pos + i < (int32_t)feature.get_start() - 1 ||
codon_start_pos + i >
(int32_t)feature.get_end() -
@@ -53,6 +90,11 @@ char *ref_antd::get_codon(int64_t pos, std::string region,
codon[i] = *(seq + codon_start_pos + i - edit_offset);
}
}
+
+ if (feature.get_strand() == '-') {
+ reverse_complement_codon(codon);
+ }
+
free(seq);
return codon;
}
@@ -73,11 +115,22 @@ char *ref_antd::get_codon(int64_t pos, std::string region, gff3_feature feature,
? edit_sequence_size
: (pos - edit_pos); // Account for edits in position of insertion
}
- codon_start_pos =
- (feature.get_start() - 1) + feature.get_phase() +
- (((pos + edit_offset - (feature.get_start() + feature.get_phase()))) /
- 3) *
- 3;
+ // TODO: Remove code duplication with function above
+ // codon_start_pos is w.r.t the reference sequence
+ if (feature.get_strand() == '-'){
+ codon_start_pos =
+ (feature.get_end() - 1) - feature.get_phase() -
+ (((feature.get_end() - feature.get_phase() - (pos + edit_offset))) /
+ 3) *
+ 3;
+ codon_start_pos -= 2; // This is to get to codon start from the 3' end and then take reverse complement
+ } else {
+ codon_start_pos =
+ (feature.get_start() - 1) + feature.get_phase() +
+ (((pos + edit_offset - (feature.get_start() + feature.get_phase()))) /
+ 3) *
+ 3;
+ }
for (i = 0; i < 3; ++i) {
if (codon_start_pos + i < edit_pos - 1 ||
edit_pos == -1) { // If before edit or with no edit
@@ -104,6 +157,11 @@ char *ref_antd::get_codon(int64_t pos, std::string region, gff3_feature feature,
}
alt_pos += edit_offset;
codon[alt_pos - 1 - codon_start_pos] = alt;
+
+ if (feature.get_strand() == '-') {
+ reverse_complement_codon(codon);
+ }
+
free(seq);
return codon;
}
@@ -168,10 +226,17 @@ int ref_antd::codon_aa_stream(std::string region,
fout << codon2aa(alt_codon[0], alt_codon[1], alt_codon[2]) << "\t";
// adding amino acid position
- int64_t start = it->get_start();
- int64_t aa_pos = ((pos - start) / 3) + 1;
- fout << aa_pos;
- fout << std::endl;
+ // factor in translation direction
+ char strand = it->get_strand();
+ int64_t aa_pos;
+ if (strand == '-') {
+ int64_t end = it->get_end();
+ aa_pos = ((end - pos) / 3) + 1;
+ } else { // when strand is equal to '+', '?', or others
+ int64_t start = it->get_start();
+ aa_pos = ((pos - start) / 3) + 1;
+ }
+ fout << aa_pos << std::endl;
delete[] ref_codon;
delete[] alt_codon;
=====================================
src/ref_seq.h
=====================================
@@ -11,6 +11,7 @@
#define ref_seq
const char UNKNOWN_BASE = 'N';
+extern const unsigned char comp_base[256];
class ref_antd {
public:
@@ -18,6 +19,7 @@ class ref_antd {
ref_antd(std::string ref_path, std::string gff_path);
~ref_antd();
char get_base(int64_t pos, std::string region);
+ void reverse_complement_codon(char* codon);
int add_gff(std::string path);
int add_seq(std::string path);
int codon_aa_stream(std::string region, std::ostringstream &line_stream,
=====================================
src/trim_primer_quality.cpp
=====================================
@@ -89,7 +89,7 @@ cigar_ quality_trim(bam1_t *r, uint8_t qual_threshold, uint8_t sliding_window) {
*cigar = bam_get_cigar(r);
uint8_t *qual = bam_get_qual(r);
int32_t start_pos;
- if (((r->core.flag & BAM_FPAIRED) != 0) && bam_is_rev(r)) {
+ if (bam_is_rev(r)) {
reverse = true;
reverse_qual(qual, r->core.l_qseq);
}
@@ -642,7 +642,6 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out,
std::vector<primer> sorted_primers = insertionSort(primers, primers.size());
std::vector<bam1_t *>::iterator aln_itr = alns.begin();
-
// Iterate through reads
while (iterate_aln(aln_itr, alns, header, in, aln) >= 0) {
unmapped_flag = false;
@@ -665,7 +664,6 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out,
continue;
}
}
-
// isize is insert size
// l_qseq is the query length
isize_flag =
@@ -715,25 +713,25 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out,
// Update read's left-most coordinate
aln->core.pos += t.start_pos;
replace_cigar(aln, t.nlength, t.cigar);
+ free(t.cigar);
// Add count to primer
cit = std::find(primers.begin(), primers.end(), cand_primer);
if (cit != primers.end()) cit->add_read_count(1);
}
-
// reverse primer unpaired read
get_overlapping_primers(aln, primers, overlapping_primers, true);
-
if (overlapping_primers.size() > 0) {
primer_trimmed = true;
cand_primer = get_min_start(overlapping_primers);
t = primer_trim(aln, isize_flag, cand_primer.get_start() - 1, true);
replace_cigar(aln, t.nlength, t.cigar);
+ free(t.cigar);
// Add count to primer
cit = std::find(primers.begin(), primers.end(), cand_primer);
if (cit != primers.end()) cit->add_read_count(1);
}
t = quality_trim(aln, min_qual, sliding_window); // Quality Trimming
- if (bam_is_rev(aln)) // if reverse strand
+ if (bam_is_rev(aln)) // if reverse strand with reverse primers trimmed
aln->core.pos = t.start_pos;
condense_cigar(&t);
replace_cigar(aln, t.nlength, t.cigar);
=====================================
tests/Makefile.am
=====================================
@@ -1,12 +1,12 @@
LIBS = -lhts -lz -lpthread
-CXXFLAGS = -g -std=c++11 -Wall -Wextra -Werror
+CXXFLAGS = -g -std=c++11 -Wall -Wextra -Werror -Wno-unused-but-set-variable
-TESTS = check_primer_trim check_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_consensus_seq_id check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases check_isize_trim check_interval_tree check_amplicon_search check_consensus_min_insert_threshold
-check_PROGRAMS = check_primer_trim check_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_consensus_seq_id check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases check_isize_trim check_interval_tree check_amplicon_search check_consensus_min_insert_threshold
+TESTS = check_primer_trim check_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_consensus_seq_id check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases check_isize_trim check_interval_tree check_amplicon_search check_consensus_min_insert_threshold check_quality_trim_unpaired check_strand_variants
+check_PROGRAMS = check_primer_trim check_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_consensus_seq_id check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases check_isize_trim check_interval_tree check_amplicon_search check_consensus_min_insert_threshold check_quality_trim_unpaired check_strand_variants
check_primer_trim_SOURCES = test_primer_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp ../src/interval_tree.cpp
check_trim_SOURCES = test_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp ../src/interval_tree.cpp
-check_quality_trim_SOURCES = check_quality_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp ../src/interval_tree.cpp
+check_quality_trim_SOURCES = test_quality_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp ../src/interval_tree.cpp
check_consensus_SOURCES = test_call_consensus_from_plup.cpp ../src/call_consensus_pileup.cpp ../src/allele_functions.cpp
check_allele_depth_SOURCES = test_allele_depth.cpp ../src/allele_functions.cpp
check_consensus_threshold_SOURCES = test_consensus_threshold.cpp ../src/call_consensus_pileup.cpp ../src/allele_functions.cpp
@@ -23,3 +23,5 @@ check_primer_trim_edge_cases_SOURCES = test_primer_trim_edge_cases.cpp ../src/tr
check_isize_trim_SOURCES = test_isize_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp ../src/interval_tree.cpp
check_interval_tree_SOURCES = test_interval_tree.cpp ../src/primer_bed.cpp ../src/interval_tree.cpp
check_amplicon_search_SOURCES = test_amplicon_search.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp ../src/interval_tree.cpp
+check_quality_trim_unpaired_SOURCES = test_quality_trim_unpaired.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp ../src/interval_tree.cpp
+check_strand_variants_SOURCES = test_variants_reverse_strand.cpp ../src/call_variants.cpp ../src/allele_functions.cpp ../src/parse_gff.cpp ../src/ref_seq.cpp
=====================================
tests/check_quality_trim.cpp → tests/test_quality_trim.cpp
=====================================
=====================================
tests/test_quality_trim_unpaired.cpp
=====================================
@@ -0,0 +1,43 @@
+#include <iostream>
+
+#include "../src/trim_primer_quality.h"
+#include "htslib/sam.h"
+
+int main() {
+ int success = 0;
+ std::string bam = "../data/test_unpaired.sorted.bam";
+ std::string region_;
+ samFile *in = hts_open(bam.c_str(), "r");
+ hts_idx_t *idx = sam_index_load(in, bam.c_str());
+ bam_hdr_t *header = sam_hdr_read(in);
+ region_.assign(header->target_name[0]);
+ std::string temp(header->text);
+ hts_itr_t *iter = NULL;
+ iter = sam_itr_querys(idx, header, region_.c_str());
+ bam1_t *aln = bam_init1();
+ cigar_ t;
+ int lengths[] = {323, 756, 320, 227}, ctr = 0;
+ int start_pos_rev[] = {456, 23, 456, 549};
+ int qual_thresholds[] = {20, 6, 20, 23};
+ while (sam_itr_next(in, iter, aln) >= 0) {
+ t = quality_trim(aln, qual_thresholds[ctr], 4);
+ std::cout << bam_get_qname(aln) << std::endl;
+ std::cout << "POS: " << t.start_pos
+ << " Length: " << bam_cigar2rlen(t.nlength, t.cigar) << std::endl;
+ if (t.start_pos != start_pos_rev[ctr]) {
+ success = -1;
+ }
+ if (bam_cigar2rlen(t.nlength, t.cigar) != lengths[ctr]) {
+ success = -1;
+ }
+ free_cigar(t);
+ ctr++;
+ }
+
+ bam_destroy1(aln);
+ bam_itr_destroy(iter);
+ sam_hdr_destroy(header);
+ hts_idx_destroy(idx);
+ hts_close(in);
+ return success;
+}
=====================================
tests/test_variants_reverse_strand.cpp
=====================================
@@ -0,0 +1,45 @@
+#include <fstream>
+#include <iostream>
+#include <string>
+
+#include "../src/allele_functions.h"
+#include "../src/call_variants.h"
+
+int call_var_check_outfile(std::string prefix, uint8_t min_qual,
+ uint8_t min_depth, double min_threshold,
+ std::string out[], int len) {
+ std::string path = "../data/test.strand.pileup";
+ std::ifstream mplp(path);
+ call_variants_from_plup(mplp, prefix, min_qual, min_threshold, min_depth,
+ "../data/db/test_ref.fa", "../data/test_strand.gff");
+ std::ifstream outFile(prefix + ".tsv");
+ std::string l;
+ getline(outFile, l); // Ignore first line
+ int comp = 0, ctr = 0;
+
+ while (ctr < len) {
+ getline(outFile, l);
+ std::cout << l << std::endl;
+ std::cout << out[ctr] << " -> CORRECT" << std::endl;
+ comp += l.compare(out[ctr]);
+ std::cout << l.compare(out[ctr]) << "\n";
+ ctr++;
+ }
+ return comp;
+}
+
+int main() {
+ int num_success = 0;
+ // Quality threshold 20. Frequency threshold: 0.03. Total_DP = 3. Indel passes
+ // filters with total_depth 4. Has two lines.
+ std::string expected_out[2] = {
+ "test\t5\tC\tT\t10\t0\t30\t10\t0\t30\t0.5\t20\t0.000290613\tTRUE\tA:test1\tGCT\tA\tGTT\tV\t2",
+ "test\t15\tG\tC\t10\t0\t30\t10\t0\t30\t0.5\t20\t0.000290613\tTRUE\tB:test2\tGTC\tV\tGTG\tV\t2"
+ };
+ num_success =
+ call_var_check_outfile("../data/test_strand", 20, 0, 0.03, expected_out, 2);
+ std::cout << num_success << std::endl;
+
+ if (num_success == 0) return 0;
+ return -1;
+}
View it on GitLab: https://salsa.debian.org/med-team/ivar/-/compare/7f9e51e273e7c7dfc476cdb9aa804e7d00359bee...de6dc239157a5d6d41c08cbebbb63eba889669ed
--
This project does not include diff previews in email notifications.
View it on GitLab: https://salsa.debian.org/med-team/ivar/-/compare/7f9e51e273e7c7dfc476cdb9aa804e7d00359bee...de6dc239157a5d6d41c08cbebbb63eba889669ed
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20240714/16ef3430/attachment-0001.htm>
More information about the debian-med-commit
mailing list