pjotrp/sambamba – sambamba – Genenetwork

10 years ago

10 years ago

10 years ago

10 years ago

10 years ago

10 years ago

10 years ago

10 years ago

10 years ago

10 years ago

  1. module sambamba.flagstat;

  2. /// port of samtools flagstat tool
  3. import bamfile;
  4. import std.stdio, std.conv, std.parallelism;

  5. ulong[2] reads, pair_all, pair_good, first, second, single, pair_map, mapped,
  6. dup, diff_chr, diff_high;

  7. void writeParam(string description, ulong[2] param) {
  8. writefln("%s + %s %s", param[0], param[1], description);
  9. }

  10. float percent(ulong a, ulong b) { return to!float(a) / b * 100.0; }

  11. void writeParamWithPercentage(string description, ulong[2] param, ulong[2] total) {
  12. writefln("%s + %s %s (%.2f%%:%.2f%%)", param[0], param[1], description,
  13. percent(param[0], total[0]), percent(param[1], total[1]));
  14. }

  15. version(standalone) {
  16. int main(string[] args) {
  17. return flagstat_main(args);
  18. }
  19. }

  20. int flagstat_main(string[] args) {
  21. if (args.length == 1 || args.length > 3) {
  22. stderr.writeln("Usage: sambamba-flagstat <input.bam> [nthreads=#cores]");
  23. return 1;
  24. }

  25. try {
  26. auto threads = args.length == 2 ? totalCPUs : to!uint(args[2]);
  27. auto task_pool = new TaskPool(threads);
  28. scope(exit) task_pool.finish();

  29. auto bam = BamFile(args[1], task_pool);

  30. foreach (read; bam.alignments) {
  31. size_t failed = read.failed_quality_control ? 1 : 0;
  32. ++reads[failed];
  33. if (!read.is_unmapped) ++mapped[failed];
  34. if (read.is_duplicate) ++dup[failed];
  35. if (read.is_paired) {
  36. ++pair_all[failed];
  37. if (read.proper_pair) ++pair_good[failed];
  38. if (read.is_first_of_pair) ++first[failed];
  39. if (read.is_second_of_pair) ++second[failed];
  40. if (read.mate_is_unmapped && !read.is_unmapped) ++single[failed];
  41. if (!read.is_unmapped && !read.mate_is_unmapped) {
  42. ++pair_map[failed];
  43. if (read.ref_id != read.next_ref_id) {
  44. ++diff_chr[failed];
  45. if (read.mapping_quality >= 5)
  46. ++diff_high[failed];
  47. }
  48. }
  49. }
  50. }

  51. scope(exit) {
  52. writeParam("in total (QC-passed reads + QC-failed reads)", reads);
  53. writeParam("duplicates", dup);
  54. writeParamWithPercentage("mapped", mapped, reads);
  55. writeParam("paired in sequencing", pair_all);
  56. writeParam("read1", first);
  57. writeParam("read2", second);
  58. writeParamWithPercentage("properly paired", pair_good, pair_all);
  59. writeParam("with itself and mate mapped", pair_map);
  60. writeParamWithPercentage("singletons", single, pair_all);
  61. writeParam("with mate mapped to a different chr", diff_chr);
  62. writeParam("with mate mapped to a different chr (mapQ>=5)", diff_high);
  63. }
  64. } catch (Throwable e) {
  65. stderr.writeln(e.msg);
  66. }
  67. return 0;
  68. }

Read more here: Source link