                                   NGS ROI
                       NGS Region of Interest Analysis

                     https://www.seqan.de/apps/ngs-roi
                                Version 0.2.2
                                October, 2013


------------------------------------------------------------------------------
Table of Contents
------------------------------------------------------------------------------

  1. Overview
  2. Examples
  3. Reference and Contact
  4. Reference and Contact
  5. File Formats

------------------------------------------------------------------------------
1. Overview
------------------------------------------------------------------------------

NGS ROI is a framework for the analysis of NGS mapping data in SAM/BAM format
based on the idea of decoupling mapping information from the strict linear
ordering of the genome.

------------------------------------------------------------------------------
2. Examples
------------------------------------------------------------------------------

You can find the binaries bam2roi, roi_feature_projection, and
roi_plot_thumbnails in the directory /usr/lib/seqan/bin/ and the example files
in /usr/share/doc/seqan-apps/ngs_roi/example/.

The following examples require a Unix environment such as Linux or Mac Os X.

------------------------------------------------------------------------------
2.1. Help
------------------------------------------------------------------------------

The following program calls print the help to the user.

  $ bam2roi --help
  $ roi_feature_projection --help
  $ roi_plot_thumbnails --help

------------------------------------------------------------------------------
2.2. Converting BAM to ROI
------------------------------------------------------------------------------

The file "example.bam" consists of ~10.000 records from the 2L chromosome from
a BAM file that was obtained by mapping the SRA dataset SRR618933 agains the
D. melanogaster genome and sorting by genomic coordinate.

We can generate a ROI file from this BAM file:

  $ bam2roi -if example.bam -of example.roi
  ...

  $ head -5 example.roi
  #ROI 0.3
  ##ref	start	end	name	length	strand	max_count	num_reads	gc_content	counts
  2L	9	60	region_0	52	+	3	4	0.4423	1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,...
  2L	62	138	region_1	77	+	15	20	0.3896	2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,8,9,...
  2L	143	596	region_2	454	+	18	148	0.4009	2,3,3,3,3,3,3,3,3,3,3,3,4,4,4,8,8,...
  2L	599	1252	region_3	654	+	21	207	0.4083	1,1,1,1,1,1,1,1,1,1,2,2,3,3,5,...
  2L	1255	2352	region_4	1098	+	26	340	0.4071	2,2,2,3,3,3,3,3,3,3,3,...
  2L	2357	2807	region_5	451	+	21	141	0.4058	1,1,1,1,1,1,1,1,1,1,6,6,8,...
  2L	2812	2888	region_6	77	+	9	16	0.3896	1,1,1,1,1,1,1,1,1,1,1,1,1,...
  2L	2890	3344	region_7	455	+	21	121	0.4110	1,1,1,1,1,2,2,2,2,3,3,5,5,...

------------------------------------------------------------------------------
2.2. Converting BAM to ROI
------------------------------------------------------------------------------

We can create a PNG thumbnail view of the ROI file:

  $ roi_plot_thumbnails -if example.roi -o plots
  ...

  $ ls plots*
  plots0.png

------------------------------------------------------------------------------
2.3. Projecting to BED
------------------------------------------------------------------------------

The file "dmel.bed" contains a subset of the gene annotations for the
D. melanogaster genome.  We had to translate some of the reference names
(e.g. chr2L to 2L) such that they match the reference name from the BAM and
ROI file.

The program roi_feature_projection reads a ROI and a BED file and create one
output ROI record for each BED record in the input.  The counts array will
have the length of the BED record and the counts of the ROI records
overlapping with the BED record will be projected to the counts of the new ROI
record.

The ROI file and the BED file have be sorted by coordinate (lexicographically
reference names, then by begin and end position).  You will find the script
bed_sort.sh in /usr/lib/seqan/bin/.  We first use it to sort the BED file
appropriately.

  $ head -5 dmel.bed
  2L	8384466	8388803	CG13384-RD	0	+
  2L	8384466	8388803	CG13384-RC	0	+
  2L	8384466	8388803	CG13384-RE	0	+
  2L	8384466	8388803	CG13384-RG	0	+
  2L	8384466	8388803	CG13384-RF	0	+

  $ bed_sort.sh -i dmel.bed -o dmel.sorted.bed

  $ head -5 dmel.sorted.bed
  2L	7528	9491	CG11023-RA	0	+
  2L	9835	18583	CG2671-RA	0	-
  2L	9835	18583	CG2671-RC	0	-
  2L	9835	21372	CG2671-RB	0	-
  2L	9835	21372	CG2671-RD	0	-

The ROI file will be sorted appropriately already after creation from BAM.
However, to demonstrate the script roi_sort.sh, here is how to sort the ROI
file by coordinate.

  $ sort_roi.sh -p -i exampe.roi -o example.sorted.roi

Then, we use the program roi_feature_projection to project the ROI records to
the BED records.

  $ roi_feature_projection -ir example.sorted.roi -if dmel.sorted.bed \
      -or projected_bed.roi
  ...

  $ head -5 projected_bed.roi
  ##ref   begin_pos       end_pos region_name     length  strand  max_count       counts
  2L      7528    9491    CG11023-RA_0    1964    +       0       0,0,0,0,0,0,0,0,0,0,...
  2L      9835    18583   CG2671-RA_0     8749    -       0       0,0,0,0,0,0,0,0,0,0,...
  2L      9835    18583   CG2671-RC_0     8749    -       55      1,1,1,1,1,1,1,1,1,1,...
  2L      9835    21372   CG2671-RB_0     11538   -       55      1,1,1,1,1,1,1,1,1,1,...

------------------------------------------------------------------------------
2.4. Projecting to GFF/GTF
------------------------------------------------------------------------------

Next, we demonstrate how to project the ROI file to a GTF file.  The procedure
is the same for GFF files.

First, we sort the GTF file by coordinate.

  $ head -5 dmel.gtf
  2L	dm3_flyBaseGene	exon	8384467	8384895	0.000000	+	.	gene_id "CG13384-RD"; transcript_id "CG13384-RD";
  2L	dm3_flyBaseGene	exon	8384997	8385097	0.000000	+	.	gene_id "CG13384-RD"; transcript_id "CG13384-RD";
  2L	dm3_flyBaseGene	start_codon	8386247	8386249	0.000000	+	.	gene_id "CG13384-RD"; transcript_id "CG13384-RD";
  2L	dm3_flyBaseGene	CDS	8386247	8386385	0.000000	+	0	gene_id "CG13384-RD"; transcript_id "CG13384-RD";
  2L	dm3_flyBaseGene	exon	8386217	8386385	0.000000	+	.	gene_id "CG13384-RD"; transcript_id "CG13384-RD";

  $ gff_sort.sh -i dmel.gtf -o dmel.sorted.gtf

  $ head -5 dmel.sorted.gtf
  2L	dm3_flyBaseGene	exon	7529	8116	0.000000	+	.	gene_id "CG11023-RA"; transcript_id "CG11023-RA";
  2L	dm3_flyBaseGene	start_codon	7680	7682	0.000000	+	.	gene_id "CG11023-RA"; transcript_id "CG11023-RA";
  2L	dm3_flyBaseGene	CDS	7680	8116	0.000000	+	0	gene_id "CG11023-RA"; transcript_id "CG11023-RA";
  2L	dm3_flyBaseGene	CDS	8229	8589	0.000000	+	1	gene_id "CG11023-RA"; transcript_id "CG11023-RA";
  2L	dm3_flyBaseGene	exon	8229	8589	0.000000	+	.	gene_id "CG11023-RA"; transcript_id "CG11023-RA";

Then, we project the ROI file to the GTF file, splicing together the exons
having the same transcript id.

  $ roi_feature_projection -ir example.roi -if sorted.gtf \
      -or projected_gtf.roi --gff-type exon --gff-group-by transcript_id
  ...
  $ head -5 projected_gtf.roi
  ##ref   begin_pos       end_pos region_name     length  strand  max_count       counts
  2L      7529    9491    CG11023-RA      1773    +       34      8,9,13,12,12,12,12,12,11...
  2L      9836    18583   CG2671-RA       5407    +       55      1,1,1,1,1,1,1,1,1,1,0,0,...
  2L      9836    18583   CG2671-RC       5264    -       55      1,1,1,1,1,1,1,1,1,1,0,0,...
  2L      9836    21372   CG2671-RB       5154    -       55      1,1,1,1,1,1,1,1,1,1,0,0,...
  2L      9836    21372   CG2671-RD       5229    -       55      1,1,1,1,1,1,1,1,1,1,0,0,...

------------------------------------------------------------------------------
2.5. Sorting and Plotting
------------------------------------------------------------------------------

We can now sort the ROI file by the the highest count (7th field, descending):

  $ roi_sort.sh -n 7 -r -i projected_gtf.roi -o projected_gtf.by_max_count.roi

  $ head -5 projected_gtf.by_max_count.roi
  ##ref   begin_pos       end_pos region_name     length  strand  max_count       counts
  2L      72388   76203   CG11372-RA      1832    +       69      0,0,0,1,1,1,1,1,1,1,1,1,1,..
  2L      9836    18583   CG2671-RA       5407    +       55      1,1,1,1,1,1,1,1,1,1,0,0,0,..
  2L      9836    18583   CG2671-RC       5264    -       55      1,1,1,1,1,1,1,1,1,1,0,0,0,..
  2L      9836    21372   CG2671-RB       5154    -       55      1,1,1,1,1,1,1,1,1,1,0,0,0,..

Let us extract the records that do not have a max_count of 0.

  $ awk '($7 > 0) { print $0) }' projected_gtf.by_max_count.roi > to_plot.roi
  $ roi_plot_9.sh -i to_plot.roi -o out.pdf

Note that the generated PDF file already contains links to a local IGV
instance listening on port 60151.  If you haven't already, you should start an
IGV instance, load the drosophila reference genome, add some feature tracks
and enable remote control via HTTP in IGV.  Clicking on a plot in the PDF file
will open this locus in IGV.

------------------------------------------------------------------------------
2.6. Computing Metrics
------------------------------------------------------------------------------

We can compute some metrics using the script roi_metrics.Rscript.  This script
can also serve as a starting point for your own scripts:

  $ roi_metrics.Rscript -i to_plot.roi -o metrics.roi

  $ head -3 metrics.roi
  # ROI written from R
  ##ref	begin_pos	end_pos	region_name	length	strand	max_count	min	median	mean	quantile75	quantile95	aoc	xreaXX	r3linear	distMax3p	distMax5p	counts
  2L	72388	76203	CG11372-RA	1832	+	69	0	7	11.8733624454148	19	44	0.172077716600215	0.328195209163977	0.00369582266667077	1286	1286	0,0,...

Let us now sort descendingly by mean coverage.

  $ roi_sort.sh -n 10 -r metrics.roi -o metrics.by_mean.roi

And we can also plot this:

  $ roi_plot_9.sh -i metric.by_mean.roi -o out.pdf

------------------------------------------------------------------------------
3. Reference and Contact
------------------------------------------------------------------------------

Reference

  Holtgrewe M., Coppee, J.-Y., Reinert, K., Jagla, B.  Novel Post-Alignment
  Visualization and Characterization of High-Throughput Sequencing
  Experiments.  Unpublished.

Contact

  Bernd Jagla <bernd.jagla@pasteur.fr>

------------------------------------------------------------------------------
4. Galaxy Integration
------------------------------------------------------------------------------

The integration of ngs_roi is available in the Galaxy Tool Shed (Testing)

  https://testtoolshed.g2.bx.psu.edu/view/holtgrewe/ngs_roi

------------------------------------------------------------------------------
5. File Formats
------------------------------------------------------------------------------

------------------------------------------------------------------------------
5.1. ROI Format
------------------------------------------------------------------------------

The file begins with a sequence of comments, followed by one column header,
followed by the records.

Comments start with a single hash (#) and are ignored in analysis.

The column header starts with two hashes (##) and is followed by the tab
separated column names. The column names must not contain spaces.

The columns 1-7 are fixed. The last column is always "counts". The columns are
as follows

  1. ref         -- The reference name.
  2. begin_pos   -- The begin position.
  3. end_pos     -- The end position.
  4. region_name -- The name of the region.
  5. length      -- The length of the region.
  6. strand      -- The strand, one of "+" and "-". "+" in case of being not
                    strand-specific.
  7. max_count   -- The largest value of the counts columns.

  additional columns
     Additional annotation data for the ROI.

  N. counts       -- Comma-separated list of length unsignd integers.

------------------------------------------------------------------------------
5.2. BED, GFF and GTF format
------------------------------------------------------------------------------

See the following external resources for a specification of these formats:

  BED https://genome.ucsc.edu/FAQ/FAQformat.html#format1
  GFF https://www.sanger.ac.uk/resources/software/gff/spec.html
      https://www.sequenceontology.org/gff3.shtml
  GTF https://mblab.wustl.edu/GTF22.html
