[med-svn] [Git][med-team/ivar][upstream] New upstream version 1.2.3+dfsg
Andreas Tille
gitlab at salsa.debian.org
Wed Oct 14 08:24:51 BST 2020
Andreas Tille pushed to branch upstream at Debian Med / ivar
Commits:
701a9ea3 by Andreas Tille at 2020-10-14T09:15:32+02:00
New upstream version 1.2.3+dfsg
- - - - -
29 changed files:
- configure.ac
- data/pair_information.tsv
- docs/MANUAL.md
- docs/html/index.html
- docs/html/manualpage.html
- docs/latex/index.tex
- docs/latex/manualpage.tex
- docs/man/man3/cookbookpage.3
- docs/man/man3/installpage.3
- docs/man/man3/manualpage.3
- src/call_consensus_pileup.cpp
- src/call_consensus_pileup.h
- src/call_variants.cpp
- src/get_masked_amplicons.cpp
- src/ivar.cpp
- src/primer_bed.cpp
- src/ref_seq.cpp
- src/ref_seq.h
- src/trim_primer_quality.cpp
- src/trim_primer_quality.h
- tests/Makefile.am
- tests/check_quality_trim.cpp
- tests/test_consensus_min_depth.cpp
- + tests/test_consensus_seq_id.cpp
- tests/test_getmasked.cpp
- tests/test_primer_trim.cpp
- tests/test_primer_trim_edge_cases.cpp
- + tests/test_trim.cpp
- tests/test_unpaired_trim.cpp
Changes:
=====================================
configure.ac
=====================================
@@ -2,7 +2,7 @@
# Process this file with autoconf to produce a configure script.
AC_PREREQ([2.63])
-AC_INIT([ivar], [1.2.2], [gkarthik at scripps.edu])
+AC_INIT([ivar], [1.2.3], [gkarthik at scripps.edu])
AM_INIT_AUTOMAKE([foreign subdir-objects])
AC_CONFIG_HEADERS([config.h])
=====================================
data/pair_information.tsv
=====================================
@@ -1,3 +1,3 @@
WNV_400_1_LEFT WNV_400_1_RIGHT
WNV_400_2_LEFT WNV_400_2_RIGHT
-WNV_400_3_LEFT WNV_400_3_RIGHT
\ No newline at end of file
+WNV_400_3_LEFT WNV_400_3_RIGHT
=====================================
docs/MANUAL.md
=====================================
@@ -27,7 +27,9 @@ Trim primer sequences with iVar
iVar uses primer positions supplied in a BED file to soft clip primer sequences from an aligned and sorted BAM file. Following this, the reads are trimmed based on a quality threshold(Default: 20). To do the quality trimming, iVar uses a sliding window approach(Default: 4). The windows slides from the 5' end to the 3' end and if at any point the average base quality in the window falls below the threshold, the remaining read is soft clipped. If after trimming, the length of the read is greater than the minimum length specified(Default: 30), the read is written to the new trimmed BAM file.
-To sort and index an aligned BAM file, the following command can be used,
+Please note that the strand is taken into account while doing the trimming so forward primers are trimmed only from forward strand and reverse primers are trimmed from reverse strand.
+
+To sort and index an aligned BAM file (OPTIONAL, if index is not present iVar will create one), the following command can be used,
```
# Input file - test.bam
@@ -115,7 +117,7 @@ Output Options Description
Example Usage:
```
-samtools mpileup -A -d 600000 -B -Q 0 test.trimmed.bam | ivar variants -p test -q 20 -t 0.03 -r test_reference.fa -g test.gff
+samtools mpileup -aa -A -d 600000 -B -Q 0 test.trimmed.bam | ivar variants -p test -q 20 -t 0.03 -r test_reference.fa -g test.gff
```
The command above will generate a test.tsv file.
@@ -262,12 +264,14 @@ If there are two nucleotides at the same frequency, both nucleotides are used to
The output of the command is a fasta file with the consensus sequence and a .txt file with the average quality of every base used to generate the consensus at each position. *For insertions, the quality is set to be the minimum quality threshold since mpileup doesn't give the quality of bases in insertions.*
Command:
+
```
+
ivar consensus
Usage: samtools mpileup -aa -A -d 0 -Q 0 <input.bam> | ivar consensus -p <prefix>
-Note : samtools mpileup output must be piped into `ivar consensus`
+Note : samtools mpileup output must be piped into ivar consensus
Input Options Description
-q Minimum quality score threshold to count base (Default: 20)
=====================================
docs/html/index.html
=====================================
@@ -74,7 +74,7 @@ $(function() {
<li><a class="el" href="manualpage.html">Manual</a></li>
<li><a class="el" href="cookbookpage.html">Cookbook</a></li>
</ul>
-<h1><a class="anchor" id="autotoc_md21"></a>
+<h1><a class="anchor" id="autotoc_md23"></a>
Publication</h1>
<p><a href="https://doi.org/10.1186/s13059-018-1618-7">An amplicon-based sequencing framework for accurately measuring intrahost virus diversity using PrimalSeq and iVar</a></p>
<p><em>Genome Biology</em> 2019 <b>20</b>:8</p>
=====================================
docs/html/manualpage.html
=====================================
@@ -75,7 +75,9 @@ $(function() {
<li class="level2"><a href="#autotoc_md17">Call variants with iVar</a></li>
<li class="level2"><a href="#autotoc_md18">Filter variants across replicates with iVar</a></li>
<li class="level2"><a href="#autotoc_md19">Generate a consensus sequences from an aligned BAM file</a></li>
-<li class="level2"><a href="#autotoc_md20">(Experimental) trimadapter</a></li>
+<li class="level2"><a href="#autotoc_md20">Get primers with mismatches to the reference sequence</a></li>
+<li class="level2"><a href="#autotoc_md21">Remove reads associated with mismatched primer indices</a></li>
+<li class="level2"><a href="#autotoc_md22">(Experimental) trimadapter</a></li>
</ul>
</li>
</ul>
@@ -108,7 +110,8 @@ Description of commands</h1>
<h2><a class="anchor" id="autotoc_md16"></a>
Trim primer sequences with iVar</h2>
<p>iVar uses primer positions supplied in a BED file to soft clip primer sequences from an aligned and sorted BAM file. Following this, the reads are trimmed based on a quality threshold(Default: 20). To do the quality trimming, iVar uses a sliding window approach(Default: 4). The windows slides from the 5' end to the 3' end and if at any point the average base quality in the window falls below the threshold, the remaining read is soft clipped. If after trimming, the length of the read is greater than the minimum length specified(Default: 30), the read is written to the new trimmed BAM file.</p>
-<p>To sort and index an aligned BAM file, the following command can be used,</p>
+<p>Please note that the strand is taken into account while doing the trimming so forward primers are trimmed only from forward strand and reverse primers are trimmed from reverse strand.</p>
+<p>To sort and index an aligned BAM file (OPTIONAL, if index is not present iVar will create one), the following command can be used,</p>
<div class="fragment"><div class="line"># Input file - test.bam</div><div class="line">samtools sort -o test.sorted.bam test.bam && samtools index test.sorted.bam.</div></div><!-- fragment --><p><b>Note</b>: All the trimming in iVar is done by soft-clipping reads in an aligned BAM file. This information is lost if reads are extracted in fastq or fasta format from the trimmed BAM file.</p>
<p>Command: </p><div class="fragment"><div class="line">ivar trim</div><div class="line"></div><div class="line">Usage: ivar trim -i <input.bam> -b <primers.bed> -p <prefix> [-m <min-length>] [-q <min-quality>] [-s <sliding-window-width>]</div><div class="line"></div><div class="line">Input Options Description</div><div class="line"> -i (Required) Sorted bam file, with aligned reads, to trim primers and quality</div><div class="line"> -b (Required) BED file with primer sequences and positions</div><div class="line"> -m Minimum length of read to retain after trimming (Default: 30)</div><div class="line"> -q Minimum quality threshold for sliding window to pass (Default: 20)</div><div class="line"> -s Width of sliding window (Default: 4)</div><div class="line"> -e Include reads with no primers. By default, reads with no primers are excluded</div><div class="line"></div><div class="line">Output Options Description</div><div class="line"> -p (Required) Prefix for the output BAM file</div></div><!-- fragment --><p>Example Usage: </p><div class="fragment"><div class="line">ivar trim -b test_primers.bed -p test.trimmed -i test.bam -q 15 -m 50 -s 4</div></div><!-- fragment --><p>The command above will produce a trimmed BAM file test.trimmed.bam after trimming the aligned reads in test.bam using the primer positions specified in test_primers.bed and a minimum quality threshold of <b>15</b>, minimum read length of <b>50</b> and a sliding window of <b>4</b>.</p>
<p>Example BED file</p>
@@ -120,7 +123,7 @@ Call variants with iVar</h2>
<h4>Account for RNA editing through polymerase slippage</h4>
<p>Some RNA viruses such as Ebola virus, might have polymerase slippage causing the insertion of a couple of nucleotides. More details can be found <a href="https://viralzone.expasy.org/857?outline=all_by_protein">here</a>. iVar can account for this editing and identify the correct open reading frames. The user will have to specify two additional parameters, <b>EditPosition</b>: Position at which edit occurs and <b>EditSequence</b>: The sequence tht is inserted at the given positon, in the "attributes" column of the GFF file to account for this. A test example is given below,</p>
<div class="fragment"><div class="line">test Genbank CDS 2 292 . + . ID=id-testedit1;Note=PinkFloyd;EditPosition=100;EditSequence=A</div><div class="line">test Genbank CDS 2 292 . + . ID=id-testedit2;Note=AnotherBrickInTheWall;EditPosition=102;EditSequence=AA</div></div><!-- fragment --><p>If a certain base is present in multiple CDSs, iVar will add a new row for each CDS frame and distinguish the rows by adding the ID (specified in attributes of GFF) of the GFF feature used for the translation. This is shown for position 42 in the example output below. There are two rows with two different GFF features: id-test3 and id-test4.</p>
-<p>Command: </p><div class="fragment"><div class="line">Usage: samtools mpileup -aa -A -d 0 -B -Q 0 --reference [<reference-fasta] <input.bam> | ivar variants -p <prefix> [-q <min-quality>] [-t <min-frequency-threshold>] [-m <minimum depth>] [-r <reference-fasta>] [-g GFF file]</div><div class="line"></div><div class="line">Note : samtools mpileup output must be piped into ivar variants</div><div class="line"></div><div class="line">Input Options Description</div><div class="line"> -q Minimum quality score threshold to count base (Default: 20)</div><div class="line"> -t Minimum frequency threshold(0 - 1) to call variants (Default: 0.03)</div><div class="line"> -m Minimum read depth to call variants (Default: 0)</div><div class="line"> -r Reference file used for alignment. This is used to translate the nucleotide sequences and identify intra host single nucleotide variants</div><div class="line"> -g A GFF file in the GFF3 format can be supplied to specify coordinates of open reading frames (ORFs). In absence of GFF file, amino acid translation will not be done.</div><div class="line"></div><div class="line">Output Options Description</div><div class="line"> -p (Required) Prefix for the output tsv variant file</div></div><!-- fragment --><p>Example Usage: </p><div class="fragment"><div class="line">samtools mpileup -A -d 600000 -B -Q 0 test.trimmed.bam | ivar variants -p test -q 20 -t 0.03 -r test_reference.fa -g test.gff</div></div><!-- fragment --><p>The command above will generate a test.tsv file.</p>
+<p>Command: </p><div class="fragment"><div class="line">Usage: samtools mpileup -aa -A -d 0 -B -Q 0 --reference [<reference-fasta] <input.bam> | ivar variants -p <prefix> [-q <min-quality>] [-t <min-frequency-threshold>] [-m <minimum depth>] [-r <reference-fasta>] [-g GFF file]</div><div class="line"></div><div class="line">Note : samtools mpileup output must be piped into ivar variants</div><div class="line"></div><div class="line">Input Options Description</div><div class="line"> -q Minimum quality score threshold to count base (Default: 20)</div><div class="line"> -t Minimum frequency threshold(0 - 1) to call variants (Default: 0.03)</div><div class="line"> -m Minimum read depth to call variants (Default: 0)</div><div class="line"> -r Reference file used for alignment. This is used to translate the nucleotide sequences and identify intra host single nucleotide variants</div><div class="line"> -g A GFF file in the GFF3 format can be supplied to specify coordinates of open reading frames (ORFs). In absence of GFF file, amino acid translation will not be done.</div><div class="line"></div><div class="line">Output Options Description</div><div class="line"> -p (Required) Prefix for the output tsv variant file</div></div><!-- fragment --><p>Example Usage: </p><div class="fragment"><div class="line">samtools mpileup -aa -A -d 600000 -B -Q 0 test.trimmed.bam | ivar variants -p test -q 20 -t 0.03 -r test_reference.fa -g test.gff</div></div><!-- fragment --><p>The command above will generate a test.tsv file.</p>
<p>Example of output .tsv file.</p>
<div class="fragment"><div class="line">REGION POS REF ALT REF_DP REF_RV REF_QUAL ALT_DP ALT_RV ALT_QUAL ALT_FREQ TOTAL_DP PVAL PASS GFF_FEATURE REF_CODON REF_AA ALT_CODON ALT_AA</div><div class="line">test 42 G T 0 0 0 1 0 49 1 1 1 FALSE id-test3 AGG R ATG M</div><div class="line">test 42 G T 0 0 0 1 0 49 1 1 1 FALSE id-test4 CAG Q CAT H</div><div class="line">test 320 A T 1 1 35 1 1 46 0.5 2 0.666667 FALSE NA NA NA NA NA</div><div class="line">test 365 A T 0 0 0 1 1 27 1 1 1 FALSE NA NA NA NA NA</div></div><!-- fragment --><p>Description</p>
<table class="markdownTable">
@@ -255,33 +258,19 @@ Generate a consensus sequences from an aligned BAM file</h2>
<td class="markdownTableBodyLeft">1 </td><td class="markdownTableBodyNone">D(A or T or G) </td></tr>
</table>
<p>The output of the command is a fasta file with the consensus sequence and a .txt file with the average quality of every base used to generate the consensus at each position. <em>For insertions, the quality is set to be the minimum quality threshold since mpileup doesn't give the quality of bases in insertions.</em></p>
-<p>Command: ``` ivar consensus</p>
-<p>Usage: samtools mpileup -aa -A -d 0 -Q 0 <input.bam> | ivar consensus -p <prefix></p>
-<p>Note : samtools mpileup output must be piped into <code>ivar consensus</code></p>
-<p>Input Options Description -q Minimum quality score threshold to count base (Default: 20) -t Minimum frequency threshold(0 - 1) to call consensus. (Default: 0) </p><table class="markdownTable">
-<tr class="markdownTableHead">
-<th class="markdownTableHeadNone">Frequently used thresholds </th><th class="markdownTableHeadNone">Description </th></tr>
-<tr class="markdownTableBody" class="markdownTableRowOdd">
-<td class="markdownTableBodyNone">0 </td><td class="markdownTableBodyNone">Majority or most common base </td></tr>
-<tr class="markdownTableBody" class="markdownTableRowEven">
-<td class="markdownTableBodyNone">0.2 </td><td class="markdownTableBodyNone">Bases that make up atleast 20% of the depth at a position </td></tr>
-<tr class="markdownTableBody" class="markdownTableRowOdd">
-<td class="markdownTableBodyNone">0.5 </td><td class="markdownTableBodyNone">Strict or bases that make up atleast 50% of the depth at a position </td></tr>
-<tr class="markdownTableBody" class="markdownTableRowEven">
-<td class="markdownTableBodyNone">0.9 </td><td class="markdownTableBodyNone">Strict or bases that make up atleast 90% of the depth at a position </td></tr>
-<tr class="markdownTableBody" class="markdownTableRowOdd">
-<td class="markdownTableBodyNone">1 </td><td class="markdownTableBodyNone">Identical or bases that make up 100% of the depth at a position. Will have highest ambiguities </td></tr>
-</table>
-<p>-m Minimum depth to call consensus(Default: 10) -k If '-k' flag is added, regions with depth less than minimum depth will not be added to the consensus sequence. Using '-k' will override any option specified using -n -n (N/-) Character to print in regions with less than minimum coverage(Default: N)</p>
-<p>Output Options Description -p (Required) Prefix for the output fasta file and quality file </p><div class="fragment"><div class="line">Example Usage:</div></div><!-- fragment --><p> samtools mpileup -d 1000 -A -Q 0 test.bam | ivar consensus -p test -q 20 -t 0 </p><div class="fragment"><div class="line">The command above will produce a test.fa fasta file with the consensus sequence and a test.qual.txt with the average quality of each base in the consensus sequence.</div><div class="line"></div><div class="line">Get primers with mismatches to the reference sequence</div><div class="line">----</div><div class="line"></div><div class="line">iVar uses a .tsv file with variants to get the zero based indices(based on the BED file) of mismatched primers. This command requires another .tsv file with each line containing the left and right primer names separated by a tab. This is used to get both the primers for an amplicon with a single mismatched primer. The output is a text file with the zero based primer indices delimited by a space. The output is written to a a text file using the prefix provided.</div><div class="line"></div><div class="line">Command:</div></div><!-- fragment --><p> ivar getmasked Usage: ivar getmasked -i <input-filtered.tsv> -b <primers.bed> -f <primer_pairs.tsv> -p <prefix> Note: This step is used only for amplicon-based sequencing.</p>
-<p>Input Options Description -i (Required) Input filtered variants tsv generated from 'ivar filtervariants' -b (Required) BED file with primer sequences and positions -f (Required) Primer pair information file containing left and right primer names for the same amplicon separated by a tab Output Options Description -p (Required) Prefix for the output text file</p>
-<div class="fragment"><div class="line">Example BED file</div></div><!-- fragment --><p> Puerto 28 52 400_1_out_L 60 + Puerto 482 504 400_1_out_R 60 - Puerto 359 381 400_2_out_L 60 + Puerto 796 818 400_2_out_R 60 - Puerto 658 680 400_3_out_L* 60 + Puerto 1054 1076 400_3_out_R* 60 -</p><div class="fragment"><div class="line">Example primer pair information file</div></div><!-- fragment --><p> 400_1_out_L 400_1_out_R 400_2_out_L 400_2_out_R 400_3_out_L 400_3_out_R</p><div class="fragment"><div class="line">Example Usage:</div></div><!-- fragment --><p> ivar getmasked -i test.filtered.tsv -b primers.bed -f pair_information.tsv -p test.masked.txt </p><div class="fragment"><div class="line">The command above produces an output file - test.masked.txt.</div><div class="line"></div><div class="line">Example Output:</div></div><!-- fragment --><p> 1 2 7 8 </p><div class="fragment"><div class="line">Remove reads associated with mismatched primer indices</div><div class="line">----</div><div class="line"></div><div class="line">This command accepts an aligned and sorted BAM file trimmed using `ivar trim` and removes the reads corresponding to the supplied primer indices, which is the output of `ivar getmasked` command. Under the hood, `ivar trim` adds the zero based primer index(based on the BED file) to the BAM auxillary data for every read. Hence, ivar removereads will only work on BAM files that have been trimmed using `ivar trim`.</div><div class="line"></div><div class="line">Command:</div></div><!-- fragment --><p> ivar removereads</p>
-<p>Usage: ivar removereads -i <input.trimmed.bam> -p <prefix> -t <text-file-with-primer-indices> Note: This step is used only for amplicon-based sequencing.</p>
-<p>Input Options Description -i (Required) Input BAM file trimmed with ivar trim. Must be sorted and indexed, which can be done using sort_index_bam.sh -t (Required) Text file with primer indices separated by spaces. This is the output of getmasked command.</p>
-<p>Output Options Description -p (Required) Prefix for the output filtered BAM file</p>
-<div class="fragment"><div class="line">Example Usage:</div></div><!-- fragment --><p> ivar trim -i test.bam -p test.trimmed ivar removereads -i test.trimmed.bam -p test.trimmed.masked.bam -t test.masked.txt ```</p>
-<p>The <code>ivar trim</code> command above trims test.bam and produced test.trimmed.bam with the primer indice data added. The <code>ivar removereads</code> command produces an output file - test.trimmed.masked.bam after removing all the reads corresponding to primer indices - 1,2,7 and 8.</p>
+<p>Command:</p>
+<div class="fragment"><div class="line">ivar consensus</div><div class="line"></div><div class="line">Usage: samtools mpileup -aa -A -d 0 -Q 0 <input.bam> | ivar consensus -p <prefix> </div><div class="line"></div><div class="line">Note : samtools mpileup output must be piped into ivar consensus</div><div class="line"></div><div class="line">Input Options Description</div><div class="line"> -q Minimum quality score threshold to count base (Default: 20)</div><div class="line"> -t Minimum frequency threshold(0 - 1) to call consensus. (Default: 0)</div><div class="line"> Frequently used thresholds | Description</div><div class="line"> ---------------------------|------------</div><div class="line"> 0 | Majority or most common base</div><div class="line"> 0.2 | Bases that make up atleast 20% of the depth at a position</div><div class="line"> 0.5 | Strict or bases that make up atleast 50% of the depth at a position</div><div class="line"> 0.9 | Strict or bases that make up atleast 90% of the depth at a position</div><div class="line"> 1 | Identical or bases that make up 100% of the depth at a position. Will have highest ambiguities</div><div class="line"> -m Minimum depth to call consensus(Default: 10)</div><div class="line"> -k If '-k' flag is added, regions with depth less than minimum depth will not be added to the consensus sequence. Using '-k' will override any option specified using -n </div><div class="line"> -n (N/-) Character to print in regions with less than minimum coverage(Default: N)</div><div class="line"></div><div class="line">Output Options Description</div><div class="line"> -p (Required) Prefix for the output fasta file and quality file</div></div><!-- fragment --><p>Example Usage: </p><div class="fragment"><div class="line">samtools mpileup -d 1000 -A -Q 0 test.bam | ivar consensus -p test -q 20 -t 0</div></div><!-- fragment --><p>The command above will produce a test.fa fasta file with the consensus sequence and a test.qual.txt with the average quality of each base in the consensus sequence.</p>
<h2><a class="anchor" id="autotoc_md20"></a>
+Get primers with mismatches to the reference sequence</h2>
+<p>iVar uses a .tsv file with variants to get the zero based indices(based on the BED file) of mismatched primers. This command requires another .tsv file with each line containing the left and right primer names separated by a tab. This is used to get both the primers for an amplicon with a single mismatched primer. The output is a text file with the zero based primer indices delimited by a space. The output is written to a a text file using the prefix provided.</p>
+<p>Command: </p><div class="fragment"><div class="line">ivar getmasked</div><div class="line">Usage: ivar getmasked -i <input-filtered.tsv> -b <primers.bed> -f <primer_pairs.tsv> -p <prefix></div><div class="line">Note: This step is used only for amplicon-based sequencing.</div><div class="line"></div><div class="line">Input Options Description</div><div class="line"> -i (Required) Input filtered variants tsv generated from 'ivar filtervariants'</div><div class="line"> -b (Required) BED file with primer sequences and positions</div><div class="line"> -f (Required) Primer pair information file containing left and right primer names for the same amplicon separated by a tab</div><div class="line">Output Options Description</div><div class="line"> -p (Required) Prefix for the output text file</div></div><!-- fragment --><p>Example BED file</p>
+<div class="fragment"><div class="line">Puerto 28 52 400_1_out_L 60 +</div><div class="line">Puerto 482 504 400_1_out_R 60 -</div><div class="line">Puerto 359 381 400_2_out_L 60 +</div><div class="line">Puerto 796 818 400_2_out_R 60 -</div><div class="line">Puerto 658 680 400_3_out_L* 60 +</div><div class="line">Puerto 1054 1076 400_3_out_R* 60 -</div><div class="line">.</div><div class="line">.</div><div class="line">.</div><div class="line">.</div></div><!-- fragment --><p>Example primer pair information file </p><div class="fragment"><div class="line">400_1_out_L 400_1_out_R</div><div class="line">400_2_out_L 400_2_out_R</div><div class="line">400_3_out_L 400_3_out_R</div><div class="line">.</div><div class="line">.</div><div class="line">.</div><div class="line">.</div></div><!-- fragment --><p>Example Usage: </p><div class="fragment"><div class="line">ivar getmasked -i test.filtered.tsv -b primers.bed -f pair_information.tsv -p test.masked.txt</div></div><!-- fragment --><p>The command above produces an output file - test.masked.txt.</p>
+<p>Example Output:</p>
+<div class="fragment"><div class="line">1 2 7 8</div></div><!-- fragment --><h2><a class="anchor" id="autotoc_md21"></a>
+Remove reads associated with mismatched primer indices</h2>
+<p>This command accepts an aligned and sorted BAM file trimmed using <code>ivar trim</code> and removes the reads corresponding to the supplied primer indices, which is the output of <code>ivar getmasked</code> command. Under the hood, <code>ivar trim</code> adds the zero based primer index(based on the BED file) to the BAM auxillary data for every read. Hence, ivar removereads will only work on BAM files that have been trimmed using <code>ivar trim</code>.</p>
+<p>Command: </p><div class="fragment"><div class="line">ivar removereads</div><div class="line"></div><div class="line">Usage: ivar removereads -i <input.trimmed.bam> -p <prefix> -t <text-file-with-primer-indices></div><div class="line">Note: This step is used only for amplicon-based sequencing.</div><div class="line"></div><div class="line">Input Options Description</div><div class="line"> -i (Required) Input BAM file trimmed with ivar trim. Must be sorted and indexed, which can be done using sort_index_bam.sh</div><div class="line"> -t (Required) Text file with primer indices separated by spaces. This is the output of getmasked command.</div><div class="line"></div><div class="line">Output Options Description</div><div class="line"> -p (Required) Prefix for the output filtered BAM file</div></div><!-- fragment --><p>Example Usage: </p><div class="fragment"><div class="line">ivar trim -i test.bam -p test.trimmed</div><div class="line">ivar removereads -i test.trimmed.bam -p test.trimmed.masked.bam -t test.masked.txt</div></div><!-- fragment --><p>The <code>ivar trim</code> command above trims test.bam and produced test.trimmed.bam with the primer indice data added. The <code>ivar removereads</code> command produces an output file - test.trimmed.masked.bam after removing all the reads corresponding to primer indices - 1,2,7 and 8.</p>
+<h2><a class="anchor" id="autotoc_md22"></a>
(Experimental) trimadapter</h2>
<p><b>Note: This feature is under active development and not completely validated yet.</b></p>
<p>trimadapter in iVar can be used to trim adapter sequences from fastq files using a supplied fasta file. </p>
=====================================
docs/latex/index.tex
=====================================
@@ -3,7 +3,7 @@
\item \mbox{\hyperlink{installpage}{Installation}}
\item \mbox{\hyperlink{manualpage}{Manual}}
\item \mbox{\hyperlink{cookbookpage}{Cookbook}}
-\end{DoxyItemize}\hypertarget{index_autotoc_md21}{}\section{Publication}\label{index_autotoc_md21}
+\end{DoxyItemize}\hypertarget{index_autotoc_md23}{}\section{Publication}\label{index_autotoc_md23}
\href{https://doi.org/10.1186/s13059-018-1618-7}{\tt An amplicon-\/based sequencing framework for accurately measuring intrahost virus diversity using Primal\+Seq and i\+Var}
{\itshape Genome Biology} 2019 {\bfseries 20}\+:8
=====================================
docs/latex/manualpage.tex
=====================================
@@ -24,7 +24,9 @@ To view detailed usage for each command type {\ttfamily ivar $<$command$>$} Note
\hypertarget{manualpage_autotoc_md16}{}\subsection{Trim primer sequences with i\+Var}\label{manualpage_autotoc_md16}
i\+Var uses primer positions supplied in a B\+ED file to soft clip primer sequences from an aligned and sorted B\+AM file. Following this, the reads are trimmed based on a quality threshold(\+Default\+: 20). To do the quality trimming, i\+Var uses a sliding window approach(\+Default\+: 4). The windows slides from the 5\textquotesingle{} end to the 3\textquotesingle{} end and if at any point the average base quality in the window falls below the threshold, the remaining read is soft clipped. If after trimming, the length of the read is greater than the minimum length specified(\+Default\+: 30), the read is written to the new trimmed B\+AM file.
-To sort and index an aligned B\+AM file, the following command can be used,
+Please note that the strand is taken into account while doing the trimming so forward primers are trimmed only from forward strand and reverse primers are trimmed from reverse strand.
+
+To sort and index an aligned B\+AM file (O\+P\+T\+I\+O\+N\+AL, if index is not present i\+Var will create one), the following command can be used,
\begin{DoxyCode}
@@ -123,7 +125,7 @@ Output Options Description
Example Usage\+:
\begin{DoxyCode}
-samtools mpileup -A -d 600000 -B -Q 0 test.trimmed.bam | ivar variants -p test -q 20 -t 0.03 -r
+samtools mpileup -aa -A -d 600000 -B -Q 0 test.trimmed.bam | ivar variants -p test -q 20 -t 0.03 -r
test\_reference.fa -g test.gff
\end{DoxyCode}
@@ -326,99 +328,136 @@ If there are two nucleotides at the same frequency, both nucleotides are used to
The output of the command is a fasta file with the consensus sequence and a .txt file with the average quality of every base used to generate the consensus at each position. {\itshape For insertions, the quality is set to be the minimum quality threshold since mpileup doesn\textquotesingle{}t give the quality of bases in insertions.}
-Command\+: \`{}\`{}\`{} ivar consensus
+Command\+:
-Usage\+: samtools mpileup -\/aa -\/A -\/d 0 -\/Q 0 $<$input.\+bam$>$ $\vert$ ivar consensus -\/p $<$prefix$>$
-Note \+: samtools mpileup output must be piped into {\ttfamily ivar consensus}
+\begin{DoxyCode}
+ivar consensus
-Input Options Description -\/q Minimum quality score threshold to count base (Default\+: 20) -\/t Minimum frequency threshold(0 -\/ 1) to call consensus. (Default\+: 0) \tabulinesep=1mm
-\begin{longtabu} spread 0pt [c]{*{2}{|X[-1]}|}
-\hline
-\rowcolor{\tableheadbgcolor}\textbf{ Frequently used thresholds }&\textbf{ Description }\\\cline{1-2}
-\endfirsthead
-\hline
-\endfoot
-\hline
-\rowcolor{\tableheadbgcolor}\textbf{ Frequently used thresholds }&\textbf{ Description }\\\cline{1-2}
-\endhead
-0 &Majority or most common base \\\cline{1-2}
-0.\+2 &Bases that make up atleast 20\% of the depth at a position \\\cline{1-2}
-0.\+5 &Strict or bases that make up atleast 50\% of the depth at a position \\\cline{1-2}
-0.\+9 &Strict or bases that make up atleast 90\% of the depth at a position \\\cline{1-2}
-1 &Identical or bases that make up 100\% of the depth at a position. Will have highest ambiguities \\\cline{1-2}
-\end{longtabu}
--\/m Minimum depth to call consensus(\+Default\+: 10) -\/k If \textquotesingle{}-\/k\textquotesingle{} flag is added, regions with depth less than minimum depth will not be added to the consensus sequence. Using \textquotesingle{}-\/k\textquotesingle{} will override any option specified using -\/n -\/n (N/-\/) Character to print in regions with less than minimum coverage(\+Default\+: N)
+Usage: samtools mpileup -aa -A -d 0 -Q 0 <input.bam> | ivar consensus -p <prefix>
-Output Options Description -\/p (Required) Prefix for the output fasta file and quality file
-\begin{DoxyCode}
-Example Usage:
+Note : samtools mpileup output must be piped into ivar consensus
+
+Input Options Description
+ -q Minimum quality score threshold to count base (Default: 20)
+ -t Minimum frequency threshold(0 - 1) to call consensus. (Default: 0)
+ Frequently used thresholds | Description
+ ---------------------------|------------
+ 0 | Majority or most common base
+ 0.2 | Bases that make up atleast 20% of the depth at a position
+ 0.5 | Strict or bases that make up atleast 50% of the depth at a
+ position
+ 0.9 | Strict or bases that make up atleast 90% of the depth at a
+ position
+ 1 | Identical or bases that make up 100% of the depth at a
+ position. Will have highest ambiguities
+ -m Minimum depth to call consensus(Default: 10)
+ -k If '-k' flag is added, regions with depth less than minimum depth will not be added to the
+ consensus sequence. Using '-k' will override any option specified using -n
+ -n (N/-) Character to print in regions with less than minimum coverage(Default: N)
+
+Output Options Description
+ -p (Required) Prefix for the output fasta file and quality file
\end{DoxyCode}
- samtools mpileup -\/d 1000 -\/A -\/Q 0 test.\+bam $\vert$ ivar consensus -\/p test -\/q 20 -\/t 0
+
+
+Example Usage\+:
\begin{DoxyCode}
-The command above will produce a test.fa fasta file with the consensus sequence and a test.qual.txt with
- the average quality of each base in the consensus sequence.
+samtools mpileup -d 1000 -A -Q 0 test.bam | ivar consensus -p test -q 20 -t 0
+\end{DoxyCode}
-Get primers with mismatches to the reference sequence
-----
-iVar uses a .tsv file with variants to get the zero based indices(based on the BED file) of mismatched
- primers. This command requires another .tsv file with each line containing the left and right primer names
- separated by a tab. This is used to get both the primers for an amplicon with a single mismatched primer. The
- output is a text file with the zero based primer indices delimited by a space. The output is written to a a
- text file using the prefix provided.
+The command above will produce a test.\+fa fasta file with the consensus sequence and a test.\+qual.\+txt with the average quality of each base in the consensus sequence.\hypertarget{manualpage_autotoc_md20}{}\subsection{Get primers with mismatches to the reference sequence}\label{manualpage_autotoc_md20}
+i\+Var uses a .tsv file with variants to get the zero based indices(based on the B\+E\+D file) of mismatched primers. This command requires another .tsv file with each line containing the left and right primer names separated by a tab. This is used to get both the primers for an amplicon with a single mismatched primer. The output is a text file with the zero based primer indices delimited by a space. The output is written to a a text file using the prefix provided.
-Command:
+Command\+:
+\begin{DoxyCode}
+ivar getmasked
+Usage: ivar getmasked -i <input-filtered.tsv> -b <primers.bed> -f <primer\_pairs.tsv> -p <prefix>
+Note: This step is used only for amplicon-based sequencing.
+
+Input Options Description
+ -i (Required) Input filtered variants tsv generated from 'ivar filtervariants'
+ -b (Required) BED file with primer sequences and positions
+ -f (Required) Primer pair information file containing left and right primer names for the
+ same amplicon separated by a tab
+Output Options Description
+ -p (Required) Prefix for the output text file
\end{DoxyCode}
- ivar getmasked Usage\+: ivar getmasked -\/i $<$input-\/filtered.\+tsv$>$ -\/b $<$primers.\+bed$>$ -\/f $<$primer\+\_\+pairs.\+tsv$>$ -\/p $<$prefix$>$ Note\+: This step is used only for amplicon-\/based sequencing.
-Input Options Description -\/i (Required) Input filtered variants tsv generated from \textquotesingle{}ivar filtervariants\textquotesingle{} -\/b (Required) B\+ED file with primer sequences and positions -\/f (Required) Primer pair information file containing left and right primer names for the same amplicon separated by a tab Output Options Description -\/p (Required) Prefix for the output text file
+
+Example B\+ED file
\begin{DoxyCode}
-Example BED file
+Puerto 28 52 400\_1\_out\_L 60 +
+Puerto 482 504 400\_1\_out\_R 60 -
+Puerto 359 381 400\_2\_out\_L 60 +
+Puerto 796 818 400\_2\_out\_R 60 -
+Puerto 658 680 400\_3\_out\_L* 60 +
+Puerto 1054 1076 400\_3\_out\_R* 60 -
+.
+.
+.
+.
\end{DoxyCode}
- Puerto 28 52 400\+\_\+1\+\_\+out\+\_\+L 60 + Puerto 482 504 400\+\_\+1\+\_\+out\+\_\+R 60 -\/ Puerto 359 381 400\+\_\+2\+\_\+out\+\_\+L 60 + Puerto 796 818 400\+\_\+2\+\_\+out\+\_\+R 60 -\/ Puerto 658 680 400\+\_\+3\+\_\+out\+\_\+\+L$\ast$ 60 + Puerto 1054 1076 400\+\_\+3\+\_\+out\+\_\+\+R$\ast$ 60 -\/
+
+
+Example primer pair information file
\begin{DoxyCode}
-Example primer pair information file
+400\_1\_out\_L 400\_1\_out\_R
+400\_2\_out\_L 400\_2\_out\_R
+400\_3\_out\_L 400\_3\_out\_R
+.
+.
+.
+.
\end{DoxyCode}
- 400\+\_\+1\+\_\+out\+\_\+L 400\+\_\+1\+\_\+out\+\_\+R 400\+\_\+2\+\_\+out\+\_\+L 400\+\_\+2\+\_\+out\+\_\+R 400\+\_\+3\+\_\+out\+\_\+L 400\+\_\+3\+\_\+out\+\_\+R
+
+
+Example Usage\+:
\begin{DoxyCode}
-Example Usage:
+ivar getmasked -i test.filtered.tsv -b primers.bed -f pair\_information.tsv -p test.masked.txt
\end{DoxyCode}
- ivar getmasked -\/i test.\+filtered.\+tsv -\/b primers.\+bed -\/f pair\+\_\+information.\+tsv -\/p test.\+masked.\+txt
-\begin{DoxyCode}
-The command above produces an output file - test.masked.txt.
-Example Output:
-\end{DoxyCode}
- 1 2 7 8
-\begin{DoxyCode}
-Remove reads associated with mismatched primer indices
-----
-This command accepts an aligned and sorted BAM file trimmed using `ivar trim` and removes the reads
- corresponding to the supplied primer indices, which is the output of `ivar getmasked` command. Under the hood,
- `ivar trim` adds the zero based primer index(based on the BED file) to the BAM auxillary data for every read.
- Hence, ivar removereads will only work on BAM files that have been trimmed using `ivar trim`.
+The command above produces an output file -\/ test.\+masked.\+txt.
-Command:
+Example Output\+:
+
+
+\begin{DoxyCode}
+1 2 7 8
\end{DoxyCode}
- ivar removereads
+\hypertarget{manualpage_autotoc_md21}{}\subsection{Remove reads associated with mismatched primer indices}\label{manualpage_autotoc_md21}
+This command accepts an aligned and sorted B\+AM file trimmed using {\ttfamily ivar trim} and removes the reads corresponding to the supplied primer indices, which is the output of {\ttfamily ivar getmasked} command. Under the hood, {\ttfamily ivar trim} adds the zero based primer index(based on the B\+E\+D file) to the B\+AM auxillary data for every read. Hence, ivar removereads will only work on B\+AM files that have been trimmed using {\ttfamily ivar trim}.
-Usage\+: ivar removereads -\/i $<$input.\+trimmed.\+bam$>$ -\/p $<$prefix$>$ -\/t $<$text-\/file-\/with-\/primer-\/indices$>$ Note\+: This step is used only for amplicon-\/based sequencing.
+Command\+:
+\begin{DoxyCode}
+ivar removereads
-Input Options Description -\/i (Required) Input B\+AM file trimmed with ivar trim. Must be sorted and indexed, which can be done using sort\+\_\+index\+\_\+bam.\+sh -\/t (Required) Text file with primer indices separated by spaces. This is the output of getmasked command.
+Usage: ivar removereads -i <input.trimmed.bam> -p <prefix> -t <text-file-with-primer-indices>
+Note: This step is used only for amplicon-based sequencing.
-Output Options Description -\/p (Required) Prefix for the output filtered B\+AM file
+Input Options Description
+ -i (Required) Input BAM file trimmed with ivar trim. Must be sorted and indexed, which can
+ be done using sort\_index\_bam.sh
+ -t (Required) Text file with primer indices separated by spaces. This is the output of
+ getmasked command.
+Output Options Description
+ -p (Required) Prefix for the output filtered BAM file
+\end{DoxyCode}
+
+Example Usage\+:
\begin{DoxyCode}
-Example Usage:
+ivar trim -i test.bam -p test.trimmed
+ivar removereads -i test.trimmed.bam -p test.trimmed.masked.bam -t test.masked.txt
\end{DoxyCode}
- ivar trim -\/i test.\+bam -\/p test.\+trimmed ivar removereads -\/i test.\+trimmed.\+bam -\/p test.\+trimmed.\+masked.\+bam -\/t test.\+masked.\+txt \`{}\`{}\`{}
-The {\ttfamily ivar trim} command above trims test.\+bam and produced test.\+trimmed.\+bam with the primer indice data added. The {\ttfamily ivar removereads} command produces an output file -\/ test.\+trimmed.\+masked.\+bam after removing all the reads corresponding to primer indices -\/ 1,2,7 and 8.\hypertarget{manualpage_autotoc_md20}{}\subsection{(\+Experimental) trimadapter}\label{manualpage_autotoc_md20}
+
+The {\ttfamily ivar trim} command above trims test.\+bam and produced test.\+trimmed.\+bam with the primer indice data added. The {\ttfamily ivar removereads} command produces an output file -\/ test.\+trimmed.\+masked.\+bam after removing all the reads corresponding to primer indices -\/ 1,2,7 and 8.\hypertarget{manualpage_autotoc_md22}{}\subsection{(\+Experimental) trimadapter}\label{manualpage_autotoc_md22}
{\bfseries Note\+: This feature is under active development and not completely validated yet.}
trimadapter in i\+Var can be used to trim adapter sequences from fastq files using a supplied fasta file.
\ No newline at end of file
=====================================
docs/man/man3/cookbookpage.3
=====================================
@@ -1,4 +1,4 @@
-.TH "cookbookpage" 3 "Fri Apr 3 2020" "iVar" \" -*- nroff -*-
+.TH "cookbookpage" 3 "Fri Jun 19 2020" "iVar" \" -*- nroff -*-
.ad l
.nh
.SH NAME
=====================================
docs/man/man3/installpage.3
=====================================
@@ -1,4 +1,4 @@
-.TH "installpage" 3 "Fri Apr 3 2020" "iVar" \" -*- nroff -*-
+.TH "installpage" 3 "Fri Jun 19 2020" "iVar" \" -*- nroff -*-
.ad l
.nh
.SH NAME
=====================================
docs/man/man3/manualpage.3
=====================================
@@ -1,4 +1,4 @@
-.TH "manualpage" 3 "Fri Apr 3 2020" "iVar" \" -*- nroff -*-
+.TH "manualpage" 3 "Fri Jun 19 2020" "iVar" \" -*- nroff -*-
.ad l
.nh
.SH NAME
@@ -14,7 +14,9 @@ To view detailed usage for each command type \fCivar <command>\fP Note : Command
.SS "Trim primer sequences with iVar"
iVar uses primer positions supplied in a BED file to soft clip primer sequences from an aligned and sorted BAM file\&. Following this, the reads are trimmed based on a quality threshold(Default: 20)\&. To do the quality trimming, iVar uses a sliding window approach(Default: 4)\&. The windows slides from the 5' end to the 3' end and if at any point the average base quality in the window falls below the threshold, the remaining read is soft clipped\&. If after trimming, the length of the read is greater than the minimum length specified(Default: 30), the read is written to the new trimmed BAM file\&.
.PP
-To sort and index an aligned BAM file, the following command can be used,
+Please note that the strand is taken into account while doing the trimming so forward primers are trimmed only from forward strand and reverse primers are trimmed from reverse strand\&.
+.PP
+To sort and index an aligned BAM file (OPTIONAL, if index is not present iVar will create one), the following command can be used,
.PP
.PP
.nf
@@ -115,7 +117,7 @@ Output Options Description
Example Usage:
.PP
.nf
-samtools mpileup -A -d 600000 -B -Q 0 test\&.trimmed\&.bam | ivar variants -p test -q 20 -t 0\&.03 -r test_reference\&.fa -g test\&.gff
+samtools mpileup -aa -A -d 600000 -B -Q 0 test\&.trimmed\&.bam | ivar variants -p test -q 20 -t 0\&.03 -r test_reference\&.fa -g test\&.gff
.fi
.PP
@@ -216,93 +218,141 @@ Minimum frequency threshold Consensus 0 T 0\&.5 T 0\&.6 T 0\&.7 D(A
.PP
The output of the command is a fasta file with the consensus sequence and a \&.txt file with the average quality of every base used to generate the consensus at each position\&. \fIFor insertions, the quality is set to be the minimum quality threshold since mpileup doesn't give the quality of bases in insertions\&.\fP
.PP
-Command: ``` ivar consensus
+Command:
.PP
-Usage: samtools mpileup -aa -A -d 0 -Q 0 <input\&.bam> | ivar consensus -p <prefix>
.PP
-Note : samtools mpileup output must be piped into \fCivar consensus\fP
+.nf
+ivar consensus
+
+Usage: samtools mpileup -aa -A -d 0 -Q 0 <input\&.bam> | ivar consensus -p <prefix>
+
+Note : samtools mpileup output must be piped into ivar consensus
+
+Input Options Description
+ -q Minimum quality score threshold to count base (Default: 20)
+ -t Minimum frequency threshold(0 - 1) to call consensus\&. (Default: 0)
+ Frequently used thresholds | Description
+ ---------------------------|------------
+ 0 | Majority or most common base
+ 0\&.2 | Bases that make up atleast 20% of the depth at a position
+ 0\&.5 | Strict or bases that make up atleast 50% of the depth at a position
+ 0\&.9 | Strict or bases that make up atleast 90% of the depth at a position
+ 1 | Identical or bases that make up 100% of the depth at a position\&. Will have highest ambiguities
+ -m Minimum depth to call consensus(Default: 10)
+ -k If '-k' flag is added, regions with depth less than minimum depth will not be added to the consensus sequence\&. Using '-k' will override any option specified using -n
+ -n (N/-) Character to print in regions with less than minimum coverage(Default: N)
+
+Output Options Description
+ -p (Required) Prefix for the output fasta file and quality file
+.fi
.PP
-Input Options Description -q Minimum quality score threshold to count base (Default: 20) -t Minimum frequency threshold(0 - 1) to call consensus\&. (Default: 0) Frequently used thresholds Description 0 Majority or most common base 0\&.2 Bases that make up atleast 20% of the depth at a position 0\&.5 Strict or bases that make up atleast 50% of the depth at a position 0\&.9 Strict or bases that make up atleast 90% of the depth at a position 1 Identical or bases that make up 100% of the depth at a position\&. Will have highest ambiguities -m Minimum depth to call consensus(Default: 10) -k If '-k' flag is added, regions with depth less than minimum depth will not be added to the consensus sequence\&. Using '-k' will override any option specified using -n -n (N/-) Character to print in regions with less than minimum coverage(Default: N)
.PP
-Output Options Description -p (Required) Prefix for the output fasta file and quality file
+Example Usage:
.PP
.nf
-Example Usage:
+samtools mpileup -d 1000 -A -Q 0 test\&.bam | ivar consensus -p test -q 20 -t 0
.fi
.PP
- samtools mpileup -d 1000 -A -Q 0 test\&.bam | ivar consensus -p test -q 20 -t 0
.PP
-.nf
The command above will produce a test\&.fa fasta file with the consensus sequence and a test\&.qual\&.txt with the average quality of each base in the consensus sequence\&.
-
-Get primers with mismatches to the reference sequence
-----
-
+.SS "Get primers with mismatches to the reference sequence"
iVar uses a \&.tsv file with variants to get the zero based indices(based on the BED file) of mismatched primers\&. This command requires another \&.tsv file with each line containing the left and right primer names separated by a tab\&. This is used to get both the primers for an amplicon with a single mismatched primer\&. The output is a text file with the zero based primer indices delimited by a space\&. The output is written to a a text file using the prefix provided\&.
+.PP
+Command:
+.PP
+.nf
+ivar getmasked
+Usage: ivar getmasked -i <input-filtered\&.tsv> -b <primers\&.bed> -f <primer_pairs\&.tsv> -p <prefix>
+Note: This step is used only for amplicon-based sequencing\&.
-Command:
+Input Options Description
+ -i (Required) Input filtered variants tsv generated from 'ivar filtervariants'
+ -b (Required) BED file with primer sequences and positions
+ -f (Required) Primer pair information file containing left and right primer names for the same amplicon separated by a tab
+Output Options Description
+ -p (Required) Prefix for the output text file
.fi
.PP
- ivar getmasked Usage: ivar getmasked -i <input-filtered\&.tsv> -b <primers\&.bed> -f <primer_pairs\&.tsv> -p <prefix> Note: This step is used only for amplicon-based sequencing\&.
.PP
-Input Options Description -i (Required) Input filtered variants tsv generated from 'ivar filtervariants' -b (Required) BED file with primer sequences and positions -f (Required) Primer pair information file containing left and right primer names for the same amplicon separated by a tab Output Options Description -p (Required) Prefix for the output text file
+Example BED file
.PP
.PP
.nf
-Example BED file
+Puerto 28 52 400_1_out_L 60 +
+Puerto 482 504 400_1_out_R 60 -
+Puerto 359 381 400_2_out_L 60 +
+Puerto 796 818 400_2_out_R 60 -
+Puerto 658 680 400_3_out_L* 60 +
+Puerto 1054 1076 400_3_out_R* 60 -
+\&.
+\&.
+\&.
+\&.
.fi
.PP
- Puerto 28 52 400_1_out_L 60 + Puerto 482 504 400_1_out_R 60 - Puerto 359 381 400_2_out_L 60 + Puerto 796 818 400_2_out_R 60 - Puerto 658 680 400_3_out_L* 60 + Puerto 1054 1076 400_3_out_R* 60 -
+.PP
+Example primer pair information file
.PP
.nf
-Example primer pair information file
+400_1_out_L 400_1_out_R
+400_2_out_L 400_2_out_R
+400_3_out_L 400_3_out_R
+\&.
+\&.
+\&.
+\&.
.fi
.PP
- 400_1_out_L 400_1_out_R 400_2_out_L 400_2_out_R 400_3_out_L 400_3_out_R
+.PP
+Example Usage:
.PP
.nf
-Example Usage:
+ivar getmasked -i test\&.filtered\&.tsv -b primers\&.bed -f pair_information\&.tsv -p test\&.masked\&.txt
.fi
.PP
- ivar getmasked -i test\&.filtered\&.tsv -b primers\&.bed -f pair_information\&.tsv -p test\&.masked\&.txt
.PP
-.nf
The command above produces an output file - test\&.masked\&.txt\&.
-
+.PP
Example Output:
-
+.PP
+.PP
+.nf
+1 2 7 8
.fi
.PP
- 1 2 7 8
+.SS "Remove reads associated with mismatched primer indices"
+This command accepts an aligned and sorted BAM file trimmed using \fCivar trim\fP and removes the reads corresponding to the supplied primer indices, which is the output of \fCivar getmasked\fP command\&. Under the hood, \fCivar trim\fP adds the zero based primer index(based on the BED file) to the BAM auxillary data for every read\&. Hence, ivar removereads will only work on BAM files that have been trimmed using \fCivar trim\fP\&.
+.PP
+Command:
.PP
.nf
-Remove reads associated with mismatched primer indices
-----
+ivar removereads
-This command accepts an aligned and sorted BAM file trimmed using `ivar trim` and removes the reads corresponding to the supplied primer indices, which is the output of `ivar getmasked` command\&. Under the hood, `ivar trim` adds the zero based primer index(based on the BED file) to the BAM auxillary data for every read\&. Hence, ivar removereads will only work on BAM files that have been trimmed using `ivar trim`\&.
+Usage: ivar removereads -i <input\&.trimmed\&.bam> -p <prefix> -t <text-file-with-primer-indices>
+Note: This step is used only for amplicon-based sequencing\&.
-Command:
+Input Options Description
+ -i (Required) Input BAM file trimmed with ivar trim\&. Must be sorted and indexed, which can be done using sort_index_bam\&.sh
+ -t (Required) Text file with primer indices separated by spaces\&. This is the output of getmasked command\&.
+
+Output Options Description
+ -p (Required) Prefix for the output filtered BAM file
.fi
.PP
- ivar removereads
-.PP
-Usage: ivar removereads -i <input\&.trimmed\&.bam> -p <prefix> -t <text-file-with-primer-indices> Note: This step is used only for amplicon-based sequencing\&.
-.PP
-Input Options Description -i (Required) Input BAM file trimmed with ivar trim\&. Must be sorted and indexed, which can be done using sort_index_bam\&.sh -t (Required) Text file with primer indices separated by spaces\&. This is the output of getmasked command\&.
-.PP
-Output Options Description -p (Required) Prefix for the output filtered BAM file
.PP
+Example Usage:
.PP
.nf
-Example Usage:
+ivar trim -i test\&.bam -p test\&.trimmed
+ivar removereads -i test\&.trimmed\&.bam -p test\&.trimmed\&.masked\&.bam -t test\&.masked\&.txt
+
.fi
.PP
- ivar trim -i test\&.bam -p test\&.trimmed ivar removereads -i test\&.trimmed\&.bam -p test\&.trimmed\&.masked\&.bam -t test\&.masked\&.txt ```
.PP
The \fCivar trim\fP command above trims test\&.bam and produced test\&.trimmed\&.bam with the primer indice data added\&. The \fCivar removereads\fP command produces an output file - test\&.trimmed\&.masked\&.bam after removing all the reads corresponding to primer indices - 1,2,7 and 8\&.
.SS "(Experimental) trimadapter"
=====================================
src/call_consensus_pileup.cpp
=====================================
@@ -130,13 +130,17 @@ ret_t get_consensus_allele(std::vector<allele> ad, uint8_t min_qual, double thre
return t;
}
-int call_consensus_from_plup(std::istream &cin, std::string out_file, uint8_t min_qual, double threshold, uint8_t min_depth, char gap, bool min_coverage_flag){
+int call_consensus_from_plup(std::istream &cin, std::string seq_id, std::string out_file, uint8_t min_qual, double threshold, uint8_t min_depth, char gap, bool min_coverage_flag){
std::string line, cell;
std::ofstream fout((out_file+".fa").c_str());
std::ofstream tmp_qout((out_file+".qual.txt").c_str());
char *o = new char[out_file.length() + 1];
strcpy(o, out_file.c_str());
- fout << ">Consensus_" << basename(o) << "_threshold_" << threshold << "_quality_" << (uint16_t) min_qual <<std::endl;
+ if(seq_id.empty()) {
+ fout << ">Consensus_" << basename(o) << "_threshold_" << threshold << "_quality_" << (uint16_t) min_qual <<std::endl;
+ } else {
+ fout << ">" << seq_id <<std::endl;
+ }
delete [] o;
int ctr = 0, mdepth = 0;
uint32_t prev_pos = 0, pos = 0;
=====================================
src/call_consensus_pileup.h
=====================================
@@ -19,7 +19,7 @@ struct ret_t {
};
void format_alleles(std::vector<allele> &ad);
-int call_consensus_from_plup(std::istream &cin, std::string out_file, uint8_t min_qual, double threshold, uint8_t min_depth, char gap, bool min_coverage_flag);
+int call_consensus_from_plup(std::istream &cin, std::string seq_id, std::string out_file, uint8_t min_qual, double threshold, uint8_t min_depth, char gap, bool min_coverage_flag);
ret_t get_consensus_allele(std::vector<allele> ad, uint8_t min_qual, double threshold, char gap);
#endif
=====================================
src/call_variants.cpp
=====================================
@@ -86,12 +86,19 @@ int call_variants_from_plup(std::istream &cin, std::string out_file, uint8_t min
}
ctr++;
}
- if(mdepth < min_depth) { // Check for minimum depth
+ ad = update_allele_depth(ref, bases, qualities, min_qual);
+ if(ad.size() == 0){
line_stream.clear();
continue;
}
- ad = update_allele_depth(ref, bases, qualities, min_qual);
- if(ad.size() == 0){
+ // Get ungapped depth
+ pdepth = 0;
+ for(std::vector<allele>::iterator it = ad.begin(); it != ad.end(); ++it) {
+ if(it->nuc[0]=='*' || it->nuc[0] == '+' || it->nuc[0] == '-')
+ continue;
+ pdepth += it->depth;
+ }
+ if(pdepth < min_depth) { // Check for minimum depth
line_stream.clear();
continue;
}
@@ -105,13 +112,6 @@ int call_variants_from_plup(std::istream &cin, std::string out_file, uint8_t min
ad.push_back(a);
ref_it = ad.end() - 1;
}
- // Get ungapped coverage
- pdepth = 0;
- for(std::vector<allele>::iterator it = ad.begin(); it != ad.end(); ++it) {
- if(it->nuc[0]=='*' || it->nuc[0] == '+' || it->nuc[0] == '-')
- continue;
- pdepth += it->depth;
- }
for(std::vector<allele>::iterator it = ad.begin(); it != ad.end(); ++it) {
if((*it == *ref_it) || it->nuc[0]=='*')
continue;
@@ -150,6 +150,7 @@ int call_variants_from_plup(std::istream &cin, std::string out_file, uint8_t min
}
out_str.str("");
out_str.clear();
+ delete[] freq_depth;
}
line_stream.clear();
}
=====================================
src/get_masked_amplicons.cpp
=====================================
@@ -35,12 +35,22 @@ int get_primers_with_mismatches(std::string bed, std::string vpath, std::string
}
pos--; // 1 based to 0 based
tmp = get_primers(primers, pos);
+ // mismatches_primers.insert(mismatches_primers.end(), tmp.begin(), tmp.end());
+ std::vector<primer>::iterator tmp_it;
for(std::vector<primer>::iterator it = tmp.begin(); it != tmp.end(); ++it) {
- std::vector<primer>::iterator tmp_it = std::find(mismatches_primers.begin(), mismatches_primers.end(), *it);
- if(tmp_it == tmp.end())
+ tmp_it = std::find(mismatches_primers.begin(), mismatches_primers.end(), *it);
+ if(tmp_it == mismatches_primers.end()){
+ std::cout << it->get_name() << std::endl;
mismatches_primers.push_back(*it);
+ }
+ // Look for primer pair
+ if(it->get_pair_indice() != -1){
+ tmp_it = std::find(mismatches_primers.begin(), mismatches_primers.end(), primers.at(it->get_pair_indice()));
+ if(tmp_it == mismatches_primers.end()){
+ mismatches_primers.push_back(primers.at(it->get_pair_indice()));
+ }
+ }
}
- mismatches_primers.insert(mismatches_primers.end(), tmp.begin(), tmp.end());
line_stream.clear();
tmp.clear();
}
=====================================
src/ivar.cpp
=====================================
@@ -16,12 +16,13 @@
#include "suffix_tree.h"
#include "get_common_variants.h"
-const std::string VERSION = "1.2.2";
+const std::string VERSION = "1.2.3";
struct args_t {
std::string bam; // -i
std::string bed; // -b
std::string text; // -t
+ std::string seq_id; // -i for consensus
std::string prefix; // -p
std::string ref; // -r
std::string region; // -R
@@ -38,6 +39,7 @@ struct args_t {
std::string primer_pair_file; // -f
std::string file_list; // -f
bool write_no_primers_flag; // -e
+ bool mark_qcfail_flag; // -f
std::string gff; // -g
} g_args;
@@ -66,7 +68,8 @@ void print_trim_usage(){
" -m Minimum length of read to retain after trimming (Default: 30)\n"
" -q Minimum quality threshold for sliding window to pass (Default: 20)\n"
" -s Width of sliding window (Default: 4)\n"
- " -e Include reads with no primers. By default, reads with no primers are excluded\n\n"
+ " -e Include reads with no primers. By default, reads with no primers are excluded\n"
+ " -f Mark reads as QCFAIL instead of excluding them\n\n"
"Output Options Description\n"
" -p (Required) Prefix for the output BAM file\n";
}
@@ -114,7 +117,8 @@ void print_consensus_usage(){
" -k If '-k' flag is added, regions with depth less than minimum depth will not be added to the consensus sequence. Using '-k' will override any option specified using -n \n"
" -n (N/-) Character to print in regions with less than minimum coverage(Default: N)\n\n"
"Output Options Description\n"
- " -p (Required) Prefix for the output fasta file and quality file\n";
+ " -p (Required) Prefix for the output fasta file and quality file\n"
+ " -i (Optional) Name of fasta header. By default, the prefix is used to create the fasta header in the following format, Consensus_<prefix>_threshold_<frequency-threshold>_quality_<minimum-quality>\n";
}
void print_removereads_usage(){
@@ -158,9 +162,9 @@ void print_version_info(){
"\nPlease raise issues and bug reports at https://github.com/andersen-lab/ivar/\n\n";
}
-static const char *trim_opt_str = "i:b:p:m:q:s:eh?";
+static const char *trim_opt_str = "i:b:p:m:q:s:efh?";
static const char *variants_opt_str = "p:t:q:m:r:g:h?";
-static const char *consensus_opt_str = "p:q:t:m:n:kh?";
+static const char *consensus_opt_str = "i:p:q:t:m:n:kh?";
static const char *removereads_opt_str = "i:p:t:b:h?";
static const char *filtervariants_opt_str = "p:t:f:h?";
static const char *getmasked_opt_str = "i:b:f:p:h?";
@@ -210,6 +214,7 @@ int main(int argc, char* argv[]){
g_args.sliding_window = 4;
g_args.min_length = 30;
g_args.write_no_primers_flag = false;
+ g_args.mark_qcfail_flag = false;
g_args.bed = "";
opt = getopt( argc, argv, trim_opt_str);
while( opt != -1 ) {
@@ -232,10 +237,13 @@ int main(int argc, char* argv[]){
case 's':
g_args.sliding_window = std::stoi(optarg);
break;
- case 'h':
case 'e':
g_args.write_no_primers_flag = true;
break;
+ case 'f':
+ g_args.mark_qcfail_flag = true;
+ break;
+ case 'h':
case '?':
print_trim_usage();
return -1;
@@ -248,7 +256,7 @@ int main(int argc, char* argv[]){
return -1;
}
g_args.prefix = get_filename_without_extension(g_args.prefix,".bam");
- res = trim_bam_qual_primer(g_args.bam, g_args.bed, g_args.prefix, g_args.region, g_args.min_qual, g_args.sliding_window, cl_cmd.str(), g_args.write_no_primers_flag, g_args.min_length);
+ res = trim_bam_qual_primer(g_args.bam, g_args.bed, g_args.prefix, g_args.region, g_args.min_qual, g_args.sliding_window, cl_cmd.str(), g_args.write_no_primers_flag, g_args.mark_qcfail_flag, g_args.min_length);
} else if (cmd.compare("variants") == 0){
g_args.min_qual = 20;
g_args.min_threshold = 0.03;
@@ -307,6 +315,7 @@ int main(int argc, char* argv[]){
res = call_variants_from_plup(std::cin, g_args.prefix, g_args.min_qual, g_args.min_threshold, g_args.min_depth, g_args.ref, g_args.gff);
} else if (cmd.compare("consensus") == 0){
opt = getopt( argc, argv, consensus_opt_str);
+ g_args.seq_id = "";
g_args.min_threshold = 0;
g_args.min_depth = 10;
g_args.gap = 'N';
@@ -317,6 +326,9 @@ int main(int argc, char* argv[]){
case 't':
g_args.min_threshold = atof(optarg);
break;
+ case 'i':
+ g_args.seq_id = optarg;
+ break;
case 'p':
g_args.prefix = optarg;
break;
@@ -360,7 +372,7 @@ int main(int argc, char* argv[]){
std::cout << "Regions with depth less than minimum depth will not added to consensus" << std::endl;
else
std::cout << "Regions with depth less than minimum depth covered by: " << g_args.gap << std::endl;
- res = call_consensus_from_plup(std::cin, g_args.prefix, g_args.min_qual, g_args.min_threshold, g_args.min_depth, g_args.gap, g_args.keep_min_coverage);
+ res = call_consensus_from_plup(std::cin, g_args.seq_id, g_args.prefix, g_args.min_qual, g_args.min_threshold, g_args.min_depth, g_args.gap, g_args.keep_min_coverage);
} else if (cmd.compare("removereads") == 0){
opt = getopt( argc, argv, removereads_opt_str);
while( opt != -1 ) {
=====================================
src/primer_bed.cpp
=====================================
@@ -80,7 +80,6 @@ void primer::add_read_count(uint32_t rc){
read_count += rc;
}
-
void print_bed_format(){
std::cout << "iVar uses the standard 6 column BED format as defined here - https://genome.ucsc.edu/FAQ/FAQformat.html#format1." << std::endl;
std::cout << "It requires the following columns delimited by a tab: chrom, chromStart, chromEnd, name, score, strand" << std::endl;
@@ -96,6 +95,7 @@ std::vector<primer> populate_from_file(std::string path){
std::string cell;
int ctr = 0;
primer p;
+ p.set_strand(0); // Set strand to NULL
while(std::getline(lineStream,cell,'\t')){
switch(ctr){
case 0:
@@ -126,9 +126,9 @@ std::vector<primer> populate_from_file(std::string path){
if(std::all_of(cell.begin(), cell.end(), ::isdigit)) {
p.set_score(stoi(cell));
} else {
- print_bed_format();
- primers.clear();
- return primers;
+ print_bed_format(); // score is missing, send warning but continue populating
+ std::cout << "\nWARNING: The BED file provided did not have the expected score column, but iVar will continue trimming\n" << std::endl;
+ p.set_score(-1);
}
break;
case 5:
@@ -142,16 +142,14 @@ std::vector<primer> populate_from_file(std::string path){
}
ctr++;
}
+ if(indice == 0 && ctr < 6)
+ std::cout << "Strand not found in primer BED file so strand will not be considered for trimming" << std::endl;
p.set_indice(indice);
p.set_pair_indice(-1);
p.set_read_count(0);
primers.push_back(p);
indice++;
}
- if(primers.size() == 0){
- print_bed_format();
- std::cout << "Found 0 primers in BED file" << std::endl;
- }
std::cout << "Found " << primers.size() << " primers in BED file" << std::endl;
return primers;
}
@@ -217,8 +215,7 @@ primer get_min_start(std::vector<primer> primers){
return *(minmax_start.first);
}
-
primer get_max_end(std::vector<primer> primers){
auto minmax_start = std::minmax_element(primers.begin(), primers.end(), [] (primer lhs, primer rhs) {return lhs.get_end() < rhs.get_end();});
return *(minmax_start.second);
-}
+}
\ No newline at end of file
=====================================
src/ref_seq.cpp
=====================================
@@ -2,12 +2,13 @@
char ref_antd::get_base(int64_t pos, std::string region){ // 1-based position
int len;
+ char base = 0;
if(!region.empty() && this->fai != NULL){
seq = fai_fetch(this->fai, region.c_str(), &len);
}
- if(seq == NULL)
- return 0;
- return *(seq + (pos - 1));
+ if (seq) base = *(seq + (pos - 1));
+ free(seq);
+ return base;
}
char* ref_antd::get_codon(int64_t pos, std::string region, gff3_feature feature){
@@ -35,6 +36,7 @@ char* ref_antd::get_codon(int64_t pos, std::string region, gff3_feature feature)
codon[i] = *(seq + codon_start_pos + i - edit_offset);
}
}
+ free(seq);
return codon;
}
@@ -68,6 +70,7 @@ 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;
+ free(seq);
return codon;
}
@@ -101,6 +104,11 @@ ref_antd::ref_antd(std::string ref_path, std::string gff_path){
this->add_gff(gff_path);
}
+ref_antd::~ref_antd()
+{
+ if (this->fai) fai_destroy(this->fai);
+}
+
int ref_antd::codon_aa_stream(std::string region, std::ostringstream &line_stream, std::ofstream &fout, int64_t pos, char alt){
std::vector<gff3_feature> features = gff.query_features(pos, "CDS");
if(features.size() == 0){ // No matching CDS
@@ -119,6 +127,9 @@ int ref_antd::codon_aa_stream(std::string region, std::ostringstream &line_strea
fout << alt_codon[0] << alt_codon[1] << alt_codon[2] << "\t";
fout << codon2aa(alt_codon[0], alt_codon[1], alt_codon[2]);
fout << std::endl;
+
+ delete[] ref_codon;
+ delete[] alt_codon;
}
return 0;
}
=====================================
src/ref_seq.h
=====================================
@@ -15,6 +15,7 @@ class ref_antd{
public:
ref_antd(std::string ref_path);
ref_antd(std::string ref_path, std::string gff_path);
+ ~ref_antd();
char get_base(int64_t pos, std::string region);
int add_gff(std::string path);
int add_seq(std::string path);
=====================================
src/trim_primer_quality.cpp
=====================================
@@ -84,8 +84,13 @@ cigar_ quality_trim(bam1_t* r, uint8_t qual_threshold, uint8_t sliding_window){
int del_len, cig, temp;
uint32_t i = 0, j = 0;
cigar_ t;
- while(m >= qual_threshold && (i < (uint32_t)r->core.l_qseq)){
+ init_cigar(&t);
+ if(0 > r->core.l_qseq - sliding_window)
+ sliding_window = (uint32_t)r->core.l_qseq;
+ while(i < (uint32_t)r->core.l_qseq){
m = mean_quality(qual, i, i+sliding_window);
+ if(m < qual_threshold)
+ break;
i++;
if(i > (uint32_t)r->core.l_qseq - sliding_window)
sliding_window--;
@@ -97,7 +102,9 @@ cigar_ quality_trim(bam1_t* r, uint8_t qual_threshold, uint8_t sliding_window){
del_len = r->core.l_qseq - i;
start_pos = get_pos_on_reference(cigar, r->core.n_cigar, del_len, r->core.pos); // For reverse reads need to set core->pos.
if(reverse && start_pos <= r->core.pos) {
+ free(ncigar);
t.cigar = cigar;
+ t.free_cig = false;
t.nlength = r->core.n_cigar;
t.start_pos = r->core.pos;
return t;
@@ -140,6 +147,7 @@ cigar_ quality_trim(bam1_t* r, uint8_t qual_threshold, uint8_t sliding_window){
}
t.cigar = ncigar;
t.nlength = j;
+ t.free_cig = true;
t.start_pos = start_pos;
return t;
}
@@ -227,6 +235,7 @@ cigar_ primer_trim(bam1_t *r, int32_t new_pos, bool unpaired_rev = false){
}
return {
ncigar,
+ true,
j,
start_pos
};
@@ -259,7 +268,7 @@ void get_overlapping_primers(bam1_t* r, std::vector<primer> primers, std::vector
start_pos = r->core.pos;
}
for(std::vector<primer>::iterator it = primers.begin(); it != primers.end(); ++it) {
- if(start_pos >= it->get_start() && start_pos <= it->get_end() && strand == it->get_strand())
+ if(start_pos >= it->get_start() && start_pos <= it->get_end() && (strand == it->get_strand() || it->get_strand() == 0))
overlapped_primers.push_back(*it);
}
}
@@ -276,31 +285,27 @@ void get_overlapping_primers(bam1_t* r, std::vector<primer> primers, std::vector
start_pos = r->core.pos;
}
for(std::vector<primer>::iterator it = primers.begin(); it != primers.end(); ++it) {
- if(start_pos >= it->get_start() && start_pos <= it->get_end() && strand == it->get_strand())
+ if(start_pos >= it->get_start() && start_pos <= it->get_end() && (strand == it->get_strand() ||it->get_strand() == 0))
overlapped_primers.push_back(*it);
}
}
-cigar_ condense_cigar(uint32_t* cigar, uint32_t n){
+void condense_cigar(cigar_ *t){
uint32_t i = 0, len = 0, cig, next_cig;
- cigar_ t;
- while(i< n -1){
- cig = bam_cigar_op(cigar[i]);
- next_cig = bam_cigar_op(cigar[i+1]);
+ while(i< t->nlength -1){
+ cig = bam_cigar_op(t->cigar[i]);
+ next_cig = bam_cigar_op(t->cigar[i+1]);
if(cig == next_cig){
- len = bam_cigar_oplen(cigar[i])+bam_cigar_oplen(cigar[i+1]);
- cigar[i] = bam_cigar_gen(len, bam_cigar_op(cigar[i]));
- for(uint32_t j = i+1; j < n - 1; j++){
- cigar[j] = cigar[j+1];
+ len = bam_cigar_oplen(t->cigar[i])+bam_cigar_oplen(t->cigar[i+1]);
+ t->cigar[i] = bam_cigar_gen(len, bam_cigar_op(t->cigar[i]));
+ for(uint32_t j = i+1; j < t->nlength - 1; j++){
+ t->cigar[j] = t->cigar[j+1];
}
- n--;
+ t->nlength--;
} else {
i++;
}
}
- t.cigar = cigar;
- t.nlength = n;
- return t;
}
void add_pg_line_to_header(bam_hdr_t** hdr, char *cmd){
@@ -315,16 +320,17 @@ void add_pg_line_to_header(bam_hdr_t** hdr, char *cmd){
(*hdr)->l_text = len-1;
}
-int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, std::string region_, uint8_t min_qual, uint8_t sliding_window, std::string cmd, bool write_no_primer_reads, int min_length = 30){
+int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, std::string region_, uint8_t min_qual, uint8_t sliding_window, std::string cmd, bool write_no_primer_reads, bool mark_qcfail_flag, int min_length = 30) {
std::vector<primer> primers;
if(!bed.empty()){
primers = populate_from_file(bed);
+ if(primers.size() == 0){
+ std::cout << "Exiting." << std::endl;
+ return -1;
+ }
}
- // if(primers.size() == 0){
- // return 0;
- // }
if(bam.empty()){
- std::cout << "Bam file in empty." << std::endl;
+ std::cout << "Bam file is empty." << std::endl;
return -1;
}
bam_out += ".bam";
@@ -399,6 +405,7 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out,
bam1_t *aln = bam_init1();
int ctr = 0;
cigar_ t;
+ init_cigar(&t);
uint32_t primer_trim_count = 0, no_primer_counter = 0, low_quality = 0;
bool unmapped_flag = false;
uint32_t unmapped_counter = 0;
@@ -423,6 +430,7 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out,
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())
@@ -431,7 +439,7 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out,
t = quality_trim(aln, min_qual, sliding_window); // Quality Trimming
if(bam_is_rev(aln))
aln->core.pos = t.start_pos;
- t = condense_cigar(t.cigar, t.nlength);
+ condense_cigar(&t);
// aln->core.pos += t.start_pos;
replace_cigar(aln, t.nlength, t.cigar);
} else { // Unpaired reads: Might be stitched reads
@@ -461,7 +469,7 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out,
cit->add_read_count(1);
}
t = quality_trim(aln, min_qual, sliding_window); // Quality Trimming
- t = condense_cigar(t.cigar, t.nlength);
+ condense_cigar(&t);
replace_cigar(aln, t.nlength, t.cigar);
}
if(primer_trimmed){
@@ -488,6 +496,7 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out,
}
} else {
if((primers.size() == 0 || write_no_primer_reads) && !unmapped_flag){ // Write mapped reads to BAM if -e flag given
+ if (mark_qcfail_flag) aln->core.flag |= BAM_FQCFAIL;
if(bam_write1(out, aln) < 0){
std::cout << "Not able to write to BAM" << std::endl;
hts_itr_destroy(iter);
@@ -503,6 +512,19 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out,
}
} else {
low_quality++;
+ if (mark_qcfail_flag) {
+ aln->core.flag |= BAM_FQCFAIL;
+ if (bam_write1(out, aln) < 0){
+ std::cout << "Not able to write to BAM" << std::endl;
+ hts_itr_destroy(iter);
+ hts_idx_destroy(idx);
+ bam_destroy1(aln);
+ bam_hdr_destroy(header);
+ sam_close(in);
+ bgzf_close(out);
+ return -1;
+ }
+ }
}
ctr++;
if(ctr % log_skip == 0){
@@ -516,9 +538,16 @@ int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out,
std::cout << cit->get_name() << "\t" << cit->get_read_count() << std::endl;
}
std::cout << std::endl << "Trimmed primers from " << round_int(primer_trim_count, mapped) << "% (" << primer_trim_count << ") of reads." << std::endl;
- std::cout << round_int( low_quality, mapped) << "% (" << low_quality << ") of reads were quality trimmed below the minimum length of " << min_length << " bp and were not writen to file." << std::endl;
+ std::cout << round_int( low_quality, mapped) << "% (" << low_quality << ") of reads were quality trimmed below the minimum length of " << min_length << " bp and were ";
+ if (mark_qcfail_flag) {
+ std::cout << "marked as failed" << std::endl;
+ } else {
+ std::cout << "not written to file." << std::endl;
+ }
if(write_no_primer_reads){
- std::cout << round_int(no_primer_counter, mapped) << "% (" << no_primer_counter << ") of reads started outside of primer regions. Since the -e flag was given, these reads were written to file." << std::endl;
+ std::cout << round_int(no_primer_counter, mapped) << "% (" << no_primer_counter << ") of reads started outside of primer regions. Since the -e flag was given, these reads were written to file";
+ if (mark_qcfail_flag) std::cout << " and the BAM_QCFAIL flag set";
+ std::cout << "." << std::endl;
} else if (primers.size() == 0) {
std::cout << round_int(no_primer_counter, mapped) << "% (" << no_primer_counter << ") of reads started outside of primer regions. Since there were no primers found in BED file, these reads were written to file." << std::endl;
} else {
=====================================
src/trim_primer_quality.h
=====================================
@@ -19,13 +19,19 @@
struct cigar_ {
uint32_t *cigar;
+ bool free_cig;
uint32_t nlength;
int32_t start_pos;
};
+inline void init_cigar(cigar_ *t) { t->cigar=NULL; t->free_cig=false; t->nlength=0; t->start_pos=0; }
+inline void free_cigar(cigar_ t) { if (t.free_cig) free(t.cigar); }
+
void add_pg_line_to_header(bam_hdr_t** hdr, char *cmd);
-int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, std::string region_, uint8_t min_qual, uint8_t sliding_window, std::string cmd, bool write_no_primer_reads, int min_length);
+
+int trim_bam_qual_primer(std::string bam, std::string bed, std::string bam_out, std::string region_, uint8_t min_qual, uint8_t sliding_window, std::string cmd, bool write_no_primer_reads, bool mark_qcfail_flag, int min_length);
+void free_cigar(cigar_ t);
int32_t get_pos_on_query(uint32_t *cigar, uint32_t ncigar, int32_t pos, int32_t ref_start);
int32_t get_pos_on_reference(uint32_t *cigar, uint32_t ncigar, uint32_t pos, uint32_t ref_start);
void reverse_qual(uint8_t *q, int l);
@@ -35,7 +41,7 @@ cigar_ quality_trim(bam1_t* r, uint8_t qual_threshold, uint8_t sliding_window);
void print_cigar(uint32_t *cigar, int nlength);
cigar_ primer_trim(bam1_t *r, int32_t new_pos, bool unpaired_rev);
void replace_cigar(bam1_t *b, uint32_t n, uint32_t *cigar);
-cigar_ condense_cigar(uint32_t* cigar, uint32_t n);
+void condense_cigar(cigar_ *t);
void get_overlapping_primers(bam1_t* r, std::vector<primer> primers, std::vector<primer> &overlapping_primers);
void get_overlapping_primers(bam1_t* r, std::vector<primer> primers, std::vector<primer> &overlapping_primers, bool unpaired_rev);
=====================================
tests/Makefile.am
=====================================
@@ -2,14 +2,16 @@ LIBS = -lhts -lz -lpthread
CXXFLAGS = -g -std=c++11 -Wall -Wextra -Werror
-TESTS = check_primer_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases
-check_PROGRAMS = check_primer_trim check_quality_trim check_consensus check_allele_depth check_consensus_threshold check_consensus_min_depth check_primer_bed check_getmasked check_removereads check_variants check_common_variants check_unpaired_trim check_primer_trim_edge_cases
+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_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_primer_trim_SOURCES = test_primer_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp
+check_trim_SOURCES = test_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.cpp
check_quality_trim_SOURCES = check_quality_trim.cpp ../src/trim_primer_quality.cpp ../src/primer_bed.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
check_consensus_min_depth_SOURCES = test_consensus_min_depth.cpp ../src/call_consensus_pileup.cpp ../src/allele_functions.cpp
+check_consensus_seq_id_SOURCES = test_consensus_seq_id.cpp ../src/call_consensus_pileup.cpp ../src/allele_functions.cpp
check_variants_SOURCES = test_variants.cpp ../src/call_variants.cpp ../src/allele_functions.cpp ../src/parse_gff.cpp ../src/ref_seq.cpp
check_common_variants_SOURCES = test_common_variants.cpp ../src/get_common_variants.cpp
check_primer_bed_SOURCES = test_primer_bed.cpp ../src/primer_bed.cpp
=====================================
tests/check_quality_trim.cpp
=====================================
@@ -15,8 +15,8 @@ int main(){
iter = sam_itr_querys(idx, header, region_.c_str());
bam1_t *aln = bam_init1();
cigar_ t;
- int lengths[6] = {150,100,100, 25,146,144}, ctr = 0;
- int start_pos_rev[6] = {19, 113, 208, 324, 199, 231};
+ int lengths[6] = {150,99,99, 24,146,144}, ctr = 0;
+ int start_pos_rev[6] = {19, 114, 209, 325, 199, 231};
while(sam_itr_next(in, iter, aln) >= 0) {
t = quality_trim(aln, 20, 4);
std::cout << bam_get_qname(aln) << std::endl;
@@ -27,7 +27,14 @@ int main(){
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_consensus_min_depth.cpp
=====================================
@@ -4,10 +4,10 @@
#include "../src/call_consensus_pileup.h"
#include "../src/allele_functions.h"
-int call_cns_check_outfile(std::string prefix, std::string cns, char gap, bool call_min_depth, int min_depth){
+int call_cns_check_outfile(std::string input_id, std::string prefix, std::string cns, char gap, bool call_min_depth, int min_depth){
std::string path = "../data/test.gap.sorted.mpileup";
std::ifstream mplp(path);
- call_consensus_from_plup(mplp, prefix, 20, 0, min_depth, gap, call_min_depth);
+ call_consensus_from_plup(mplp, input_id, prefix, 20, 0, min_depth, gap, call_min_depth);
std::ifstream outFile(prefix+".fa");
std::string l;
getline(outFile, l); // Ignore first line
@@ -20,11 +20,11 @@ int main() {
std::string c_ = "CTGCTGGGTCATGGGCCCATCATGATGGTCTTGGCGATTCTAGCCTTTTTGAGATTCACGGCAATCAAGCCATCACTGGGTCTCATCAATAGATGGGGTTCAGTGGGGAAAAAAGAGGCTATGGAAACAATAAAGAAGTTCAAGAAAGAT------------------------------AGGAAGGAGAAGAAGAGACGWGGCGCAGATACTAGTGTCGGAATTGTTGGMCTCCTGCTGACCACAGCTATGGMAGCGGAGGTCACKAGACGTGGGAGTGCATACTATATGTACTTGGACWGAAACGATGCKGGGGAGGCCATATCTTTTCCAACCACATTGGGGTTGAATAAGTG";
std::string cN = "CTGCTGGGTCATGGGCCCATCATGATGGTCTTGGCGATTCTAGCCTTTTTGAGATTCACGGCAATCAAGCCATCACTGGGTCTCATCAATAGATGGGGTTCAGTGGGGAAAAAAGAGGCTATGGAAACAATAAAGAAGTTCAAGAAAGATNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNAGGAAGGAGAAGAAGAGACGWGGCGCAGATACTAGTGTCGGAATTGTTGGMCTCCTGCTGACCACAGCTATGGMAGCGGAGGTCACKAGACGTGGGAGTGCATACTATATGTACTTGGACWGAAACGATGCKGGGGAGGCCATATCTTTTCCAACCACATTGGGGTTGAATAAGTG";
std::string ck = "CTGCTGGGTCATGGGCCCATCATGATGGTCTTGGCGATTCTAGCCTTTTTGAGATTCACGGCAATCAAGCCATCACTGGGTCTCATCAATAGATGGGGTTCAGTGGGGAAAAAAGAGGCTATGGAAACAATAAAGAAGTTCAAGAAAGATAGGAAGGAGAAGAAGAGACGWGGCGCAGATACTAGTGTCGGAATTGTTGGMCTCCTGCTGACCACAGCTATGGMAGCGGAGGTCACKAGACGTGGGAGTGCATACTATATGTACTTGGACWGAAACGATGCKGGGGAGGCCATATCTTTTCCAACCACATTGGGGTTGAATAAGTG";
- num_success = call_cns_check_outfile("../data/test.gap", c_, '-', true, 0);
+ num_success = call_cns_check_outfile("", "../data/test.gap", c_, '-', true, 0);
std::cout << num_success << std::endl;
- num_success += call_cns_check_outfile("../data/test.gap", cN, 'N', true, 0);
+ num_success += call_cns_check_outfile("TESTID", "../data/test.gap", cN, 'N', true, 0);
std::cout << num_success << std::endl;
- num_success += call_cns_check_outfile("../data/test.gap", ck, 'N', false, 0);
+ num_success += call_cns_check_outfile("TESTID", "../data/test.gap", ck, 'N', false, 0);
std::cout << num_success << std::endl;
if(num_success == 0)
return 0;
=====================================
tests/test_consensus_seq_id.cpp
=====================================
@@ -0,0 +1,39 @@
+#include<iostream>
+#include<fstream>
+#include<string>
+#include "../src/call_consensus_pileup.h"
+#include "../src/allele_functions.h"
+
+int call_cns_check_outfile(std::string seq_id, std::string prefix, char gap, bool call_min_depth, int min_depth){
+ std::string path = "../data/test.gap.sorted.mpileup";
+ std::string expctd_hdr = "";
+ std::ifstream mplp(path);
+ call_consensus_from_plup(mplp, seq_id, prefix, 20, 0, min_depth, gap, call_min_depth);
+ std::ifstream outFile(prefix+".fa");
+ std::string l;
+ getline(outFile, l); // header
+ if(seq_id.empty()) {
+ char *o = new char[prefix.length() + 1];
+ strcpy(o, prefix.c_str());
+ expctd_hdr = ">Consensus_" + std::string(basename(o)) + "_threshold_0_quality_20";
+
+ } else {
+ expctd_hdr = ">" + seq_id;
+ }
+
+ return l.compare(expctd_hdr);
+}
+
+int main() {
+ int num_success = 0;
+ std::string c_ = "CTGCTGGGTCATGGGCCCATCATGATGGTCTTGGCGATTCTAGCCTTTTTGAGATTCACGGCAATCAAGCCATCACTGGGTCTCATCAATAGATGGGGTTCAGTGGGGAAAAAAGAGGCTATGGAAACAATAAAGAAGTTCAAGAAAGAT------------------------------AGGAAGGAGAAGAAGAGACGWGGCGCAGATACTAGTGTCGGAATTGTTGGMCTCCTGCTGACCACAGCTATGGMAGCGGAGGTCACKAGACGTGGGAGTGCATACTATATGTACTTGGACWGAAACGATGCKGGGGAGGCCATATCTTTTCCAACCACATTGGGGTTGAATAAGTG";
+ num_success = call_cns_check_outfile("TESTID", "../data/test.gap", '-', true, 0);
+ std::cout << num_success << std::endl;
+ num_success += call_cns_check_outfile("", "../data/test.gap", 'N', true, 0);
+ std::cout << num_success << std::endl;
+ num_success += call_cns_check_outfile("TESTID2", "../data/test.gap", 'N', false, 0);
+ std::cout << num_success << std::endl;
+ if(num_success == 0)
+ return 0;
+ return -1;
+}
\ No newline at end of file
=====================================
tests/test_getmasked.cpp
=====================================
@@ -10,7 +10,7 @@ int main()
std::ifstream masked_indices_file("../data/test.masked_primer_indices.txt");
std::string indices;
getline(masked_indices_file, indices);
- if(indices.compare("WNV_400_2_LEFT\tWNV_400_2_LEFT_alt\tWNV_400_1_LEFT\tWNV_400_1_LEFT_alt\tWNV_400_3_LEFT") == 0)
+ if(indices.compare("WNV_400_2_LEFT\tWNV_400_2_RIGHT\tWNV_400_2_LEFT_alt\tWNV_400_1_LEFT\tWNV_400_1_RIGHT\tWNV_400_1_LEFT_alt\tWNV_400_3_LEFT") == 0)
success += 1;
return (num_tests == success) ? 0 : -1;
}
=====================================
tests/test_primer_trim.cpp
=====================================
@@ -105,7 +105,7 @@ int main(){
}
// Condense cigar
std::cout << std::endl << "Condensing cigar ... " << std::endl;
- t = condense_cigar(t.cigar, t.nlength);
+ condense_cigar(&t);
replace_cigar(aln, t.nlength, t.cigar);
cigar = bam_get_cigar(aln);
print_cigar(cigar, t.nlength);
@@ -120,11 +120,19 @@ int main(){
}
}
primer_ctr++;
+ free_cigar(t);
}
ctr++;
std::cout << " ---- " << std::endl;
}
// Check if primers found at all
success = (primer_ctr > 0) ? success : -1;
+
+ bam_destroy1(aln);
+ bam_itr_destroy(iter);
+ sam_hdr_destroy(header);
+ hts_idx_destroy(idx);
+ hts_close(in);
+
return success;
}
=====================================
tests/test_primer_trim_edge_cases.cpp
=====================================
@@ -66,7 +66,7 @@ int main(){
print_cigar(cigar, aln->core.n_cigar);
// Condense cigar
std::cout << std::endl << "Condensing cigar ... " << std::endl;
- t = condense_cigar(t.cigar, t.nlength);
+ condense_cigar(&t);
replace_cigar(aln, t.nlength, t.cigar);
cigar = bam_get_cigar(aln);
print_cigar(cigar, aln->core.n_cigar);
@@ -82,8 +82,16 @@ int main(){
primer_ctr++;
ctr++;
std::cout << " ---- " << std::endl;
+ free_cigar(t);
}
// Check if primers found at all
success = (primer_ctr > 0) ? success : -1;
+
+ bam_destroy1(aln);
+ bam_itr_destroy(iter);
+ sam_hdr_destroy(header);
+ hts_idx_destroy(idx);
+ hts_close(in);
+
return success;
}
=====================================
tests/test_trim.cpp
=====================================
@@ -0,0 +1,93 @@
+#include<iostream>
+#include "../src/trim_primer_quality.h"
+#include "../src/primer_bed.h"
+#include "htslib/sam.h"
+
+typedef struct {
+ uint16_t flag;
+ uint32_t n_cigar;
+ hts_pos_t isize;
+} result_t;
+
+int test_trim(uint8_t min_qual, uint8_t sliding_window, bool no_write_flag, bool mark_qcfail_flag, int min_length, std::string testname, int nexpected, result_t *expected)
+{
+ int success = 0;
+ int res;
+
+ bam1_t *aln = bam_init1();
+
+ std::string bam = "../data/test.unmapped.sorted.bam";
+ std::string bed = "../data/test.bed";
+ std::string prefix = "/tmp/trim";
+ std::string bam_out = "/tmp/trim.bam";
+
+ std::string region_ = "";
+ std::string cmd = "@PG\tID:ivar-trim\tPN:ivar\tVN:1.0.0\tCL:ivar trim\n";
+
+ // Test with default paramaters
+ res = trim_bam_qual_primer(bam, bed, prefix, region_, min_qual, sliding_window, cmd, no_write_flag, mark_qcfail_flag, min_length);
+ if (res) {
+ success = -1;
+ std::cerr << "default paramater test failed: trim_bam_qual_primer() returned " << res << std::endl;
+ } else {
+ samFile *in = hts_open(bam_out.c_str(), "r");
+ if (!in) {
+ success = -1;
+ std::cerr << "default paramater test failed: Can't open " << bam_out << std::endl;
+ } else {
+ sam_hdr_t *hdr = sam_hdr_read(in);
+ if (!hdr) {
+ success = -1;
+ std::cerr << "default paramater test failed: Can't read header from " << bam_out << std::endl;
+ } else {
+ int n = 0;
+ while (sam_read1(in, hdr, aln) >= 0) {
+std::cerr << aln->core.n_cigar << "\t" << aln->core.l_qseq << "\t" << aln->core.mtid << "\t" << aln->core.mpos << "\t" << aln->core.isize << std::endl;
+ if (aln->core.isize != expected[n].isize) {
+ success = -1;
+ std::cerr << testname << " test failed: found isize " << aln->core.isize << " at record " << n << ": expected " << expected[n].isize << std::endl;
+ }
+ n++;
+ }
+ if (n != nexpected) {
+ success = -1;
+ std::cerr << testname << " test failed: found " << n << " records: expected " << nexpected << std::endl;
+ }
+ sam_hdr_destroy(hdr);
+ }
+ sam_close(in);
+ }
+ }
+ return success;
+}
+
+int main() {
+ int success = 0;
+
+ // Default paramaters
+ {
+ result_t expected[5] = { { 163, 2, 364 }, {163, 2, 382 }, {83, 2, -356}, {83, 2, -382}, {67, 2, -398} };
+ if (test_trim(20, 4, false, false, 30, "default paramaters", 5, expected)) success = -1;
+ }
+
+ // -e (write missing primers)
+ {
+ result_t expected[9] = { {131, 7, 403}, { 163, 2, 364 }, {163, 2, 382 }, {163, 2, 242}, {83, 2, -242}, {147, 2, -150}, {83, 2, -356}, {83, 2, -382}, {67, 2, -398} };
+ if (test_trim(20, 4, true, false, 30, "write missing primers", 9, expected)) success = -1;
+ }
+
+ // -f (mark quality trimmed as failed)
+ {
+ result_t expected[6] = { { 163, 2, 364 }, {163, 2, 382 }, {611, 2, -150}, {83, 2, -356}, {83, 2, -382}, {67, 2, -398} };
+ if (test_trim(20, 4, false, true, 30, "default paramaters", 6, expected)) success = -1;
+ }
+
+ // -e -f (write everything)
+ {
+ result_t expected[10] = { {643, 7, 403}, { 163, 2, 364 }, {163, 2, 382 }, {675, 2, 242}, {595, 2, -242}, {611, 2, -150}, {659, 2, -150}, {83, 2, -356}, {83, 2, -382}, {67, 2, -398} };
+ if (test_trim(20, 4, true, true, 30, "write everything", 10, expected)) success = -1;
+ }
+
+ return success;
+}
+
=====================================
tests/test_unpaired_trim.cpp
=====================================
@@ -30,6 +30,7 @@ int main(){
bam1_t *aln = bam_init1();
cigar_ t = {
NULL,
+ false,
0,
0
};
@@ -99,6 +100,7 @@ int main(){
get_overlapping_primers(aln, primers, overlapping_primers, true);
if(overlapping_primers.size() > 0){
cand_primer = get_min_start(overlapping_primers);
+ free_cigar(t);
t = primer_trim(aln, cand_primer.get_start() - 1, true);
replace_cigar(aln, t.nlength, t.cigar);
if(overlapping_primers.size() != rev_overlapping_primer_sizes[primer_ctr]){
@@ -128,7 +130,7 @@ int main(){
}
// Condense cigar
std::cout << std::endl << "Condensing cigar ... " << std::endl;
- t = condense_cigar(t.cigar, t.nlength);
+ condense_cigar(&t);
replace_cigar(aln, t.nlength, t.cigar);
cigar = bam_get_cigar(aln);
for (uint i = 0; i < t.nlength; ++i){
@@ -142,8 +144,16 @@ int main(){
}
}
primer_ctr++;
+ free_cigar(t);
}
// Check if primers found at all
success = (primer_ctr > 0) ? success : -1;
+
+ bam_destroy1(aln);
+ bam_itr_destroy(iter);
+ sam_hdr_destroy(header);
+ hts_idx_destroy(idx);
+ hts_close(in);
+
return success;
}
View it on GitLab: https://salsa.debian.org/med-team/ivar/-/commit/701a9ea32e3c15bf4d0fb7adeacea2451c117c82
--
View it on GitLab: https://salsa.debian.org/med-team/ivar/-/commit/701a9ea32e3c15bf4d0fb7adeacea2451c117c82
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/20201014/55747ccb/attachment-0001.html>
More information about the debian-med-commit
mailing list