iVar
ivar.cpp
Go to the documentation of this file.
1 
3 #include <iostream>
4 #include <unistd.h>
5 
6 #include "remove_reads_from_amplicon.h"
7 #include "call_consensus_pileup.h"
8 #include "call_variants.h"
9 #include "trim_primer_quality.h"
10 #include "get_masked_amplicons.h"
11 #include "suffix_tree.h"
12 
13 const std::string VERSION = "1.0";
14 
15 struct args_t {
16  std::string bam; // -i
17  std::string bed; // -b
18  std::string prefix; // -p
19  std::string ref; // -r
20  std::string region; // -R
21  uint8_t min_qual; // -q
22  uint8_t sliding_window; // -s
23  double min_threshold; // -t
24  int min_length; // -m
25  std::string f1; // -1
26  std::string f2; // -2
27  std::string adp_path; // -a
28 } g_args;
29 
30 void print_usage(){
31  std::cout <<
32  "Usage: ivar [command <trim|callvariants|filtervariants|consensus|createbed|getmasked|removereads|version|help>]\n"
33  "\n"
34  " Command Description\n"
35  " trim Trim reads in aligned BAMs\n"
36  " variants Call variants from aligned BAM file\n"
37  " filtervariants Filter variants across replicates\n"
38  " consensus Call consensus from aligned BAM file\n"
39  " getmasked Get amplicons with primer mismatches\n"
40  " removereads Remove reads from trimmed BAM file\n"
41  " version Show version information\n"
42  " trimadapter (EXPERIMENTAL) Trim adapter sequences from reads\n"
43  "\n"
44  "To view detailed usage for each command type `ivar <command>` \n";
45 }
46 
47 void print_trim_usage(){
48  std::cout <<
49  "Usage: ivar trim -i <input.bam> -b <primers.bed> -p <prefix> [-m <min-length>] [-q <min-quality>] [-s <sliding-window-width>]\n\n"
50  "Input Options Description\n"
51  " -i (Required) Indexed aligned bam file to trim primers and quality\n"
52  " -b (Required) BED file with primer sequences and positions\n"
53  " -m Minimum length of read to retain after trimming (Default: 30)\n"
54  " -q Minimum quality threshold for sliding window to pass (Default: 20)\n"
55  " -s Width of sliding window (Default: 4)\n\n"
56  "Output Options Description\n"
57  " -p (Required) Prefix for the output BAM file\n";
58 }
59 
60 void print_variants_usage(){
61  std::cout <<
62  "Usage: samtools mpileup -A -d 300000 --reference <reference-fasta> -Q 0 -F 0 <input.bam> | ivar variants -p <prefix> [-q <min-quality>] [-t <min-frequency-threshold>]\n\n"
63  "Note : samtools mpileup output must be piped into ivar variants\n\n"
64  "Input Options Description\n"
65  " -q Minimum quality score threshold to count base (Default: 20)\n"
66  " -t Minimum frequency threshold(0 - 1) to call variants (Default: 0.03)\n\n"
67  "Output Options Description\n"
68  " -p (Required) Prefix for the output tsv variant file\n\n";
69 }
70 
71 void print_filtervariants_usage(){
72  std::cout <<
73  "Usage: ivar filtervariants -p <prefix> replicate-one.tsv replicate-two.tsv ... \n\n"
74  "Input: Variant tsv files for each replicate\n\n"
75  "Output Options Description\n"
76  " -p (Required) Prefix for the output filtered tsv file\n";
77 }
78 
79 void print_consensus_usage(){
80  std::cout <<
81  "Usage: samtools mpileup -A -d 300000 -Q 0 -F 0 <input.bam> | ivar consensus -p <prefix> \n\n"
82  "Note : samtools mpileup output must be piped into `ivar consensus`\n\n"
83  "Input Options Description\n"
84  " -q Minimum quality score threshold to count base (Default: 20)\n"
85  " -t Minimum frequency threshold(0 - 1) to call consensus. (Default: 0)\n"
86  " Frequently used thresholds | Description\n"
87  " ---------------------------|------------\n"
88  " 0 | Majority or most common base\n"
89  " 0.2 | Bases that make up atleast 20% of the depth at a position\n"
90  " 0.5 | Strict or bases that make up atleast 50% of the depth at a position\n"
91  " 0.9 | Strict or bases that make up atleast 90% of the depth at a position\n"
92  " 1 | Identical or bases that make up 100% of the depth at a position. Will have highest ambiguities\n"
93  "Output Options Description\n"
94  " -p (Required) Prefix for the output fasta file and quality file\n";
95 }
96 
97 void print_removereads_usage(){
98  std::cout <<
99  "Usage: ivar removereads -i <input.trimmed.bam> -p <prefix> primer-index-1 primer-index-2 primer-index-3 primer-index-4 ... \n\n"
100  "Input Options Description\n"
101  " -i (Required) Input BAM file trimmed with ‘ivar trim’. Must be sorted and indexed, which can be done using sort_index_bam.sh\n"
102  "Output Options Description\n"
103  " -p (Required) Prefix for the output filtered BAM file\n";
104 }
105 
106 void print_getmasked_usage(){
107  std::cout <<
108  "Usage: ivar getmasked -i <input-filtered.tsv> -b <primers.bed>\n\n"
109  "Input Options Description\n"
110  " -i (Required) Input filtered variants tsv generated from `ivar filtervariants`\n"
111  " -b (Required) BED file with primer sequences and positions\n";
112 }
113 
114 void print_trimadapter_usage(){
115  std::cout <<
116  "Usage: ivar trimadapter [-f1 <input-fastq>] [-f2 <input-fastq-2>] [-p prefix] [-a <adapter-fasta-file>]\n\n"
117  "Input Options Description\n"
118  " -1 (Required) Input fastq file\n"
119  " -2 Input fastq file 2 (for pair ended reads)\n"
120  " -a (Required) Adapter Fasta File\n"
121  " -p (Required) Prefix of output fastq files\n";
122 }
123 
124 void print_version_info(){
125  std::cout << "iVar version " << VERSION << std::endl <<
126  "\nPlease raise issues and bug reports at https://github.com/andersen-lab/ivar/\n";
127 }
128 
129 static const char *trim_opt_str = "i:b:p:m::q::s::h?";
130 static const char *variants_opt_str = "p:t::q::h?";
131 static const char *consensus_opt_str = "p:q::t::h?";
132 static const char *removereads_opt_str = "i:p:h?";
133 static const char *filtervariants_opt_str = "p:h?";
134 static const char *getmasked_opt_str = "i:b:h?";
135 static const char *trimadapter_opt_str = "1:2::p:a:h?";
136 
144 int main(int argc, char* argv[]){
145  if(argc == 1){
146  print_usage();
147  return -1;
148  }
149  std::string cmd(argv[1]);
150  if(cmd.compare("-v") == 0){
151 
152  return 0;
153  }
154  int opt = 0, res = 0;
155  // Sift arg by 1 for getopt
156  argv[1] = argv[0];
157  argv++;
158  argc--;
159  g_args.min_qual = 255;
160  g_args.sliding_window = 255;
161  g_args.min_threshold = -1;
162  g_args.min_length = -1;
163  if (cmd.compare("trim") == 0){
164  opt = getopt( argc, argv, trim_opt_str);
165  while( opt != -1 ) {
166  switch( opt ) {
167  case 'i':
168  g_args.bam = optarg;
169  break;
170  case 'b':
171  g_args.bed = optarg;
172  break;
173  case 'p':
174  g_args.prefix = optarg;
175  break;
176  case 'm':
177  g_args.min_length = atoi(optarg);
178  break;
179  case 'q':
180  g_args.min_qual = atoi(optarg);
181  break;
182  case 's':
183  g_args.sliding_window = atoi(optarg);
184  break;
185  case 'h':
186  case '?':
187  print_trim_usage();
188  return -1;
189  break;
190  }
191  opt = getopt( argc, argv, trim_opt_str);
192  }
193  if(g_args.bam.empty() || g_args.bed.empty() || g_args.prefix.empty()){
194  print_trim_usage();
195  return -1;
196  }
197  g_args.min_qual = (g_args.min_qual == 255) ? 20 : g_args.min_qual;
198  g_args.sliding_window = (g_args.sliding_window == 255) ? 4 : g_args.sliding_window;
199  g_args.min_length = (g_args.min_length == -1) ? 30 : g_args.min_length;
200  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, g_args.min_length);
201  } else if (cmd.compare("variants") == 0){
202  opt = getopt( argc, argv, variants_opt_str);
203  while( opt != -1 ) {
204  switch( opt ) {
205  case 'p':
206  g_args.prefix = optarg;
207  break;
208  case 't':
209  g_args.min_threshold = atof(optarg);
210  break;
211  case 'q':
212  g_args.min_qual = atoi(optarg);
213  break;
214  case 'h':
215  case '?':
216  print_variants_usage();
217  return 0;
218  break;
219  }
220  opt = getopt( argc, argv, variants_opt_str);
221  }
222  if(g_args.prefix.empty()){
223  print_variants_usage();
224  return -1;
225  }
226  g_args.min_threshold = (g_args.min_threshold == -1 || g_args.min_threshold < 0 || g_args.min_threshold >= 1) ? 0.03: g_args.min_threshold;
227  g_args.min_qual = (g_args.min_qual == 255) ? 20 : g_args.min_qual;
228  res = call_variants_from_plup(std::cin, g_args.prefix, g_args.min_qual, g_args.min_threshold);
229  } else if (cmd.compare("consensus") == 0){
230  opt = getopt( argc, argv, consensus_opt_str);
231  g_args.min_threshold = 0;
232  while( opt != -1 ) {
233  switch( opt ) {
234  case 't':
235  g_args.min_threshold = atof(optarg);
236  break;
237  case 'p':
238  g_args.prefix = optarg;
239  break;
240  case 'q':
241  g_args.min_qual = atoi(optarg);
242  break;
243  case 'h':
244  case '?':
245  print_consensus_usage();
246  return 0;
247  break;
248  }
249  opt = getopt( argc, argv, consensus_opt_str);
250  }
251  if(g_args.prefix.empty()){
252  print_consensus_usage();
253  return -1;
254  }
255  g_args.min_qual = (g_args.min_qual == 255) ? 20 : g_args.min_qual;
256  std::cout <<"Min Quality" << g_args.min_qual << std::endl;
257  std::cout << "Threshold: " << g_args.min_threshold << std::endl;
258  res = call_consensus_from_plup(std::cin, g_args.prefix, g_args.min_qual, g_args.min_threshold);
259  } else if (cmd.compare("removereads") == 0){
260  opt = getopt( argc, argv, removereads_opt_str);
261  while( opt != -1 ) {
262  switch( opt ) {
263  case 'i':
264  g_args.bam = optarg;
265  break;
266  case 'p':
267  g_args.prefix = optarg;
268  break;
269  case 'h':
270  case '?':
271  print_removereads_usage();
272  return 0;
273  break;
274  }
275  opt = getopt( argc, argv, removereads_opt_str);
276  }
277  if(optind >= argc){
278  print_removereads_usage();
279  return -1;
280  }
281  uint16_t amp[100];
282  for(int i = optind; i<argc;i++){
283  amp[i] = atoi(argv[i]);
284  std::cout << amp[i];
285  }
286  if(g_args.bam.empty() || g_args.prefix.empty()){
287  print_removereads_usage();
288  return -1;
289  }
290  res = rmv_reads_from_amplicon(g_args.bam, g_args.region, g_args.prefix, amp, argc - optind);
291  } else if(cmd.compare("filtervariants") == 0){
292  opt = getopt( argc, argv, filtervariants_opt_str);
293  while( opt != -1 ) {
294  switch( opt ) {
295  case 'p':
296  g_args.prefix = optarg;
297  break;
298  case 'h':
299  case '?':
300  print_filtervariants_usage();
301  return 0;
302  break;
303  }
304  opt = getopt( argc, argv, filtervariants_opt_str);
305  }
306  if(optind >= argc){
307  print_filtervariants_usage();
308  return -1;
309  }
310  if(g_args.prefix.empty()){
311  print_filtervariants_usage();
312  return -1;
313  }
314  std::string rep = "get_common_variants.sh ";
315  for(int i = optind; i<argc;i++){
316  rep += argv[i];
317  rep += " ";
318  }
319  rep += " | sort -s -n -k 2 > "+g_args.prefix+".tsv";
320  system(rep.c_str());
321  } else if(cmd.compare("getmasked") == 0){
322  opt = getopt( argc, argv, getmasked_opt_str);
323  while( opt != -1 ) {
324  switch( opt ) {
325  case 'i':
326  g_args.bam = optarg;
327  break;
328  case 'b':
329  g_args.bed = optarg;
330  break;
331  case 'h':
332  case '?':
333  print_getmasked_usage();
334  return 0;
335  break;
336  }
337  opt = getopt( argc, argv, getmasked_opt_str);
338  }
339  if(g_args.bed.empty() || g_args.bam.empty()){
340  print_getmasked_usage();
341  return -1;
342  }
343  res = get_primers_with_mismatches(g_args.bed, g_args.bam);
344  } else if (cmd.compare("trimadapter") == 0){
345  opt = getopt( argc, argv, trimadapter_opt_str);
346  while( opt != -1 ) {
347  switch( opt ) {
348  case '1':
349  g_args.f1 = optarg;
350  break;
351  case '2':
352  g_args.f2 = optarg;
353  break;
354  case 'p':
355  g_args.prefix = optarg;
356  break;
357  case 'a':
358  g_args.adp_path = optarg;
359  break;
360  case 'h':
361  case '?':
362  print_trimadapter_usage();
363  return 0;
364  break;
365  }
366  opt = getopt( argc, argv, trimadapter_opt_str);
367  }
368  if(g_args.f1.empty() || g_args.prefix.empty() || g_args.adp_path.empty()){
369  print_trimadapter_usage();
370  return -1;
371  }
372  res = trim_adapter(g_args.f1, g_args.f2, g_args.adp_path, g_args.prefix);
373  } else {
374  print_usage();
375  }
376  return res;
377 }
args_t
Definition: ivar.cpp:15
main
int main(int argc, char *argv[])
Definition: ivar.cpp:144
Generated by   doxygen 1.8.14