annotate bismark_methylation_extractor @ 5:b100248c35b8 draft

Uploaded

1 #!/usr/bin/perl
2 use warnings;
3 use strict;
4 $|++;
5 use Getopt::Long;
6 use Cwd;
7 use Carp;
8 use FindBin qw($Bin);
9 use lib "$Bin/../lib";
10 
11 
12 ## This program is Copyright (C) 2010-13, Felix Krueger ([e-mail])
13 
14 ## This program is free software: you can redistribute it and/or modify
15 ## it under the terms of the GNU General Public License as published by
16 ## the Free Software Foundation, either version 3 of the License, or
17 ## (at your option) any later version.
18 
19 ## This program is distributed in the hope that it will be useful,
20 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
21 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
22 ## GNU General Public License for more details.
23 
24 ## You should have received a copy of the GNU General Public License
25 ## along with this program. If not, see <http://www.gnu.org/licenses/>.
26 
27 my @filenames; # input files
28 my %counting;
29 my $parent_dir = getcwd();
30 
31 my %fhs;
32 
33 my $version = 'v0.10.1';
34 my ($ignore,$genomic_fasta,$single,$paired,$full,$report,$no_overlap,$merge_non_CpG,$vanilla,$output_dir,$no_header,$bedGraph,$remove,$coverage_threshold,$counts,$cytosine_report,$genome_folder,$zero,$CpG_only,$CX_context,$split_by_chromosome,$sort_size,$samtools_path,$gzip,$ignore_r2,$mbias_only,$gazillion,$ample_mem) = process_commandline();
35 
36 
37 ### only needed for bedGraph output
38 my @sorting_files; # if files are to be written to bedGraph format, these are the methylation extractor output files
39 my @methylcalls = qw (0 0 0); # [0] = methylated, [1] = unmethylated, [2] = total
40 my @bedfiles;
41 
42 ### only needed for genome-wide cytosine methylation report
43 my %chromosomes;
44 
45 my %mbias_1;
46 my %mbias_2;
47 
48 ##############################################################################################
49 ### Summarising Run Parameters
50 ##############################################################################################
51 
52 ### METHYLATION EXTRACTOR
53 
54 warn "Summarising Bismark methylation extractor parameters:\n";
55 warn '='x63,"\n";
56 
57 if ($single){
58   if ($vanilla){
59     warn "Bismark single-end vanilla format specified\n";
60   }
61   else{
62     warn "Bismark single-end SAM format specified (default)\n"; # default
63   }
64 }
65 elsif ($paired){
66   if ($vanilla){
67     warn "Bismark paired-end vanilla format specified\n";
68   }
69   else{
70     warn "Bismark paired-end SAM format specified (default)\n"; # default
71   }
72 }
73 
74 if ($single){
75   if ($ignore){
76     warn "First $ignore bp will be disregarded when processing the methylation call string\n";
77   }
78 }
79 else{ ## paired-end
80   if ($ignore){
81     warn "First $ignore bp will be disregarded when processing the methylation call string of Read 1\n";
82   }
83   if ($ignore_r2){
84     warn "First $ignore_r2 bp will be disregarded when processing the methylation call string of Read 2\n";
85   }
86 }
87 
88 
89 if ($full){
90   warn "Strand-specific outputs will be skipped. Separate output files for cytosines in CpG, CHG and CHH context will be generated\n";
91 }
92 if ($merge_non_CpG){
93   warn "Merge CHG and CHH context to non-CpG context specified\n";
94 }
95 ### output directory
96 if ($output_dir eq ''){
97   warn "Output will be written to the current directory ('$parent_dir')\n";
98 }
99 else{
00   warn "Output path specified as: $output_dir\n";
01 }
02 
03 
04 sleep (1);
05 
06 ### BEDGRAPH
07 
08 if ($bedGraph){
09   warn "\n\nSummarising bedGraph parameters:\n";
10   warn '='x63,"\n";
11 
12   if ($counts){
13     warn "Generating additional output in bedGraph and coverage format\nbedGraph format:\t<Chromosome> <Start Position> <End Position> <Methylation Percentage>\ncoverage format:\t<Chromosome> <Start Position> <End Position> <Methylation Percentage> <count methylated> <count non-methylated>\n\n";
14   }
15   else{
16     warn "Generating additional sorted output in bedGraph format (output format: <Chromosome> <Start Position> <End Position> <Methylation Percentage>)\n";
17   }
18 
19   warn "Using a cutoff of $coverage_threshold read(s) to report cytosine positions\n";
20 
21   if ($CX_context){
22     warn "Reporting and sorting methylation information for all cytosine context (sorting may take a long time, you have been warned ...)\n";
23   }
24   else{ # default
25     $CpG_only = 1;
26     warn "Reporting and sorting cytosine methylation information in CpG context only (default)\n";
27   }
28 
29   if ($remove){
30     warn "White spaces in read ID names will be removed prior to sorting\n";
31   }
32 
33   if ($ample_mem){
34     warn "Sorting chromosomal postions for the bedGraph step using arrays instead of using UNIX sort\n";
35   }
36   elsif (defined $sort_size){
37     warn "The bedGraph UNIX sort command will use the following memory setting:\t'$sort_size'. Temporary directory used for sorting is the output directory\n";
38   }
39   else{
40     warn "Setting a default memory usage for the bedGraph UNIX sort command to 2GB\n";
41   }
42 
43 
44 
45   sleep (1);
46 
47   if ($cytosine_report){
48     warn "\n\nSummarising genome-wide cytosine methylation report parameters:\n";
49     warn '='x63,"\n";
50     warn "Generating comprehensive genome-wide cytosine report\n(output format: <Chromosome> <Position> <Strand> <count methylated> <count non-methylated>  <C-context>  <trinucleotide context> )\n";
51 
52 
53     if ($CX_context){
54       warn "Reporting methylation for all cytosine contexts. Be aware that this will generate enormous files\n";
55     }
56     else{ # default
57       $CpG_only = 1;
58       warn "Reporting cytosine methylation in CpG context only (default)\n";
59     }
60 
61     if ($split_by_chromosome){
62       warn "Splitting the cytosine report output up into individual files for each chromosome\n";
63     }
64 
65     ### Zero-based coordinates
66     if ($zero){
67       warn "Using zero-based genomic coordinates (user-defined)\n";
68     }
69     else{ # default, 1-based coords
70       warn "Using 1-based genomic coordinates (default)\n";
71     }
72 
73     ### GENOME folder
74     if ($genome_folder){
75       unless ($genome_folder =~/\/$/){
76 	$genome_folder =~ s/$/\//;
77       }
78       warn "Genome folder was specified as $genome_folder\n";
79     }
80     else{
81       $genome_folder  = '/data/public/Genomes/Mouse/NCBIM37/';
82       warn "Using the default genome folder /data/public/Genomes/Mouse/NCBIM37/\n";
83     }
84     sleep (1);
85   }
86 }
87 
88 warn "\n";
89 sleep (5);
90 
91 ######################################################
92 ### PROCESSING FILES
93 ######################################################
94 
95 foreach my $filename (@filenames){
96   # resetting counters and filehandles
97   %fhs = ();
98   %counting =(
99 	      total_meCHG_count => 0,
00 	      total_meCHH_count => 0,
01 	      total_meCpG_count => 0,
02 	      total_unmethylated_CHG_count => 0,
03 	      total_unmethylated_CHH_count => 0,
04 	      total_unmethylated_CpG_count => 0,
05 	      sequences_count => 0,
06 	     );
07 
08   @sorting_files = ();
09   @bedfiles = ();
10 
11   %mbias_1 = ();
12   %mbias_2 = ();
13 
14   ### performing a quick check to see if a paired-end SAM file has been sorted by positions which does interfere with the logic used by the extractor
15   unless ($vanilla){
16     if ($paired){
17       test_positional_sorting($filename);
18     }
19   }
20 
21   process_Bismark_results_file($filename);
22 
23   ### Closing all filehandles so that the Bismark methylation extractor output doesn't get truncated due to buffering issues
24   foreach my $fh (keys %fhs) {
25     if ($fh =~ /^[1230]$/) {
26       foreach my $context (keys %{$fhs{$fh}}) {
27 	close $fhs{$fh}->{$context} or die $!;
28       }
29     }
30     else{
31       close $fhs{$fh} or die $!;
32     }
33   }
34 
35   ### printing out all M-Bias data
36   produce_mbias_plots ($filename);
37 
38   delete_unused_files();
39 
40   if ($bedGraph){
41 
42     my $out = (split (/\//,$filename))[-1]; # extracting the filename if a full path was specified
43     $out =~ s/gz$//;
44     $out =~ s/sam$//;
45     $out =~ s/bam$//;
46     $out =~ s/txt$//;
47     $out =~ s/$/bedGraph/;
48 
49     my $bedGraph_output = $out;
50     my @args;
51 
52     if ($remove){
53       push @args, '--remove';
54     }
55     if ($CX_context){
56       push @args, '--CX_context';
57     }
58     if ($no_header){
59       push @args, '--no_header';
60     }
61     if ($gazillion){
62       push @args, '--gazillion';
63     }
64     if ($ample_mem){
65       push @args, '--ample_memory';
66     }
67 
68 
69     #   if ($counts){
70     #      push @args, "--counts";
71     #   }
72 
73     push @args, "--buffer_size $sort_size";
74     push @args, "--cutoff $coverage_threshold";
75     push @args, "--output $bedGraph_output";
76     push @args, "--dir '$output_dir'";
77 
78     ### adding all files to be sorted to @args
79     foreach my $f (@sorting_files){
80       push @args, $f;
81     }
82 
83     #  print join "\t",@args,"\n";
84 
85     system ("$Bin/bismark2bedGraph @args");
86 
87     warn "Finished BedGraph conversion ...\n\n";
88     sleep(3);
89 
90     # open (OUT,'>',$output_dir.$bedGraph_output) or die "Problems with the bedGraph output filename detected: file path: '$output_dir'\tfile name: '$bedGraph_output' $!";
91     # warn "Writing bedGraph to file: $bedGraph_output\n";
92     # process_bedGraph_output();
93     # close OUT or die $!;
94 
95     ### genome-wide cytosine methylation report requires bedGraph processing anyway
96     if ($cytosine_report){
97 
98       @args = (); # resetting @args
99       my $cytosine_out = $out;
00       $cytosine_out =~ s/bedGraph$//;
01 
02       if ($CX_context){
03 	$cytosine_out =~ s/$/CX_report.txt/;
04       }
05       else{
06 	$cytosine_out =~ s/$/CpG_report.txt/;
07       }
08 
09       push @args, "--output $cytosine_out";
10       push @args, "--dir '$output_dir'";
11       push @args, "--genome '$genome_folder'";
12       push @args, "--parent_dir '$parent_dir'";
13 
14       if ($zero){
15 	push @args, "--zero";
16       }
17       if ($CX_context){
18 	push @args, '--CX_context';
19       }
20       if ($split_by_chromosome){
21 	push @args, '--split_by_chromosome';
22       }
23 
24       my $coverage_output = $bedGraph_output;
25       $coverage_output =~ s/bedGraph$/bismark.cov/;
26 
27       push @args, $output_dir . $coverage_output; # this will be the infile
28 
29       system ("$Bin/coverage2cytosine @args");
30       # generate_genome_wide_cytosine_report($bedGraph_output,$cytosine_out);
31       warn "\n\nFinished generating genome-wide cytosine report\n\n";
32     }
33   }
34 }
35 
36 sub delete_unused_files{
37 
38   warn "Deleting unused files ...\n\n"; sleep(1);
39 
40   my $index = 0;
41 
42   while ($index <= $#sorting_files){
43     if ($sorting_files[$index] =~ /gz$/){
44       open (USED,"zcat $sorting_files[$index] |") or die "Failed to read from methylation extractor output file $sorting_files[$index]: $!\n";
45     }
46     else{
47       open (USED,$sorting_files[$index]) or die "Failed to read from methylation extractor output file $sorting_files[$index]: $!\n";
48     }
49 
50     my $used = 0;
51 
52     while (<USED>){
53       next if (/^Bismark/);
54       if ($_){
55 	$used = 1;
56 	last;
57       }
58     }
59 
60     if ($used){
61       warn "$sorting_files[$index] contains data ->\tkept\n";
62       ++$index;
63     }
64     else{
65 
66       my $delete = unlink $sorting_files[$index];
67 
68       if ($delete){
69 	warn "$sorting_files[$index] was empty ->\tdeleted\n";
70       }
71       else{
72 	warn "$sorting_files[$index] was empty, however deletion was unsuccessful: $!\n"
73       }
74 
75       ### we also need to remove the element from @sorting_files
76       splice @sorting_files, $index, 1;
77     }
78   }
79   warn "\n\n"; ## can't close the piped filehandles at this point because it will die (unfortunately)
80 }
81 
82 sub produce_mbias_plots{
83 
84   my $filename = shift;
85 
86   my $mbias = (split (/\//,$filename))[-1]; # extracting the filename if a full path was specified
87   $mbias =~ s/gz$//;
88   $mbias =~ s/sam$//;
89   $mbias =~ s/bam$//;
90   $mbias =~ s/txt$//;
91    my $mbias_graph_1 = my $mbias_graph_2 = $mbias;
92   $mbias_graph_1 = $output_dir . $mbias_graph_1 . 'M-bias_R1.png';
93   $mbias_graph_2 = $output_dir . $mbias_graph_2 . 'M-bias_R2.png';
94 
95   $mbias =~ s/$/M-bias.txt/;
96 
97   open (MBIAS,'>',"$output_dir$mbias") or die "Failed to open file for the M-bias data\n\n";
98 
99   # determining maximum read length
00   my $max_length_1 = 0;
01   my $max_length_2 = 0;
02 
03   foreach my $context (keys %mbias_1){
04     foreach my $pos (sort {$a<=>$b} keys %{$mbias_1{$context}}){
05       $max_length_1 = $pos unless ($max_length_1 >= $pos);
06     }
07   }
08   if ($paired){
09     foreach my $context (keys %mbias_2){
10       foreach my $pos (sort {$a<=>$b} keys %{$mbias_2{$context}}){
11 	$max_length_2 = $pos unless ($max_length_2 >= $pos);
12       }
13     }
14   }
15 
16   if ($single){
17     warn "Determining maximum read length for M-Bias plot\n";
18     warn "Maximum read length of Read 1: $max_length_1\n\n";
19   }
20   else{
21     warn "Determining maximum read lengths for M-Bias plots\n";
22     warn "Maximum read length of Read 1: $max_length_1\n";
23     warn "Maximum read length of Read 2: $max_length_2\n\n";
24   }
25   # sleep(3);
26 
27   my @mbias_read1;
28   my @mbias_read2;
29 
30   #Check whether the module GD::Graph:lines is installed
31   my $gd_graph_installed = 0;
32   eval{
33     require GD::Graph::lines;
34     GD::Graph::lines->import();
35   };
36 
37   unless($@) { # syntax or routine error variable, set if something goes wron in the last eval{ require ...}
38     $gd_graph_installed = 1;
39 
40     #Check whether the module GD::Graph::colour is installed
41     eval{
42       require GD::Graph::colour;
43       GD::Graph::colour->import(qw(:colours :lists :files :convert));
44     };
45 
46     if ($@) {
47       warn "Perl module GD::Graph::colour not found, skipping drawing M-bias plots (only writing out M-bias plot table)\n";
48       sleep(2);
49       $gd_graph_installed = 0;
50     }
51 
52 
53   }
54   else{
55     warn "Perl module GD::Graph::lines is not installed, skipping drawing M-bias plots (only writing out M-bias plot table)\n";
56     sleep(2);
57   }
58 
59 
60   my $graph_title;
61   my $graph1;
62   my $graph2;
63 
64   if ( $gd_graph_installed){
65     $graph1 = GD::Graph::lines->new(800,600);
66     if ($paired){
67       $graph2 = GD::Graph::lines->new(800,600);
68     }
69   }
70 
71   foreach my $context (qw(CpG CHG CHH)){
72     @{$mbias_read1[0]} = ();
73 
74     if ($paired){
75       print MBIAS "$context context (R1)\n================\n";
76       $graph_title = 'M-bias (Read 1)';
77     }
78     else{
79       print MBIAS "$context context\n===========\n";
80       $graph_title = 'M-bias';
81     }
82     print MBIAS "position\tcount methylated\tcount unmethylated\t% methylation\tcoverage\n";
83 
84     foreach my $pos (1..$max_length_1){
85 
86       unless (defined $mbias_1{$context}->{$pos}->{meth}){
87 	$mbias_1{$context}->{$pos}->{meth} = 0;
88       }
89       unless (defined $mbias_1{$context}->{$pos}->{un}){
90 	$mbias_1{$context}->{$pos}->{un} = 0;
91       }
92 
93       my $percent = '';
94       if (($mbias_1{$context}->{$pos}->{meth} + $mbias_1{$context}->{$pos}->{un}) > 0){
95 	$percent = sprintf("%.2f",$mbias_1{$context}->{$pos}->{meth} * 100/ ( $mbias_1{$context}->{$pos}->{meth} + $mbias_1{$context}->{$pos}->{un}) );
96       }
97       my $coverage = $mbias_1{$context}->{$pos}->{un} + $mbias_1{$context}->{$pos}->{meth};
98 
99       print MBIAS "$pos\t$mbias_1{$context}->{$pos}->{meth}\t$mbias_1{$context}->{$pos}->{un}\t$percent\t$coverage\n";
00       push @{$mbias_read1[0]},$pos;
01 
02       if ($context eq 'CpG'){
03 	push @{$mbias_read1[1]},$percent;
04 	push @{$mbias_read1[4]},$coverage;
05       }
06       elsif ($context eq 'CHG'){
07 	push @{$mbias_read1[2]},$percent;
08 	push @{$mbias_read1[5]},$coverage;
09       }
10       elsif ($context eq 'CHH'){
11     	push @{$mbias_read1[3]},$percent;
12 	push @{$mbias_read1[6]},$coverage;
13       }
14     }
15     print MBIAS "\n";
16   }
17 
18   if ( $gd_graph_installed){
19 
20     add_colour(nice_blue => [31,120,180]);
21     add_colour(nice_orange => [255,127,0]);
22     add_colour(nice_green => [51,160,44]);
23     add_colour(pale_blue => [153,206,227]);
24     add_colour(pale_orange => [253,204,138]);
25     add_colour(pale_green => [191,230,207]);
26 
27     $graph1->set(
28 		 x_label              => 'position (bp)',
29 		 y1_label              => '% methylation',
30 		 y2_label              => '# methylation calls',
31 		 title                => $graph_title,
32 		 line_width           => 2,
33 		 x_max_value          => $max_length_1,
34 		 x_min_value          => 0,
35 		 y_tick_number        => 10,
36 		 y_label_skip         => 2,
37 		 y1_max_value          => 100,
38 		 y1_min_value          => 0,
39 		 y_label_skip         => 2,
40 		 y2_min_value          => 0,
41 		 x_label_skip         => 5,
42 		 x_label_position     => 0.5,
43 		 x_tick_offset        => -1,
44 		 bgclr                => 'white',
45 		 transparent          => 0,
46 		 two_axes             => 1,
47 		 use_axis             => [1,1,1,2,2,2],
48 		 legend_placement     => 'RC',
49 		 legend_spacing       => 6,
50 		 legend_marker_width  => 24,
51 		 legend_marker_height => 18,
52 		 dclrs              => [ qw(nice_blue nice_orange nice_green pale_blue pale_orange pale_green)],
53 		) or die $graph1->error;
54 
55     $graph1->set_legend('CpG methylation','CHG methylation','CHH methylation','CpG total calls','CHG total calls','CHH total calls');
56 
57     my $gd1 = $graph1->plot(\@mbias_read1) or die $graph1->error;
58 
59     open (MBIAS_G1,'>',$mbias_graph_1) or die "Failed to write to file for M-bias plot 1: $!\n\n";
60     binmode MBIAS_G1;
61     print MBIAS_G1 $gd1->png;
62   }
63 
64   if ($paired){
65 
66     foreach my $context (qw(CpG CHG CHH)){
67       @{$mbias_read2[0]} = ();
68 
69       print MBIAS "$context context (R2)\n================\n";
70       print MBIAS "position\tcount methylated\tcount unmethylated\t% methylation\tcoverage\n";
71       foreach my $pos (1..$max_length_2){
72 	
73 	unless (defined $mbias_2{$context}->{$pos}->{meth}){
74 	  $mbias_2{$context}->{$pos}->{meth} = 0;
75 	}
76 	unless (defined $mbias_2{$context}->{$pos}->{un}){
77 	  $mbias_2{$context}->{$pos}->{un} = 0;
78 	}
79 
80 	my $percent = '';
81 	if (($mbias_2{$context}->{$pos}->{meth} + $mbias_2{$context}->{$pos}->{un}) > 0){
82 	  $percent = sprintf("%.2f",$mbias_2{$context}->{$pos}->{meth} * 100/ ($mbias_2{$context}->{$pos}->{meth} + $mbias_2{$context}->{$pos}->{un}) );
83 	}
84 	my $coverage = $mbias_2{$context}->{$pos}->{un} + $mbias_2{$context}->{$pos}->{meth};
85 
86 	print MBIAS "$pos\t$mbias_2{$context}->{$pos}->{meth}\t$mbias_2{$context}->{$pos}->{un}\t$percent\t$coverage\n";
87 	
88 	push @{$mbias_read2[0]},$pos;
89 	
90 	if ($context eq 'CpG'){
91 	  push @{$mbias_read2[1]},$percent;
92 	  push @{$mbias_read2[4]},$coverage;
93 	}
94 	elsif ($context eq 'CHG'){
95 	  push @{$mbias_read2[2]},$percent;
96 	  push @{$mbias_read2[5]},$coverage;
97 	}
98 	elsif ($context eq 'CHH'){
99 	  push @{$mbias_read2[3]},$percent;
00 	  push @{$mbias_read2[6]},$coverage;
01 	}
02       }	
03       print MBIAS "\n";
04     }
05 
06     if ( $gd_graph_installed){
07 
08       add_colour(nice_blue => [31,120,180]);
09       add_colour(nice_orange => [255,127,0]);
10       add_colour(nice_green => [51,160,44]);
11       add_colour(pale_blue => [153,206,227]);
12       add_colour(pale_orange => [253,204,138]);
13       add_colour(pale_green => [191,230,207]);
14 
15       $graph2->set(
16 		   x_label              => 'position (bp)',
17 		   line_width           => 2,
18 		   x_max_value          => $max_length_1,
19 		   x_min_value          => 0,
20 		   y_tick_number        => 10,
21 		   y_label_skip         => 2,
22 		   y1_max_value          => 100,
23 		   y1_min_value          => 0,
24 		   y_label_skip         => 2,
25 		   y2_min_value          => 0,
26 		   x_label_skip         => 5,
27 		   x_label_position     => 0.5,
28 		   x_tick_offset        => -1,
29 		   bgclr                => 'white',
30 		   transparent          => 0,
31 		   two_axes             => 1,
32 		   use_axis             => [1,1,1,2,2,2],
33 		   legend_placement     => 'RC',
34 		   legend_spacing       => 6,
35 		   legend_marker_width  => 24,
36 		   legend_marker_height => 18,
37 		   dclrs                => [ qw(nice_blue nice_orange nice_green pale_blue pale_orange pale_green)],
38 		   x_label              => 'position (bp)',
39 		   y1_label             => '% methylation',
40 		   y2_label             => '# calls',
41 		   title                => 'M-bias (Read 2)',
42 		  ) or die $graph2->error;
43 
44       $graph2->set_legend('CpG methylation','CHG methylation','CHH methylation','CpG total calls','CHG total calls','CHH total calls');
45       my $gd2 = $graph2->plot(\@mbias_read2) or die $graph2->error;
46 
47       open (MBIAS_G2,'>',$mbias_graph_2) or die "Failed to write to file for M-bias plot 2: $!\n\n";
48       binmode MBIAS_G2;
49       print MBIAS_G2 $gd2->png;
50 
51     }
52   }
53 }
54 
55 sub process_commandline{
56   my $help;
57   my $single_end;
58   my $paired_end;
59   my $ignore;
60   my $ignore_r2;
61   my $genomic_fasta;
62   my $full;
63   my $report;
64   my $extractor_version;
65   my $no_overlap;
66   my $merge_non_CpG;
67   my $vanilla;
68   my $output_dir;
69   my $no_header;
70   my $bedGraph;
71   my $coverage_threshold = 1; # Minimum number of reads covering before calling methylation status
72   my $remove;
73   my $counts;
74   my $cytosine_report;
75   my $genome_folder;
76   my $zero;
77   my $CpG_only;
78   my $CX_context;
79   my $split_by_chromosome;
80   my $sort_size;
81   my $samtools_path;
82   my $gzip;
83   my $mbias_only;
84   my $gazillion;
85   my $ample_mem;
86 
87   my $command_line = GetOptions ('help|man'             => \$help,
88 				 'p|paired-end'         => \$paired_end,
89 				 's|single-end'         => \$single_end,
90 				 'fasta'                => \$genomic_fasta,
91 				 'ignore=i'             => \$ignore,
92 				 'ignore_r2=i'          => \$ignore_r2,
93 				 'comprehensive'        => \$full,
94 				 'report'               => \$report,
95 				 'version'              => \$extractor_version,
96 				 'no_overlap'           => \$no_overlap,
97 				 'merge_non_CpG'        => \$merge_non_CpG,
98 				 'vanilla'              => \$vanilla,
99 				 'o|output=s'           => \$output_dir,
00 				 'no_header'            => \$no_header,
01 				 'bedGraph'             => \$bedGraph,
02 				 "cutoff=i"             => \$coverage_threshold,
03 				 "remove_spaces"        => \$remove,
04 				 "counts"               => \$counts,
05 				 "cytosine_report"      => \$cytosine_report,
06 				 'g|genome_folder=s'    => \$genome_folder,
07 				 "zero_based"           => \$zero,	
08 				 "CX|CX_context"        => \$CX_context,
09 				 "split_by_chromosome"  => \$split_by_chromosome,
10 				 "buffer_size=s"        => \$sort_size,
11 				 'samtools_path=s'      => \$samtools_path,
12 				 "gzip"                 => \$gzip,
13 				 "mbias_only"           => \$mbias_only,			
14 				 "gazillion|scaffolds"  => \$gazillion,
15 				 "ample_memory"         => \$ample_mem,
16 	);
17 
18   ### EXIT ON ERROR if there were errors with any of the supplied options
19   unless ($command_line){
20     die "Please respecify command line options\n";
21   }
22 
23   ### HELPFILE
24   if ($help){
25     print_helpfile();
26     exit;
27   }
28 
29   if ($extractor_version){
30     print << "VERSION";
31 
32 
33                            Bismark Methylation Extractor
34 
35                       Bismark Extractor Version: $version
36               Copyright 2010-13 Felix Krueger, Babraham Bioinformatics
37                 www.bioinformatics.babraham.ac.uk/projects/bismark/
38 
39 
40 VERSION
41     exit;
42   }
43 
44 
45   ### no files provided
46   unless (@ARGV){
47     die "You need to provide one or more Bismark files to create an individual C methylation output. Please respecify!\n";
48   }
49   @filenames = @ARGV;
50 
51   warn "\n *** Bismark methylation extractor version $version ***\n\n";
52 
53   ### M-BIAS ONLY
54   if ($mbias_only){
55     if ($bedGraph){
56       die "Option '--mbias_only' skips all sorts of methylation extraction, including the bedGraph generation. Please respecify!\n";
57     }
58     if ($cytosine_report){
59       die "Option '--mbias_only' skips all sorts of methylation extraction, including the genome-wide cytosine methylation report generation. Please respecify!\n";
60     }
61     if ($merge_non_CpG){
62       warn "Option '--mbias_only' skips all sorts of methylation extraction, thus '--merge' won't have any effect\n";
63     }
64     if ($full){
65       warn "Option '--mbias_only' skips all sorts of methylation extraction, thus '--comprehensive' won't have any effect\n";
66     }
67     sleep(3);
68   }
69 
70   ### PRINT A REPORT
71   unless ($report){
72     $report = 0;
73   }
74 
75   ### OUTPUT DIR PATH
76   if ($output_dir){
77     unless ($output_dir =~ /\/$/){
78       $output_dir =~ s/$/\//;
79     }
80   }
81   else{
82     $output_dir = '';
83   }
84 
85   ### NO HEADER
86   unless ($no_header){
87     $no_header = 0;
88   }
89 
90   ### OLD (VANILLA) OUTPUT FORMAT
91   unless ($vanilla){
92     $vanilla = 0;
93   }
94 
95   if ($single_end){
96     $paired_end = 0;   ### SINGLE END ALIGNMENTS
97   }
98   elsif ($paired_end){
99     $single_end = 0;   ### PAIRED-END ALIGNMENTS
00   }
01   else{
02 
03     ### we will try to determine whether the input file was a single-end or paired-end sequencing run from the SAM header
04 
05     if ($vanilla){
06       die "Please specify whether the supplied file(s) are in Bismark single-end or paired-end format with '-s' or '-p'\n\n";
07     }
08     else{ # SAM/BAM format
09 
10       my $file = $filenames[0];
11       warn "Trying to determine the type of mapping from the SAM header line of file $file\n"; sleep(1);
12 	
13       ### if the user did not specify whether the alignment file was single-end or paired-end we are trying to get this information from the @PG header line in the SAM/BAM file
14       if ($file =~ /\.gz$/){
15 	open (DETERMINE,"zcat $file |") or die "Unable to read from gzipped file $file: $!\n";
16       }
17       elsif ($file =~ /\.bam$/ ||  `file -b $file` =~ /^gzip/){
18 	open (DETERMINE,"samtools view -h $file |") or die "Unable to read from BAM file $file: $!\n";
19       }
20       else{
21 	open (DETERMINE,$file) or die "Unable to read from $file: $!\n";
22       }
23 
24       while (<DETERMINE>){
25 	last unless (/^\@/);
26 	if ($_ =~ /^\@PG/){
27 	  # warn "found the \@PG line:\n";
28 	  # warn "$_";
29 	
30 	  if ($_ =~ /-1/ and $_ =~ /-2/){
31 	    warn "Treating file(s) as paired-end data (as extracted from \@PG line)\n\n"; sleep(1);
32 	    $paired_end = 1;
33 	    $single_end = 0;
34 	  }
35 	  else{
36 	    warn "Treating file(s) as single-end data (as extracted from \@PG line)\n\n"; sleep(1);
37 	    $paired_end = 0;
38 	    $single_end = 1;
39 	  }
40 	}
41       }
42 
43       close DETERMINE or warn $!;
44 	
45     }
46   }
47 
48   ### IGNORING <INT> bases at the start of the read when processing the methylation call string
49   unless ($ignore){
50     $ignore = 0;
51   }
52 
53   if (defined $ignore_r2){
54     die "You can only specify --ignore_r2 for paired-end result files\n" unless ($paired_end);
55   }
56   else{
57     $ignore_r2 = 0;
58   }
59 
60 
61   ### NO OVERLAP
62   if ($no_overlap){
63     die "The option '--no_overlap' can only be specified for paired-end input!\n" unless ($paired_end);
64   }
65   else{
66     $no_overlap = 0;
67   }
68 
69   ### COMPREHENSIVE OUTPUT
70   unless ($full){
71     $full = 0;
72   }
73 
74   ### MERGE NON-CpG context
75   unless ($merge_non_CpG){
76     $merge_non_CpG = 0;
77   }
78 
79   ### remove white spaces in read ID (needed for sorting using the sort command
80   unless ($remove){
81     $remove = 0;
82   }
83 
84   ### COVERAGE THRESHOLD FOR bedGraph OUTPUT
85   if (defined $coverage_threshold){
86     unless ($coverage_threshold > 0){
87       die "Please select a coverage greater than 0 (positive integers only)\n";
88     }
89   }
90   else{
91     $coverage_threshold = 1;
92   }
93 
94   ### SORT buffer size
95   if (defined $sort_size){
96     unless ($sort_size =~ /^\d+\%$/ or $sort_size =~ /^\d+(K|M|G|T)$/){
97       die "Please select a buffer size as percentage (e.g. --buffer_size 20%) or a number to be multiplied with K, M, G, T etc. (e.g. --buffer_size 20G). For more information on sort type 'info sort' on a command line\n";
98     }
99   }
00   else{
01     $sort_size = '2G';
02   }
03 
04   if ($zero){
05     die "Option '--zero' is only available if  '--cytosine_report' is specified as well. Please respecify\n" unless ($cytosine_report);
06   }
07 
08   if ($CX_context){
09     die "Option '--CX_context' is only available if  '--cytosine_report' or '--bedGraph' is specified as well. Please respecify\n" unless ($cytosine_report or $bedGraph);
10   }
11   else{
12     $CX_context = 0;
13   }
14 
15   unless ($counts){
16     $counts = 1; # counts will always be set
17   }
18 
19   if ($cytosine_report){
20 
21     ### GENOME folder
22     if ($genome_folder){
23       unless ($genome_folder =~/\/$/){
24 	$genome_folder =~ s/$/\//;
25       }
26     }
27     else{
28       die "Please specify a genome folder to proceed (full path only)\n";
29     }
30 
31     unless ($bedGraph){
32       warn "Setting the option '--bedGraph' since this is required for the genome-wide cytosine report\n";
33       $bedGraph = 1;
34     }
35     unless ($counts){
36       # warn "Setting the option '--counts' since this is required for the genome-wide cytosine report\n";
37       $counts = 1;
38     }
39     warn "\n";
40   }
41 
42   ### PATH TO SAMTOOLS
43   if (defined $samtools_path){
44     # if Samtools was specified as full command
45     if ($samtools_path =~ /samtools$/){
46       if (-e $samtools_path){
47 	# Samtools executable found
48       }
49       else{
50 	die "Could not find an installation of Samtools at the location $samtools_path. Please respecify\n";
51       }
52     }
53     else{
54       unless ($samtools_path =~ /\/$/){
55 	$samtools_path =~ s/$/\//;
56       }
57       $samtools_path .= 'samtools';
58       if (-e $samtools_path){
59 	# Samtools executable found
60       }
61       else{
62 	die "Could not find an installation of Samtools at the location $samtools_path. Please respecify\n";
63       }
64     }
65   }
66   # Check whether Samtools is in the PATH if no path was supplied by the user
67   else{
68     if (!system "which samtools >/dev/null 2>&1"){ # STDOUT is binned, STDERR is redirected to STDOUT. Returns 0 if Samtools is in the PATH
69       $samtools_path = `which samtools`;
70       chomp $samtools_path;
71     }
72   }
73 
74   unless (defined $samtools_path){
75     $samtools_path = '';
76   }
77 
78 
79   if ($gazillion){
80     if ($ample_mem){
81       die "You can't currently select '--ample_mem' together with '--gazillion'. Make your pick!\n\n";
82     }
83   }
84 
85   return ($ignore,$genomic_fasta,$single_end,$paired_end,$full,$report,$no_overlap,$merge_non_CpG,$vanilla,$output_dir,$no_header,$bedGraph,$remove,$coverage_threshold,$counts,$cytosine_report,$genome_folder,$zero,$CpG_only,$CX_context,$split_by_chromosome,$sort_size,$samtools_path,$gzip,$ignore_r2,$mbias_only,$gazillion,$ample_mem);
86 }
87 
88 
89 sub test_positional_sorting{
90 
91   my $filename = shift;
92 
93   print "\nNow testing Bismark result file $filename for positional sorting (which would be bad...)\t";
94   sleep(1);
95 
96   if ($filename =~ /\.gz$/) {
97     open (TEST,"zcat $filename |") or die "Can't open gzipped file $filename: $!\n";
98   }
99   elsif ($filename =~ /bam$/ || `file -b $filename` =~ /^gzip/) {
00     if ($samtools_path){
01       open (TEST,"$samtools_path view -h $filename |") or die "Can't open BAM file $filename: $!\n";
02     }
03     else{
04       die "Sorry couldn't find an installation of Samtools. Either specifiy an alternative path using the option '--samtools_path /your/path/', or use a SAM file instead\n\n";
05     }
06   }
07   else {
08     open (TEST,$filename) or die "Can't open file $filename: $!\n";
09   }
10 
11   my $count = 0;
12 
13   while (<TEST>) {
14     if (/^\@/) {	     # testing header lines if they contain the @SO flag (for being sorted)
15       if (/^\@SO/) {
16 	die "SAM/BAM header line '$_' indicates that the Bismark aligment file has been sorted by chromosomal positions which is is incompatible with correct methylation extraction. Please use an unsorted file instead\n\n";
17       }
18       next;
19     }
20     $count++;
21 
22     last if ($count > 100000); # else we test the first 100000 sequences if they start with the same read ID
23 
24     my ($id_1) = (split (/\t/));
25 
26     ### reading the next line which should be read 2
27     $_ = <TEST>;
28     my ($id_2) = (split (/\t/));
29     last unless ($id_2);
30     ++$count;
31 
32     if ($id_1 eq $id_2){
33       ### ids are the same
34       next;
35     }
36     else{ ### in previous versions of Bismark we appended /1 and /2 to the read IDs for easier eyeballing which read is which. These tags need to be removed first
37       my $id_1_trunc = $id_1;
38       $id_1_trunc =~ s/\/1$//;
39       my $id_2_trunc = $id_2;
40       $id_2_trunc =~ s/\/2$//;
41 
42       unless ($id_1_trunc eq $id_2_trunc){
43 	die "The IDs of Read 1 ($id_1) and Read 2 ($id_2) are not the same. This might be a result of sorting the paired-end SAM/BAM files by chromosomal position which is not compatible with correct methylation extraction. Please use an unsorted file instead\n\n";
44       }
45     }
46   }
47   #  close TEST or die $!; somehow fails on our cluster...
48   ### If it hasen't died so far then it seems the file is in the correct Bismark format (read 1 and read 2 of a pair directly following each other)
49   warn "...passed!\n";
50   sleep(1);
51 
52 }
53 
54 
55 sub process_Bismark_results_file{
56   my $filename = shift;
57 
58   warn "\nNow reading in Bismark result file $filename\n\n";
59 
60   if ($filename =~ /\.gz$/) {
61     open (IN,"zcat $filename |") or die "Can't open gzipped file $filename: $!\n";
62   }
63   elsif ($filename =~ /bam$/ || `file -b $filename` =~ /^gzip/) {
64     if ($samtools_path){
65       open (IN,"$samtools_path view -h $filename |") or die "Can't open BAM file $filename: $!\n";
66     }
67     else{
68       die "Sorry couldn't find an installation of Samtools. Either specifiy an alternative path using the option '--samtools_path /your/path/', or use a SAM file instead\n\n";
69     }
70   }
71   else {
72     open (IN,$filename) or die "Can't open file $filename: $!\n";
73   }
74 
75   ### Vanilla and SAM output need to read different numbers of header lines
76   if ($vanilla) {
77     my $bismark_version = <IN>; ## discarding the Bismark version info
78     chomp $bismark_version;
79     $bismark_version =~ s/\r//; # replaces \r line feed
80     $bismark_version =~  s/Bismark version: //;
81     if ($bismark_version =~ /^\@/) {
82       warn "Detected \@ as the first character of the version information. Is it possible that the file is in SAM format?\n\n";
83       sleep (2);
84     }
85 
86     unless ($version eq $bismark_version){
87       die "The methylation extractor and Bismark itself need to be of the same version!\n\nVersions used:\nmethylation extractor: '$version'\nBismark: '$bismark_version'\n";
88     }
89   } else {
90     # If the read is in SAM format (default) it can either start with @ header lines or start with alignments directly.
91     # We are reading from it further down
92   }
93 
94   my $output_filename = (split (/\//,$filename))[-1];
95 
96   ### OPENING OUTPUT-FILEHANDLES
97   if ($report) {
98     my $report_filename = $output_filename;
99     $report_filename =~ s/\.sam$//;
00     $report_filename =~ s/\.txt$//;
01     $report_filename =~ s/$/_splitting_report.txt/;
02     $report_filename = $output_dir . $report_filename;
03     open (REPORT,'>',$report_filename) or die "Failed to write to file $report_filename $!\n";
04   }
05 
06   if ($report) {
07     print REPORT "$output_filename\n\n";
08     print REPORT "Parameters used to extract methylation information:\n";
09     if ($paired) {
10       if ($vanilla) {
11 	print REPORT "Bismark result file: paired-end (vanilla Bismark format)\n";
12       } else {
13 	print REPORT "Bismark result file: paired-end (SAM format)\n"; # default
14       }
15     }
16 
17     if ($single) {
18       if ($vanilla) {
19 	print REPORT "Bismark result file: single-end (vanilla Bismark format)\n";
20       } else {
21 	print REPORT "Bismark result file: single-end (SAM format)\n"; # default
22       }
23     }
24     if ($single){
25       if ($ignore) {
26 	print REPORT "Ignoring first $ignore bp\n";
27       }
28     }
29     else{ # paired-end
30       if ($ignore) {
31 	print REPORT "Ignoring first $ignore bp of Read 1\n";
32       }
33       if ($ignore_r2){
34 	print REPORT "Ignoring first $ignore_r2 bp of Read 2\n";
35       }
36     }
37 
38     if ($full) {
39       print REPORT "Output specified: comprehensive\n";
40     } else {
41       print REPORT "Output specified: strand-specific (default)\n";
42     }
43 
44     if ($no_overlap) {
45       print REPORT "No overlapping methylation calls specified\n";
46     }
47     if ($genomic_fasta) {
48       print REPORT "Genomic equivalent sequences will be printed out in FastA format\n";
49     }
50     if ($merge_non_CpG) {
51       print REPORT "Methylation in CHG and CHH context will be merged into \"non-CpG context\" output\n";
52     }
53 
54     print REPORT "\n";
55   }
56 
57 #####   open (OUT,"| gzip -c - > $output_dir$outfile") or die "Failed to write to $outfile: $!\n";
58 
59   ### CpG-context and non-CpG context. THIS SECTION IS OPTIONAL
60   ### if --comprehensive AND --merge_non_CpG was specified we are only writing out one CpG-context and one Any-Other-context result file
61   if ($full and $merge_non_CpG) {
62     my $cpg_output = my $other_c_output = $output_filename;
63     ### C in CpG context
64     $cpg_output =~ s/^/CpG_context_/;
65     $cpg_output =~ s/sam$/txt/;
66     $cpg_output =~ s/bam$/txt/;
67     $cpg_output =~ s/$/.txt/ unless ($cpg_output =~ /\.txt$/);
68     $cpg_output = $output_dir . $cpg_output;
69 
70     if ($gzip){
71       $cpg_output .= '.gz';
72       open ($fhs{CpG_context},"| gzip -c - > $cpg_output") or die "Failed to write to $cpg_output $! \n" unless($mbias_only);
73     }
74     else{
75       open ($fhs{CpG_context},'>',$cpg_output) or die "Failed to write to $cpg_output $! \n" unless($mbias_only);
76     }
77 
78     warn "Writing result file containing methylation information for C in CpG context to $cpg_output\n" unless($mbias_only);
79     push @sorting_files,$cpg_output;
80 
81     unless ($no_header) {
82       print {$fhs{CpG_context}} "Bismark methylation extractor version $version\n" unless($mbias_only);
83     }
84 
85     ### C in any other context than CpG
86     $other_c_output =~ s/^/Non_CpG_context_/;
87     $other_c_output =~ s/sam$/txt/;
88     $other_c_output =~ s/bam$/txt/;
89     $other_c_output =~ s/$/.txt/ unless ($other_c_output =~ /\.txt$/);
90     $other_c_output = $output_dir . $other_c_output;
91 
92     if ($gzip){
93       $other_c_output .= '.gz';
94       open ($fhs{other_context},"| gzip -c - > $other_c_output") or die "Failed to write to $other_c_output $! \n" unless($mbias_only);
95     }
96     else{
97       open ($fhs{other_context},'>',$other_c_output) or die "Failed to write to $other_c_output $!\n" unless($mbias_only);
98     }
99 
00     warn "Writing result file containing methylation information for C in any other context to $other_c_output\n" unless($mbias_only);
01     push @sorting_files,$other_c_output;
02 
03 
04     unless ($no_header) {
05       print {$fhs{other_context}} "Bismark methylation extractor version $version\n" unless($mbias_only);
06     }
07   }
08 
09   ### if only --merge_non_CpG was specified we will write out 8 different output files, depending on where the (first) unique best alignment has been found
10   elsif ($merge_non_CpG) {
11 
12     my $cpg_ot = my $cpg_ctot = my $cpg_ctob = my $cpg_ob = $output_filename;
13 
14     ### For cytosines in CpG context
15     $cpg_ot =~ s/^/CpG_OT_/;
16     $cpg_ot =~ s/sam$/txt/;
17     $cpg_ot =~ s/bam$/txt/;
18     $cpg_ot =~ s/$/.txt/ unless ($cpg_ot =~ /\.txt$/);
19     $cpg_ot = $output_dir . $cpg_ot;
20 
21     if ($gzip){
22       $cpg_ot .= '.gz';
23       open ($fhs{0}->{CpG},"| gzip -c - > $cpg_ot") or die "Failed to write to $cpg_ot $!\n" unless($mbias_only);
24     }
25     else{
26       open ($fhs{0}->{CpG},'>',$cpg_ot) or die "Failed to write to $cpg_ot $!\n" unless($mbias_only);
27     }
28 
29     warn "Writing result file containing methylation information for C in CpG context from the original top strand to $cpg_ot\n" unless($mbias_only);
30     push @sorting_files,$cpg_ot;
31 
32     unless($no_header){
33       print {$fhs{0}->{CpG}} "Bismark methylation extractor version $version\n" unless($mbias_only);
34     }
35 
36     $cpg_ctot =~ s/^/CpG_CTOT_/;
37     $cpg_ctot =~ s/sam$/txt/;
38     $cpg_ctot =~ s/bam$/txt/;
39     $cpg_ctot =~ s/$/.txt/ unless ($cpg_ctot =~ /\.txt$/);
40     $cpg_ctot = $output_dir . $cpg_ctot;
41 
42     if ($gzip){
43       $cpg_ctot .= '.gz';
44       open ($fhs{1}->{CpG},"| gzip -c - > $cpg_ctot") or die "Failed to write to $cpg_ctot $!\n" unless($mbias_only);
45     }
46     else{
47       open ($fhs{1}->{CpG},'>',$cpg_ctot) or die "Failed to write to $cpg_ctot $!\n" unless($mbias_only);
48     }
49 
50     warn "Writing result file containing methylation information for C in CpG context from the complementary to original top strand to $cpg_ctot\n" unless($mbias_only);
51     push @sorting_files,$cpg_ctot;
52 
53     unless($no_header){
54       print {$fhs{1}->{CpG}} "Bismark methylation extractor version $version\n" unless($mbias_only);
55     }
56 
57     $cpg_ctob =~ s/^/CpG_CTOB_/;
58     $cpg_ctob =~ s/sam$/txt/;
59     $cpg_ctob =~ s/bam$/txt/;
60     $cpg_ctob =~ s/$/.txt/ unless ($cpg_ctob =~ /\.txt$/);
61     $cpg_ctob = $output_dir . $cpg_ctob;
62 
63     if ($gzip){
64       $cpg_ctob .= '.gz';
65       open ($fhs{2}->{CpG},"| gzip -c - > $cpg_ctob") or die "Failed to write to $cpg_ctob $!\n" unless($mbias_only);
66     }
67     else{
68       open ($fhs{2}->{CpG},'>',$cpg_ctob) or die "Failed to write to $cpg_ctob $!\n" unless($mbias_only);
69     }
70 
71     warn "Writing result file containing methylation information for C in CpG context from the complementary to original bottom strand to $cpg_ctob\n" unless($mbias_only);
72     push @sorting_files,$cpg_ctob;
73 
74     unless($no_header){
75       print {$fhs{2}->{CpG}}  "Bismark methylation extractor version $version\n" unless($mbias_only);
76     }
77 
78     $cpg_ob =~ s/^/CpG_OB_/;
79     $cpg_ob =~ s/sam$/txt/;
80     $cpg_ob =~ s/bam$/txt/;
81     $cpg_ob =~ s/$/.txt/ unless ($cpg_ob =~ /\.txt$/);
82     $cpg_ob = $output_dir . $cpg_ob;
83 
84     if ($gzip){
85       $cpg_ob .= '.gz';
86       open ($fhs{3}->{CpG},"| gzip -c - > $cpg_ob") or die "Failed to write to $cpg_ob $!\n" unless($mbias_only);
87     }
88     else{
89       open ($fhs{3}->{CpG},'>',$cpg_ob) or die "Failed to write to $cpg_ob $!\n" unless($mbias_only);
90     }
91 
92     warn "Writing result file containing methylation information for C in CpG context from the original bottom strand to $cpg_ob\n\n" unless($mbias_only);
93     push @sorting_files,$cpg_ob;
94 
95     unless($no_header){
96       print {$fhs{3}->{CpG}}  "Bismark methylation extractor version $version\n" unless($mbias_only);
97     }
98 
99     ### For cytosines in Non-CpG (CC, CT or CA) context
00     my $other_c_ot = my $other_c_ctot = my $other_c_ctob = my $other_c_ob = $output_filename;
01 
02     $other_c_ot =~ s/^/Non_CpG_OT_/;
03     $other_c_ot =~ s/sam$/txt/;
04     $other_c_ot =~ s/bam$/txt/;
05     $other_c_ot =~ s/$/.txt/ unless ($other_c_ot =~ /\.txt$/);
06     $other_c_ot = $output_dir . $other_c_ot;
07 
08     if ($gzip){
09       $other_c_ot .= '.gz';
10       open ($fhs{0}->{other_c},"| gzip -c - > $other_c_ot") or die "Failed to write to $other_c_ot $!\n" unless($mbias_only);
11     }
12     else{
13       open ($fhs{0}->{other_c},'>',$other_c_ot) or die "Failed to write to $other_c_ot $!\n" unless($mbias_only);
14     }
15 
16     warn "Writing result file containing methylation information for C in any other context from the original top strand to $other_c_ot\n" unless($mbias_only);
17     push @sorting_files,$other_c_ot;
18 
19     unless($no_header){
20       print {$fhs{0}->{other_c}} "Bismark methylation extractor version $version\n" unless($mbias_only);
21     }
22 
23     $other_c_ctot =~ s/^/Non_CpG_CTOT_/;
24     $other_c_ctot =~ s/sam$/txt/;
25     $other_c_ctot =~ s/bam$/txt/;
26     $other_c_ctot =~ s/$/.txt/ unless ($other_c_ctot =~ /\.txt$/);
27     $other_c_ctot = $output_dir . $other_c_ctot;
28 
29     if ($gzip){
30       $other_c_ctot .= '.gz';
31       open ($fhs{1}->{other_c},"| gzip -c - > $other_c_ctot") or die "Failed to write to $other_c_ctot $!\n" unless($mbias_only);
32     }
33     else{
34       open ($fhs{1}->{other_c},'>',$other_c_ctot) or die "Failed to write to $other_c_ctot $!\n" unless($mbias_only);
35     }
36 
37     warn "Writing result file containing methylation information for C in any other context from the complementary to original top strand to $other_c_ctot\n" unless($mbias_only);
38     push @sorting_files,$other_c_ctot;
39 
40     unless($no_header){
41       print {$fhs{1}->{other_c}} "Bismark methylation extractor version $version\n" unless($mbias_only);
42     }
43 
44     $other_c_ctob =~ s/^/Non_CpG_CTOB_/;
45     $other_c_ctob =~ s/sam$/txt/;
46     $other_c_ctob =~ s/bam$/txt/;
47     $other_c_ctob =~ s/$/.txt/ unless ($other_c_ctob =~ /\.txt$/);
48     $other_c_ctob = $output_dir . $other_c_ctob;
49 
50     if ($gzip){
51       $other_c_ctob .= '.gz';
52       open ($fhs{2}->{other_c},"| gzip -c - > $other_c_ctob") or die "Failed to write to $other_c_ctob $!\n" unless($mbias_only);
53     }
54     else{
55       open ($fhs{2}->{other_c},'>',$other_c_ctob) or die "Failed to write to $other_c_ctob $!\n" unless($mbias_only);
56     }
57 
58     warn "Writing result file containing methylation information for C in any other context from the complementary to original bottom strand to $other_c_ctob\n" unless($mbias_only);
59     push @sorting_files,$other_c_ctob;
60 
61     unless($no_header){
62       print {$fhs{2}->{other_c}} "Bismark methylation extractor version $version\n" unless($mbias_only);
63     }
64 
65     $other_c_ob =~ s/^/Non_CpG_OB_/;
66     $other_c_ob =~ s/sam$/txt/;
67     $other_c_ob =~ s/sam$/txt/;
68     $other_c_ob =~ s/$/.txt/ unless ($other_c_ob =~ /\.txt$/);
69     $other_c_ob = $output_dir . $other_c_ob;
70 
71     if ($gzip){
72       $other_c_ob .= '.gz';
73       open ($fhs{3}->{other_c},"| gzip -c - > $other_c_ob") or die "Failed to write to $other_c_ob $!\n" unless($mbias_only);
74     }
75     else{
76       open ($fhs{3}->{other_c},'>',$other_c_ob) or die "Failed to write to $other_c_ob $!\n" unless($mbias_only);
77     }
78 
79     warn "Writing result file containing methylation information for C in any other context from the original bottom strand to $other_c_ob\n\n" unless($mbias_only);
80     push @sorting_files,$other_c_ob;
81 
82     unless($no_header){
83       print {$fhs{3}->{other_c}} "Bismark methylation extractor version $version\n" unless($mbias_only);
84     }
85   }
86   ### THIS SECTION IS THE DEFAULT (CpG, CHG and CHH context)
87 
88   ### if --comprehensive was specified we are only writing one file per context
89   elsif ($full) {
90     my $cpg_output = my $chg_output =  my $chh_output = $output_filename;
91     ### C in CpG context
92     $cpg_output =~ s/^/CpG_context_/;
93     $cpg_output =~ s/sam$/txt/;
94     $cpg_output =~ s/bam$/txt/;
95     $cpg_output =~ s/$/.txt/ unless ($cpg_output =~ /\.txt$/);
96     $cpg_output = $output_dir . $cpg_output;
97 
98     if ($gzip){
99       $cpg_output .= '.gz';
00       open ($fhs{CpG_context},"| gzip -c - > $cpg_output") or die "Failed to write to $cpg_output $! \n" unless($mbias_only);
01     }
02     else{
03       open ($fhs{CpG_context},'>',$cpg_output) or die "Failed to write to $cpg_output $! \n" unless($mbias_only);
04     }
05 
06     warn "Writing result file containing methylation information for C in CpG context to $cpg_output\n" unless($mbias_only);
07     push @sorting_files,$cpg_output;
08 
09     unless($no_header){
10       print {$fhs{CpG_context}} "Bismark methylation extractor version $version\n" unless($mbias_only);
11     }
12 
13     ### C in CHG context
14     $chg_output =~ s/^/CHG_context_/;
15     $chg_output =~ s/sam$/txt/;
16     $chg_output =~ s/bam$/txt/;
17     $chg_output =~ s/$/.txt/ unless ($chg_output =~ /\.txt$/);
18     $chg_output = $output_dir . $chg_output;
19 
20     if ($gzip){
21       $chg_output .= '.gz';
22       open ($fhs{CHG_context},"| gzip -c - > $chg_output") or die "Failed to write to $chg_output $!\n" unless($mbias_only);
23     }
24     else{
25       open ($fhs{CHG_context},'>',$chg_output) or die "Failed to write to $chg_output $!\n" unless($mbias_only);
26     }
27 
28     warn "Writing result file containing methylation information for C in CHG context to $chg_output\n" unless($mbias_only);
29     push @sorting_files,$chg_output;
30 
31     unless($no_header){
32       print {$fhs{CHG_context}} "Bismark methylation extractor version $version\n" unless($mbias_only);
33     }
34 
35     ### C in CHH context
36     $chh_output =~ s/^/CHH_context_/;
37     $chh_output =~ s/sam$/txt/;
38     $chh_output =~ s/bam$/txt/;
39     $chh_output =~ s/$/.txt/ unless ($chh_output =~ /\.txt$/);
40     $chh_output = $output_dir . $chh_output;
41 
42     if ($gzip){
43       $chh_output .= '.gz';
44       open ($fhs{CHH_context},"| gzip -c - > $chh_output") or die "Failed to write to $chh_output $!\n" unless($mbias_only);
45     }
46     else{
47       open ($fhs{CHH_context},'>',$chh_output) or die "Failed to write to $chh_output $!\n" unless($mbias_only);
48     }
49 
50     warn "Writing result file containing methylation information for C in CHH context to $chh_output\n" unless($mbias_only);
51     push @sorting_files, $chh_output;
52 
53     unless($no_header){
54       print {$fhs{CHH_context}} "Bismark methylation extractor version $version\n" unless($mbias_only);
55     }
56   }
57   ### else we will write out 12 different output files, depending on where the (first) unique best alignment was found
58   else {
59     my $cpg_ot = my $cpg_ctot = my $cpg_ctob = my $cpg_ob = $output_filename;
60 
61     ### For cytosines in CpG context
62     $cpg_ot =~ s/^/CpG_OT_/;
63     $cpg_ot =~ s/sam$/txt/;
64     $cpg_ot =~ s/bam$/txt/;
65     $cpg_ot =~ s/$/.txt/ unless ($cpg_ot =~ /\.txt$/);
66     $cpg_ot = $output_dir . $cpg_ot;
67 
68     if ($gzip){
69       $cpg_ot .= '.gz';
70       open ($fhs{0}->{CpG},"| gzip -c - > $cpg_ot") or die "Failed to write to $cpg_ot $!\n" unless($mbias_only);
71     }
72     else{
73       open ($fhs{0}->{CpG},'>',$cpg_ot) or die "Failed to write to $cpg_ot $!\n" unless($mbias_only);
74     }
75 
76     warn "Writing result file containing methylation information for C in CpG context from the original top strand to $cpg_ot\n" unless($mbias_only);
77     push @sorting_files,$cpg_ot;
78 
79     unless($no_header){
80       print {$fhs{0}->{CpG}} "Bismark methylation extractor version $version\n" unless($mbias_only);
81     }
82 
83     $cpg_ctot =~ s/^/CpG_CTOT_/;
84     $cpg_ctot =~ s/sam$/txt/;
85     $cpg_ctot =~ s/bam$/txt/;
86     $cpg_ctot =~ s/$/.txt/ unless ($cpg_ctot =~ /\.txt$/);
87     $cpg_ctot = $output_dir . $cpg_ctot;
88 
89     if ($gzip){
90       $cpg_ctot .= '.gz';
91       open ($fhs{1}->{CpG},"| gzip -c - > $cpg_ctot") or die "Failed to write to $cpg_ctot $!\n" unless($mbias_only);
92     }
93     else{
94       open ($fhs{1}->{CpG},'>',$cpg_ctot) or die "Failed to write to $cpg_ctot $!\n" unless($mbias_only);
95     }
96 
97     warn "Writing result file containing methylation information for C in CpG context from the complementary to original top strand to $cpg_ctot\n" unless($mbias_only);
98     push @sorting_files,$cpg_ctot;
99 
00     unless($no_header){
01       print {$fhs{1}->{CpG}} "Bismark methylation extractor version $version\n" unless($mbias_only);
02     }
03 
04     $cpg_ctob =~ s/^/CpG_CTOB_/;
05     $cpg_ctob =~ s/sam$/txt/;
06     $cpg_ctob =~ s/bam$/txt/;
07     $cpg_ctob =~ s/$/.txt/ unless ($cpg_ctob =~ /\.txt$/);
08     $cpg_ctob = $output_dir . $cpg_ctob;
09 
10     if ($gzip){
11       $cpg_ctob .= '.gz';
12       open ($fhs{2}->{CpG},"| gzip -c - > $cpg_ctob") or die "Failed to write to $cpg_ctob $!\n" unless($mbias_only);
13     }
14     else{
15       open ($fhs{2}->{CpG},'>',$cpg_ctob) or die "Failed to write to $cpg_ctob $!\n" unless($mbias_only);
16     }
17 
18     warn "Writing result file containing methylation information for C in CpG context from the complementary to original bottom strand to $cpg_ctob\n" unless($mbias_only);
19     push @sorting_files,$cpg_ctob;
20 
21     unless($no_header){
22       print {$fhs{2}->{CpG}}  "Bismark methylation extractor version $version\n" unless($mbias_only);
23     }
24 
25     $cpg_ob =~ s/^/CpG_OB_/;
26     $cpg_ob =~ s/sam$/txt/;
27     $cpg_ob =~ s/bam$/txt/;
28     $cpg_ob =~ s/$/.txt/ unless ($cpg_ob =~ /\.txt$/);
29     $cpg_ob = $output_dir . $cpg_ob;
30 
31     if ($gzip){
32       $cpg_ob .= '.gz';
33       open ($fhs{3}->{CpG},"| gzip -c - > $cpg_ob") or die "Failed to write to $cpg_ob $!\n" unless($mbias_only);
34     }
35     else{
36       open ($fhs{3}->{CpG},'>',$cpg_ob) or die "Failed to write to $cpg_ob $!\n" unless($mbias_only);
37     }
38 
39     warn "Writing result file containing methylation information for C in CpG context from the original bottom strand to $cpg_ob\n\n" unless($mbias_only);
40     push @sorting_files,$cpg_ob;
41 
42     unless($no_header){
43       print {$fhs{3}->{CpG}}  "Bismark methylation extractor version $version\n" unless($mbias_only);
44     }
45 
46     ### For cytosines in CHG context
47     my $chg_ot = my $chg_ctot = my $chg_ctob = my $chg_ob = $output_filename;
48 
49     $chg_ot =~ s/^/CHG_OT_/;
50     $chg_ot =~ s/sam$/txt/;
51     $chg_ot =~ s/bam$/txt/;
52     $chg_ot =~ s/$/.txt/ unless ($chg_ot =~ /\.txt$/);
53     $chg_ot = $output_dir . $chg_ot;
54 
55     if ($gzip){
56       $chg_ot .= '.gz';
57       open ($fhs{0}->{CHG},"| gzip -c - > $chg_ot") or die "Failed to write to $chg_ot $!\n" unless($mbias_only);
58     }
59     else{
60       open ($fhs{0}->{CHG},'>',$chg_ot) or die "Failed to write to $chg_ot $!\n" unless($mbias_only);
61     }
62 
63     warn "Writing result file containing methylation information for C in CHG context from the original top strand to $chg_ot\n" unless($mbias_only);
64     push @sorting_files,$chg_ot;
65 
66     unless($no_header){
67       print {$fhs{0}->{CHG}} "Bismark methylation extractor version $version\n" unless($mbias_only);
68     }
69 
70     $chg_ctot =~ s/^/CHG_CTOT_/;
71     $chg_ctot =~ s/sam$/txt/;
72     $chg_ctot =~ s/bam$/txt/;
73     $chg_ctot =~ s/$/.txt/ unless ($chg_ctot =~ /\.txt$/);
74     $chg_ctot = $output_dir . $chg_ctot;
75 
76     if ($gzip){
77       $chg_ctot .= '.gz';
78       open ($fhs{1}->{CHG},"| gzip -c - > $chg_ctot") or die "Failed to write to $chg_ctot $!\n" unless($mbias_only);
79     }
80     else{
81       open ($fhs{1}->{CHG},'>',$chg_ctot) or die "Failed to write to $chg_ctot $!\n" unless($mbias_only);
82     }
83 
84     warn "Writing result file containing methylation information for C in CHG context from the complementary to original top strand to $chg_ctot\n" unless($mbias_only);
85     push @sorting_files,$chg_ctot;
86 
87     unless($no_header){
88       print {$fhs{1}->{CHG}} "Bismark methylation extractor version $version\n" unless($mbias_only);
89     }
90 
91     $chg_ctob =~ s/^/CHG_CTOB_/;
92     $chg_ctob =~ s/sam$/txt/;
93     $chg_ctob =~ s/bam$/txt/;
94     $chg_ctob =~ s/$/.txt/ unless ($chg_ctob =~ /\.txt$/);
95     $chg_ctob = $output_dir . $chg_ctob;
96 
97     if ($gzip){
98       $chg_ctob .= '.gz';
99       open ($fhs{2}->{CHG},"| gzip -c - > $chg_ctob") or die "Failed to write to $chg_ctob $!\n" unless($mbias_only);
00     }
01     else{
02       open ($fhs{2}->{CHG},'>',$chg_ctob) or die "Failed to write to $chg_ctob $!\n" unless($mbias_only);
03     }
04 
05     warn "Writing result file containing methylation information for C in CHG context from the complementary to original bottom strand to $chg_ctob\n" unless($mbias_only);
06     push @sorting_files,$chg_ctob;
07 
08     unless($no_header){
09       print {$fhs{2}->{CHG}} "Bismark methylation extractor version $version\n" unless($mbias_only);
10     }
11 
12     $chg_ob =~ s/^/CHG_OB_/;
13     $chg_ob =~ s/sam$/txt/;
14     $chg_ob =~ s/bam$/txt/;
15     $chg_ob =~ s/$/.txt/ unless ($chg_ob =~ /\.txt$/);
16     $chg_ob = $output_dir . $chg_ob;
17 
18     if ($gzip){
19       $chg_ob .= '.gz';
20       open ($fhs{3}->{CHG},"| gzip -c - > $chg_ob") or die "Failed to write to $chg_ob $!\n" unless($mbias_only);
21     }
22     else{
23       open ($fhs{3}->{CHG},'>',$chg_ob) or die "Failed to write to $chg_ob $!\n" unless($mbias_only);
24     }
25 
26     warn "Writing result file containing methylation information for C in CHG context from the original bottom strand to $chg_ob\n\n" unless($mbias_only);
27     push @sorting_files,$chg_ob;
28 
29     unless($no_header){
30       print {$fhs{3}->{CHG}} "Bismark methylation extractor version $version\n" unless($mbias_only);
31     }
32 
33     ### For cytosines in CHH context
34     my $chh_ot = my $chh_ctot = my $chh_ctob = my $chh_ob = $output_filename;
35 
36     $chh_ot =~ s/^/CHH_OT_/;
37     $chh_ot =~ s/sam$/txt/;
38     $chh_ot =~ s/bam$/txt/;
39     $chh_ot =~ s/$/.txt/ unless ($chh_ot =~ /\.txt$/);
40     $chh_ot = $output_dir . $chh_ot;
41 
42     if ($gzip){
43       $chh_ot .= '.gz';
44       open ($fhs{0}->{CHH},"| gzip -c - > $chh_ot") or die "Failed to write to $chh_ot $!\n" unless($mbias_only);
45     }
46     else{
47       open ($fhs{0}->{CHH},'>',$chh_ot) or die "Failed to write to $chh_ot $!\n" unless($mbias_only);
48     }
49 
50     warn "Writing result file containing methylation information for C in CHH context from the original top strand to $chh_ot\n" unless($mbias_only);
51     push @sorting_files,$chh_ot;
52 
53     unless($no_header){
54       print {$fhs{0}->{CHH}} "Bismark methylation extractor version $version\n" unless($mbias_only);
55     }
56 
57     $chh_ctot =~ s/^/CHH_CTOT_/;
58     $chh_ctot =~ s/sam$/txt/;
59     $chh_ctot =~ s/bam$/txt/;
60     $chh_ctot =~ s/$/.txt/ unless ($chh_ctot =~ /\.txt$/);
61     $chh_ctot = $output_dir . $chh_ctot;
62 
63     if ($gzip){
64       $chh_ctot .= '.gz';
65       open ($fhs{1}->{CHH},"| gzip -c - > $chh_ctot") or die "Failed to write to $chh_ctot $!\n" unless($mbias_only);
66     }
67     else{
68       open ($fhs{1}->{CHH},'>',$chh_ctot) or die "Failed to write to $chh_ctot $!\n" unless($mbias_only);
69     }
70 
71     warn "Writing result file containing methylation information for C in CHH context from the complementary to original top strand to $chh_ctot\n" unless($mbias_only);
72     push @sorting_files,$chh_ctot;
73 
74     unless($no_header){
75       print {$fhs{1}->{CHH}} "Bismark methylation extractor version $version\n" unless($mbias_only);
76     }
77 
78     $chh_ctob =~ s/^/CHH_CTOB_/;
79     $chh_ctob =~ s/sam$/txt/;
80     $chh_ctob =~ s/bam$/txt/;
81     $chh_ctob =~ s/$/.txt/ unless ($chh_ctob =~ /\.txt$/);
82     $chh_ctob = $output_dir . $chh_ctob;
83 
84     if ($gzip){
85       $chh_ctob .= '.gz';
86       open ($fhs{2}->{CHH},"| gzip -c - > $chh_ctob") or die "Failed to write to $chh_ctob $!\n" unless($mbias_only);
87     }
88     else{
89       open ($fhs{2}->{CHH},'>',$chh_ctob) or die "Failed to write to $chh_ctob $!\n" unless($mbias_only);
90     }
91 
92     warn "Writing result file containing methylation information for C in CHH context from the complementary to original bottom strand to $chh_ctob\n" unless($mbias_only);
93     push @sorting_files,$chh_ctob;
94 
95     unless($no_header){
96       print {$fhs{2}->{CHH}} "Bismark methylation extractor version $version\n" unless($mbias_only);
97     }
98 
99     $chh_ob =~ s/^/CHH_OB_/;
00     $chh_ob =~ s/sam$/txt/;
01     $chh_ob =~ s/bam$/txt/;
02     $chh_ob =~ s/$/.txt/ unless ($chh_ob =~ /\.txt$/);
03     $chh_ob = $output_dir . $chh_ob;
04 
05     if ($gzip){
06       $chh_ob .= '.gz';
07       open ($fhs{3}->{CHH},"| gzip -c - > $chh_ob") or die "Failed to write to $chh_ob $!\n" unless($mbias_only);
08     }
09     else{
10       open ($fhs{3}->{CHH},'>',$chh_ob) or die "Failed to write to $chh_ob $!\n" unless($mbias_only);
11     }
12 
13     warn "Writing result file containing methylation information for C in CHH context from the original bottom strand to $chh_ob\n\n" unless($mbias_only);
14     push @sorting_files,$chh_ob;
15 
16     unless($no_header){
17       print {$fhs{3}->{CHH}} "Bismark methylation extractor version $version\n" unless($mbias_only);
18     }
19   }
20 
21   my $methylation_call_strings_processed = 0;
22   my $line_count = 0;
23 
24   ### proceeding differently now for single-end or paired-end Bismark files
25 
26   ### PROCESSING SINGLE-END RESULT FILES
27   if ($single) {
28 
29     ### also proceeding differently now for SAM format or vanilla Bismark format files
30     if ($vanilla) {		# old vanilla Bismark output format
31       while (<IN>) {
32 	++$line_count;
33 	warn "Processed lines: $line_count\n" if ($line_count%500000==0);
34 	
35 	### $seq here is the chromosomal sequence (to use for the repeat analysis for example)
36 	my ($id,$strand,$chrom,$start,$seq,$meth_call,$read_conversion,$genome_conversion) = (split("\t"))[0,1,2,3,6,7,8,9];
37 	
38 	### we need to remove 2 bp of the genomic sequence as we were extracting read + 2bp long fragments to make a methylation call at the first or
39 	### last position
40 	chomp $genome_conversion;
41 
42 	my $index;
43 	if ($meth_call) {
44 
45 	  if ($read_conversion eq 'CT' and $genome_conversion eq 'CT') { ## original top strand
46 	    $index = 0;
47 	  } elsif ($read_conversion eq 'GA' and $genome_conversion eq 'CT') { ## complementary to original top strand
48 	    $index = 1;
49 	  } elsif ($read_conversion eq 'CT' and $genome_conversion eq 'GA') { ## original bottom strand
50 	    $index = 3;
51 	  } elsif ($read_conversion eq 'GA' and $genome_conversion eq 'GA') { ## complementary to original bottom strand
52 	    $index = 2;
53 	  } else {
54 	    die "Unexpected combination of read and genome conversion: '$read_conversion' / '$genome_conversion'\n";
55 	  }
56 	
57 	  ### Clipping off the first <int> number of bases from the methylation call string as specified with --ignore <int>
58 	  if ($ignore) {
59 	    $meth_call = substr($meth_call,$ignore,length($meth_call)-$ignore);	
60 	
61 	    ### If we are clipping off some bases at the start we need to adjust the start position of the alignments accordingly!
62 	    if ($strand eq '+') {
63 	      $start += $ignore;
64 	    } elsif ($strand eq '-') {
65 	      $start += length($meth_call)-1; ## $meth_call is already shortened!
66 	    } else {
67 	      die "Alignment did not have proper strand information: $strand\n";
68 	    }
69 	  }
70 	  ### printing out the methylation state of every C in the read
71 	  print_individual_C_methylation_states_single_end($meth_call,$chrom,$start,$id,$strand,$index);
72 	
73 	  ++$methylation_call_strings_processed; # 1 per single-end result
74 	}
75       }
76     } else {		  # processing single-end SAM format (default)
77       while (<IN>) {
78 	### SAM format can either start with header lines (starting with @) or start with alignments directly
79 	if (/^\@/) {	     # skipping header lines (starting with @)
80 	  warn "skipping SAM header line:\t$_";
81 	  next;
82 	}
83 
84 	++$line_count;
85 	warn "Processed lines: $line_count\n" if ($line_count%500000==0);
86 	
87 	# example read in SAM format
88 	# 1_R1/1	67	5	103172224	255	40M	=	103172417	233	AATATTTTTTTTATTTTAAAATGTGTATTGATTTAAATTT	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:4	XX:Z:4T1T24TT7	XM:Z:....h.h........................hh.......	XR:Z:CT	XG:Z:CT
89 	###
90 
91 	# < 0.7.6 my ($id,$chrom,$start,$meth_call,$read_conversion,$genome_conversion) = (split("\t"))[0,2,3,13,14,15];
92 	# < 0.7.6 $meth_call =~ s/^XM:Z://;
93 	# < 0.7.6 $read_conversion =~ s/^XR:Z://;
94 	# < 0.7.6 $genome_conversion =~ s/^XG:Z://;	
95 
96 	my ($id,$chrom,$start,$cigar) = (split("\t"))[0,2,3,5];
97 
98 	### detecting the following SAM flags in case the SAM entry was shuffled by CRAM or Goby compression/decompression
99 	my $meth_call;	  ### Thanks to Zachary Zeno for this solution
00 	my $read_conversion;
01 	my $genome_conversion;
02 	
03 	while ( /(XM|XR|XG):Z:([^\t]+)/g ) {
04 	  my $tag = $1;
05 	  my $value = $2;
06 
07 	  if ($tag eq "XM") {
08 	    $meth_call = $value;
09 	    $meth_call =~ s/\r//;
10 	  } elsif ($tag eq "XR") {
11 	    $read_conversion = $value;
12 	    $read_conversion =~ s/\r//;
13 	  } elsif ($tag eq "XG") {
14 	    $genome_conversion = $value;
15 	    $genome_conversion =~ s/\r//;
16 	  }
17 	}
18 
19 	my $strand;
20 	chomp $genome_conversion;
21 	# print "$meth_call\n$read_conversion\n$genome_conversion\n";
22 	
23 	my $index;
24 	if ($meth_call) {
25 	  if ($read_conversion eq 'CT' and $genome_conversion eq 'CT') { ## original top strand
26 	    $index = 0;
27 	    $strand = '+';
28 	  } elsif ($read_conversion eq 'GA' and $genome_conversion eq 'CT') { ## complementary to original top strand
29 	    $index = 1;
30 	    $strand = '-';
31 	  } elsif ($read_conversion eq 'GA' and $genome_conversion eq 'GA') { ## complementary to original bottom strand
32 	    $index = 2;
33 	    $strand = '+';
34 	  } elsif ($read_conversion eq 'CT' and $genome_conversion eq 'GA') { ## original bottom strand
35 	    $index = 3;
36 	    $strand = '-';
37 	  } else {
38 	    die "Unexpected combination of read and genome conversion: '$read_conversion' / '$genome_conversion'\n";
39 	  }
40 	
41 	  ### If the read is in SAM format we need to reverse the methylation call if the read has been reverse-complemented for the output
42 	  if ($strand eq '-') {
43 	    $meth_call = reverse $meth_call;
44 	  }
45 	
46 	  ### Clipping off the first <int> number of bases from the methylation call string as specified with --ignore <int>
47 	  if ($ignore) {
48 	    # print "\n\n$meth_call\n";
49 	    $meth_call = substr($meth_call,$ignore,length($meth_call)-$ignore);	
50 	    # print "$meth_call\n";
51 
52 	    ### If we are ignoring a part of the sequence we also need to adjust the cigar string accordingly
53 
54 	    my @len = split (/\D+/,$cigar); # storing the length per operation
55 	    my @ops = split (/\d+/,$cigar); # storing the operation
56 	    shift @ops;		# remove the empty first element
57 	    die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len == scalar @ops);
58 		
59 	    my @comp_cigar; # building an array with all CIGAR operations
60 	    foreach my $index (0..$#len) {
61 	      foreach (1..$len[$index]) {
62 		# print  "$ops[$index]";
63 		push @comp_cigar, $ops[$index];
64 	      }
65 	    }
66 	    # print "original CIGAR: $cigar\n";
67 	    # print "original CIGAR: @comp_cigar\n";
68 
69 	    ### If we are clipping off some bases at the start we need to adjust the start position of the alignments accordingly!
70 	    if ($strand eq '+') {
71 	
72 	      my $D_count = 0; # counting all deletions that affect the ignored genomic position, i.e. Deletions and insertions
73 	      my $I_count = 0;
74 
75 	      for (1..$ignore) {
76 		my $op = shift @comp_cigar; # adjusting composite CIGAR string by removing $ignore operations from the start
77 		# print "$_ deleted $op\n";
78 
79 		while ($op eq 'D') { # repeating this for deletions (D)
80 		  $D_count++;
81 		  $op = shift @comp_cigar;
82 		  # print "$_ deleted $op\n";
83 		}
84 		if ($op eq 'I') { # adjusting the genomic position for insertions (I)
85 		  $I_count++;
86 		}
87 	      }
88 	      $start += $ignore + $D_count - $I_count;
89 	      # print "start $start\t ignore: $ignore\t D count: $D_count I_count: $I_count\n";
90 	    } elsif ($strand eq '-') {
91 
92 	      for (1..$ignore) {
93 		my $op = pop @comp_cigar; # adjusting composite CIGAR string by removing $ignore operations, here the last value of the array
94 		while ($op eq 'D') { # repeating this for deletions (D)
95 		  $op = pop @comp_cigar;
96 		}
97 	      }
98 
99 	      ### For reverse strand alignments we need to determine the number of matching bases (M) or deletions (D) in the read from the CIGAR
00 	      ### string to be able to work out the starting position of the read which is on the 3' end of the sequence
01 	      my $MD_count = 0;	# counting all operations that affect the genomic position, i.e. M and D. Insertions do not affect the start position
02 	      foreach (@comp_cigar) {
03 		++$MD_count if ($_ eq 'M' or $_ eq 'D');
04 	      }
05 	      $start += $MD_count - 1;
06 	    }
07 	
08 	    ### reconstituting shortened CIGAR string
09 	    my $new_cigar;
10 	    my $count = 0;
11 	    my $last_op;
12 	    # print "ignore adjusted: @comp_cigar\n";
13 	    foreach my $op (@comp_cigar) {
14 	      unless (defined $last_op){
15 		$last_op = $op;
16 		++$count;
17 		next;
18 	      }
19 	      if ($last_op eq $op) {
20 		++$count;
21 	      } else {
22 		$new_cigar .= "$count$last_op";
23 		$last_op = $op;
24 		$count = 1;
25 	      }
26 	    }
27 	    $new_cigar .= "$count$last_op"; # appending the last operation and count
28 	    $cigar = $new_cigar;
29 	    # print "ignore adjusted scalar: $cigar\n";
30 	  }
31 	}
32 	### printing out the methylation state of every C in the read
33 	print_individual_C_methylation_states_single_end($meth_call,$chrom,$start,$id,$strand,$index,$cigar);
34 	
35 	++$methylation_call_strings_processed; # 1 per single-end result
36       }
37     }
38   }
39 
40   ### PROCESSING PAIRED-END RESULT FILES
41   elsif ($paired) {
42 
43     ### proceeding differently now for SAM format or vanilla Bismark format files
44     if ($vanilla) {	# old vanilla Bismark paired-end output format
45       while (<IN>) {
46 	++$line_count;
47 	warn "processed line: $line_count\n" if ($line_count%500000==0);
48 
49 	### $seq here is the chromosomal sequence (to use for the repeat analysis for example)
50 	my ($id,$strand,$chrom,$start_read_1,$end_read_2,$seq_1,$meth_call_1,$seq_2,$meth_call_2,$first_read_conversion,$genome_conversion) = (split("\t"))[0,1,2,3,4,6,7,9,10,11,12,13];
51 
52 	my $index;
53 	chomp $genome_conversion;
54 	
55 	if ($first_read_conversion eq 'CT' and $genome_conversion eq 'CT') {
56 	  $index = 0;		## this is OT
57 	} elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'GA') {
58 	  $index = 2;		## this is CTOB!!!
59 	} elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'CT') {
60 	  $index = 1;		## this is CTOT!!!
61 	} elsif ($first_read_conversion eq 'CT' and $genome_conversion eq 'GA') {
62 	  $index = 3;		## this is OB
63 	} else {
64 	  die "Unexpected combination of read and genome conversion: $first_read_conversion / $genome_conversion\n";
65 	}
66 	
67 	if ($meth_call_1 and $meth_call_2) {
68 	  ### Clipping off the first <int> number of bases from the methylation call strings as specified with '--ignore <int>'
69 
70 	  if ($ignore) {
71 	    $meth_call_1 = substr($meth_call_1,$ignore,length($meth_call_1)-$ignore);
72 	
73 	    ### we also need to adjust the start and end positions of the alignments accordingly if '--ignore' was specified
74 	    $start_read_1 += $ignore;
75 	  }
76 	  if ($ignore_r2) {
77 	    $meth_call_2 = substr($meth_call_2,$ignore_r2,length($meth_call_2)-$ignore_r2);
78 	
79 	    ### we also need to adjust the start and end positions of the alignments accordingly if '--ignore_r2' was specified
80 	    $end_read_2   -= $ignore_r2;
81 	  }
82 
83 	  my $end_read_1;
84 	  my $start_read_2;
85 
86 	  if ($strand eq '+') {
87 
88 	    $end_read_1 = $start_read_1+length($meth_call_1)-1;
89 	    $start_read_2 = $end_read_2-length($meth_call_2)+1;
90 	
91 	    ## we first pass the first read which is in + orientation on the forward strand
92 	    print_individual_C_methylation_states_paired_end_files($meth_call_1,$chrom,$start_read_1,$id,'+',$index,0,0,undef,1); # the last two values are CIGAR string and read identity
93 	
94 	    # we next pass the second read which is in - orientation on the reverse strand
95 	    ### if --no_overlap was specified we also pass the end of read 1. If read 2 starts to overlap with read 1 we can stop extracting methylation calls from read 2
96 	    print_individual_C_methylation_states_paired_end_files($meth_call_2,$chrom,$end_read_2,$id,'-',$index,$no_overlap,$end_read_1,undef,2);
97 	  }
98 	  else {
99 	
00 	    $end_read_1 = $start_read_1+length($meth_call_2)-1;	# read 1 is the second reported read!
01 	    $start_read_2 = $end_read_2-length($meth_call_1)+1;	# read 2 is the first reported read!
02 
03 	    ## we first pass the first read which is in - orientation on the reverse strand
04 	    print_individual_C_methylation_states_paired_end_files($meth_call_1,$chrom,$end_read_2,$id,'-',$index,0,0,undef,1);
05 
06 	    # we next pass the second read which is in + orientation on the forward strand
07 	    ### if --no_overlap was specified we also pass the end of read 2. If read 2 starts to overlap with read 1 we will stop extracting methylation calls from read 2
08 	    print_individual_C_methylation_states_paired_end_files($meth_call_2,$chrom,$start_read_1,$id,'+',$index,$no_overlap,$start_read_2,undef,2);
09 	  }
10 	
11 	  $methylation_call_strings_processed += 2; # paired-end = 2 methylation calls
12 	}	
13       }
14     }
15     else {	      # Bismark paired-end SAM output format (default)
16       while (<IN>) {
17 	### SAM format can either start with header lines (starting with @) or start with alignments directly
18 	if (/^\@/) {	     # skipping header lines (starting with @)
19 	  warn "skipping SAM header line:\t$_";
20 	  next;
21 	}
22 	
23 	++$line_count;
24 	warn "Processed lines: $line_count\n" if ($line_count%500000==0);
25 	
26 	# example paired-end reads in SAM format (2 consecutive lines)
27 	# 1_R1/1	67	5	103172224	255	40M	=	103172417	233	AATATTTTTTTTATTTTAAAATGTGTATTGATTTAAATTT	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:4	XX:Z:4T1T24TT7	XM:Z:....h.h........................hh.......	XR:Z:CT	XG:Z:CT
28 	# 1_R1/2	131	5	103172417	255	40M	=	103172224	-233	TATTTTTTTTTAGAGTATTTTTTAATGGTTATTAGATTTT	IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII	NM:i:6	XX:Z:T5T1T9T9T7T3	XM:Z:h.....h.h.........h.........h.......h...	XR:Z:GA	XG:Z:CT
29 	
30 	#  < version 0.7.6 my ($id_1,$chrom,$start_read_1,$meth_call_1,$first_read_conversion,$genome_conversion) = (split("\t"))[0,2,3,13,14,15];
31 
32 	my ($id_1,$chrom,$start_read_1,$cigar_1) = (split("\t"))[0,2,3,5]; ### detecting the following SAM flags in case the SAM entry was shuffled by CRAM or Goby compression/decompression
33 	my $meth_call_1;
34 	my $first_read_conversion;
35 	my $genome_conversion;
36 	
37 	while ( /(XM|XR|XG):Z:([^\t]+)/g ) {
38 	  my $tag = $1;
39 	  my $value = $2;
40 
41 	  if ($tag eq "XM") {
42 	    $meth_call_1 = $value;
43 	    $meth_call_1 =~ s/\r//;
44 	  } elsif ($tag eq "XR") {
45 	    $first_read_conversion = $value;
46 	    $first_read_conversion =~ s/\r//;
47 	  } elsif ($tag eq "XG") {
48 	    $genome_conversion = $value;
49 	    $genome_conversion =~ s/\r//;
50 	  }
51 	}
52 
53 	$_ = <IN>;		# reading in the paired read
54 
55 	# < version 0.7.6 my ($id_2,$start_read_2,$meth_call_2,$second_read_conversion) = (split("\t"))[0,3,13,14];
56 	# < version 0.7.6 $meth_call_1 =~ s/^XM:Z://;
57 	# < version 0.7.6 $meth_call_2 =~ s/^XM:Z://;
58 	# < version 0.7.6 $first_read_conversion =~ s/^XR:Z://;
59 	# < version 0.7.6 $second_read_conversion =~ s/^XR:Z://;
60 
61 	my ($id_2,$start_read_2,$cigar_2) = (split("\t"))[0,3,5]; ### detecting the following SAM flags in case the SAM entry was shuffled by CRAM or Goby compression/decompression
62 
63 	my $meth_call_2;
64 	my $second_read_conversion;
65 	
66 	while ( /(XM|XR):Z:([^\t]+)/g ) {
67 	  my $tag = $1;
68 	  my $value = $2;
69 
70 	  if ($tag eq "XM") {
71 	    $meth_call_2 = $value;
72 	    $meth_call_2 =~ s/\r//;
73 	  } elsif ($tag eq "XR") {
74 	    $second_read_conversion = $value;
75 	    $second_read_conversion = s/\r//;
76 	  }
77 	}
78 	
79 	# < version 0.7.6 $genome_conversion =~ s/^XG:Z://;	
80 	chomp $genome_conversion; # in case it captured a new line character	
81 
82 	# print join ("\t",$meth_call_1,$meth_call_2,$first_read_conversion,$second_read_conversion,$genome_conversion),"\n";
83 
84 	my $index;
85 	my $strand;
86 
87 	if ($first_read_conversion eq 'CT' and $genome_conversion eq 'CT') {
88 	  $index = 0;		## this is OT
89 	  $strand = '+';
90 	} elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'CT') {
91 	  $index = 1;		## this is CTOT
92 	  $strand = '-';
93 	} elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'GA') {
94 	  $index = 2;		## this is CTOB
95 	  $strand = '+';
96 	} elsif ($first_read_conversion eq 'CT' and $genome_conversion eq 'GA') {
97 	  $index = 3;		## this is OB
98 	  $strand = '-';
99 	} else {
00 	  die "Unexpected combination of read and genome conversion: $first_read_conversion / $genome_conversion\n";
01 	}
02 
03 	### reversing the methylation call of the read that was reverse-complemented
04 	if ($strand eq '+') {
05 	  $meth_call_2 = reverse $meth_call_2;
06 	} else {
07 	  $meth_call_1 = reverse $meth_call_1;
08 	}
09 
10 	if ($meth_call_1 and $meth_call_2) {
11 
12 	  my $end_read_1;
13 	
14 	  ### READ 1
15 	  my @len_1 = split (/\D+/,$cigar_1); # storing the length per operation
16 	  my @ops_1 = split (/\d+/,$cigar_1); # storing the operation
17 	  shift @ops_1;		# remove the empty first element
18 
19 	  die "CIGAR string contained a non-matching number of lengths and operations: $cigar_1\n".join(" ",@len_1)."\n".join(" ",@ops_1)."\n" unless (scalar @len_1 == scalar @ops_1);
20 	
21 	  my @comp_cigar_1; # building an array with all CIGAR operations
22 	  foreach my $index (0..$#len_1) {
23 	    foreach (1..$len_1[$index]) {
24 	      # print  "$ops_1[$index]";
25 	      push @comp_cigar_1, $ops_1[$index];
26 	    }
27 	  }
28 	  # print "original CIGAR read 1: $cigar_1\n";
29 	  # print "original CIGAR read 1: @comp_cigar_1\n";
30 
31 	  ### READ 2
32 	  my @len_2 = split (/\D+/,$cigar_2); # storing the length per operation
33 	  my @ops_2 = split (/\d+/,$cigar_2); # storing the operation
34 	  shift @ops_2;		# remove the empty first element
35 	  die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len_2 == scalar @ops_2);
36 	  my @comp_cigar_2; # building an array with all CIGAR operations for read 2
37 	  foreach my $index (0..$#len_2) {
38 	    foreach (1..$len_2[$index]) {
39 	      # print  "$ops_2[$index]";
40 	      push @comp_cigar_2, $ops_2[$index];
41 	    }
42 	  }
43 	  # print "original CIGAR read 2: $cigar_2\n";
44 	  # print "original CIGAR read 2: @comp_cigar_2\n";
45 	
46 
47 
48 	  if ($ignore) {
49 	    ### Clipping off the first <int> number of bases from the methylation call strings as specified with '--ignore <int>' for read 1	
50 	    ### the methylation calls have already been reversed where necessary
51 	    $meth_call_1 = substr($meth_call_1,$ignore,length($meth_call_1)-$ignore);
52 
53 	    if ($strand eq '+') {
54 		
55 	      ### if the (read 1) strand information is '+', read 1 needs to be trimmed from the start
56 	      my $D_count_1 = 0; # counting all deletions that affect the ignored genomic position for read 1, i.e. Deletions and insertions
57 	      my $I_count_1 = 0;
58 	      
59 	      for (1..$ignore) {
60 		my $op = shift @comp_cigar_1; # adjusting composite CIGAR string of read 1 by removing $ignore operations from the start
61 		# print "$_ deleted $op\n";
62 		
63 		while ($op eq 'D') { # repeating this for deletions (D)
64 		  $D_count_1++;
65 		  $op = shift @comp_cigar_1;
66 		  # print "$_ deleted $op\n";
67 		}
68 		if ($op eq 'I') { # adjusting the genomic position for insertions (I)
69 		  $I_count_1++;
70 		}
71 	      }
72 		
73 	      $start_read_1 += $ignore + $D_count_1 - $I_count_1;
74 	      # print "start read 1 $start_read_1\t ignore: $ignore\t D count 1: $D_count_1\tI_count 1: $I_count_1\n";
75 
76 	      # the start position of reads mapping to the reverse strand is being adjusted further below
77 	    }
78 	    elsif ($strand eq '-') {
79 	
80 	      ### if the (read 1) strand information is '-', read 1 needs to be trimmed from the back
81 	      for (1..$ignore) {
82 		my $op = pop @comp_cigar_1; # adjusting composite CIGAR string by removing $ignore operations, here the last value of the array
83 		while ($op eq 'D') { # repeating this for deletions (D)
84 		  $op = pop @comp_cigar_1;
85 		}
86 	      }
87 	      # the start position of reads mapping to the reverse strand is being adjusted further below
88 
89 	    }
90 	  }
91 
92 	  if ($ignore_r2) {
93 	    ### Clipping off the first <int> number of bases from the methylation call string as specified with '--ignore_r2 <int>' for read 2	
94 	    ### the methylation calls have already been reversed where necessary
95 	    $meth_call_2 = substr($meth_call_2,$ignore_r2,length($meth_call_2)-$ignore_r2);
96 	
97 	    ### If we are ignoring a part of the sequence we also need to adjust the cigar string accordingly
98 	
99 	    if ($strand eq '+') {
00 		
01 	      ### if the (read 1) strand information is '+', read 2 needs to be trimmed from the back
02 
03 	      for (1..$ignore_r2) {
04 		my $op = pop @comp_cigar_2; # adjusting composite CIGAR string by removing $ignore operations, here the last value of the array
05 		while ($op eq 'D') { # repeating this for deletions (D)
06 		  $op = pop @comp_cigar_2;
07 		}
08 	      }
09 	      # the start position of reads mapping to the reverse strand is being adjusted further below
10 	    }
11 	    elsif ($strand eq '-') {
12 	
13 	      ### if the (read 1) strand information is '-', read 2 needs to be trimmed from the start
14 	      my $D_count_2 = 0; # counting all deletions that affect the ignored genomic position for read 2, i.e. Deletions and insertions
15 		      my $I_count_2 = 0;
16 
17 	      for (1..$ignore_r2) {
18 		my $op = shift @comp_cigar_2; # adjusting composite CIGAR string of read 2 by removing $ignore operations from the start
19 		# print "$_ deleted $op\n";
20 		
21 		while ($op eq 'D') { # repeating this for deletions (D)
22 		  $D_count_2++;
23 		  $op = shift @comp_cigar_2;
24 		  # print "$_ deleted $op\n";
25 		}
26 		if ($op eq 'I') { # adjusting the genomic position for insertions (I)
27 		  $I_count_2++;
28 		}
29 	      }
30 		
31 	      $start_read_2 += $ignore_r2 + $D_count_2 - $I_count_2;
32 	      # print "start read 2 $start_read_2\t ignore R2: $ignore_r2\t D count 2: $D_count_2\tI_count 2: $I_count_2\n";
33 	    }
34 	  }
35 	
36 	  if ($ignore){
37 	    ### reconstituting shortened CIGAR string 1
38 	    my $new_cigar_1;
39 	    my $count_1 = 0;
40 	    my $last_op_1;
41 	    # print "ignore adjusted CIGAR 1: @comp_cigar_1\n";
42 	    foreach my $op (@comp_cigar_1) {
43 	      unless (defined $last_op_1){
44 		$last_op_1 = $op;
45 		++$count_1;
46 		next;
47 	      }
48 	      if ($last_op_1 eq $op) {
49 		++$count_1;
50 	      } else {
51 		$new_cigar_1 .= "$count_1$last_op_1";
52 		$last_op_1 = $op;
53 		$count_1 = 1;
54 	      }
55 	    }
56 	    $new_cigar_1 .= "$count_1$last_op_1"; # appending the last operation and count
57 	    $cigar_1 = $new_cigar_1;
58 	    # print "ignore adjusted CIGAR 1 scalar: $cigar_1\n";
59 	  }
60 
61 	  if ($ignore_r2){
62 
63 	    ### reconstituting shortened CIGAR string 2
64 	    my $new_cigar_2;
65 	    my $count_2 = 0;
66 	    my $last_op_2;
67 	    # print "ignore adjusted CIGAR 2: @comp_cigar_2\n";
68 	    foreach my $op (@comp_cigar_2) {
69 	      unless (defined $last_op_2){
70 		$last_op_2 = $op;
71 		++$count_2;
72 		next;
73 	      }
74 	      if ($last_op_2 eq $op) {
75 		++$count_2;
76 	      }
77 	      else {
78 		$new_cigar_2 .= "$count_2$last_op_2";
79 		$last_op_2 = $op;
80 		$count_2 = 1;
81 	      }
82 	    }
83 	    $new_cigar_2 .= "$count_2$last_op_2"; # appending the last operation and count
84 	    $cigar_2 = $new_cigar_2;
85 	    # print "ignore_r2 adjusted CIGAR 2 scalar: $cigar_2\n";
86 	  }
87 	
88 	  ### Adjusting CIGAR string and starting position of reads in reverse orientation which we will pass to the extraction subroutine later on
89 	
90 	  if ($strand eq '+') {
91 	    ### adjusting the start position for all reads mapping to the reverse strand, in this case read 2
92 	    @comp_cigar_2  = reverse@comp_cigar_2; # the CIGAR string needs to be reversed for all reads aligning to the reverse strand, too
93 	    # print "reverse: @comp_cigar_2\n";
94 	
95 	    my $MD_count_1 = 0;
96 	    foreach (@comp_cigar_1) {
97 	      ++$MD_count_1 if ($_ eq 'M' or $_ eq 'D'); # Matching bases or deletions affect the genomic position of the 3' ends of reads, insertions don't
98 	    }
99 
00 	    my $MD_count_2 = 0;
01 	    foreach (@comp_cigar_2) {
02 	      ++$MD_count_2 if ($_ eq 'M' or $_ eq 'D'); # Matching bases or deletions affect the genomic position of the 3' ends of reads, insertions don't
03 	    }
04 
05 	    $end_read_1 = $start_read_1 + $MD_count_1 - 1;
06 	    $start_read_2 += $MD_count_2 - 1; ## Passing on the start position on the reverse strand
07 	  }
08 	  else {
09 	    ### adjusting the start position for all reads mapping to the reverse strand, in this case read 1
10 	
11 	    @comp_cigar_1  = reverse@comp_cigar_1; # the CIGAR string needs to be reversed for all reads aligning to the reverse strand, too
12 	    # print "reverse: @comp_cigar_1\n";
13 
14 	    my $MD_count_1 = 0;
15 	    foreach (@comp_cigar_1) {
16 	      ++$MD_count_1 if ($_ eq 'M' or $_ eq 'D'); # Matching bases or deletions affect the genomic position of the 3' ends of reads, insertions don't
17 	    }
18 
19 	    $end_read_1 = $start_read_1;	
20 	    $start_read_1 +=  $MD_count_1 - 1; ### Passing on the start position on the reverse strand
21 	  }
22 
23 	  if ($strand eq '+') {
24 	    ## we first pass the first read which is in + orientation on the forward strand; the last value is the read identity
25 	    print_individual_C_methylation_states_paired_end_files($meth_call_1,$chrom,$start_read_1,$id_1,'+',$index,0,0,$cigar_1,1);
26 	
27 	    # we next pass the second read which is in - orientation on the reverse strand
28 	    ### if --no_overlap was specified we also pass the end of read 1. If read 2 starts to overlap with read 1 we can stop extracting methylation calls from read 2
29 	    print_individual_C_methylation_states_paired_end_files($meth_call_2,$chrom,$start_read_2,$id_2,'-',$index,$no_overlap,$end_read_1,$cigar_2,2);
30 	  } else {
31 	    ## we first pass the first read which is in - orientation on the reverse strand
32 	    print_individual_C_methylation_states_paired_end_files($meth_call_1,$chrom,$start_read_1,$id_1,'-',$index,0,0,$cigar_1,1);
33 	
34 	    # we next pass the second read which is in + orientation on the forward strand
35 	    ### if --no_overlap was specified we also pass the end of read 1. If read 2 starts to overlap with read 1 we will stop extracting methylation calls from read 2
36 	    print_individual_C_methylation_states_paired_end_files($meth_call_2,$chrom,$start_read_2,$id_2,'+',$index,$no_overlap,$end_read_1,$cigar_2,2);
37 	  }
38 	
39 	  $methylation_call_strings_processed += 2; # paired-end = 2 methylation calls
40 	}	
41       }
42     }
43   } else {
44     die "Single-end or paired-end reads not specified properly\n";
45   }
46 
47   warn "\n\nProcessed $line_count lines from $filename in total\n";
48   warn "Total number of methylation call strings processed: $methylation_call_strings_processed\n\n";
49   if ($report) {
50     print REPORT "\n\nProcessed $line_count lines from $filename in total\n";
51     print REPORT "Total number of methylation call strings processed: $methylation_call_strings_processed\n\n";
52   }
53   print_splitting_report ();
54 }
55 
56 
57 
58 sub print_splitting_report{
59 
60   ### Calculating methylation percentages if applicable
61 
62   my $percent_meCpG;
63   if (($counting{total_meCpG_count}+$counting{total_unmethylated_CpG_count}) > 0){
64     $percent_meCpG = sprintf("%.1f",100*$counting{total_meCpG_count}/($counting{total_meCpG_count}+$counting{total_unmethylated_CpG_count}));
65   }
66 
67   my $percent_meCHG;
68   if (($counting{total_meCHG_count}+$counting{total_unmethylated_CHG_count}) > 0){
69     $percent_meCHG = sprintf("%.1f",100*$counting{total_meCHG_count}/($counting{total_meCHG_count}+$counting{total_unmethylated_CHG_count}));
70   }
71 
72   my $percent_meCHH;
73   if (($counting{total_meCHH_count}+$counting{total_unmethylated_CHH_count}) > 0){
74     $percent_meCHH = sprintf("%.1f",100*$counting{total_meCHH_count}/($counting{total_meCHH_count}+$counting{total_unmethylated_CHH_count}));
75   }
76 
77   my $percent_non_CpG_methylation;
78   if ($merge_non_CpG){
79     if ( ($counting{total_meCHH_count}+$counting{total_unmethylated_CHH_count}+$counting{total_meCHG_count}+$counting{total_unmethylated_CHG_count}) > 0){
80       $percent_non_CpG_methylation = sprintf("%.1f",100* ( $counting{total_meCHH_count}+$counting{total_meCHG_count} ) / ( $counting{total_meCHH_count}+$counting{total_unmethylated_CHH_count}+$counting{total_meCHG_count}+$counting{total_unmethylated_CHG_count} ) );
81     }
82   }
83 
84   if ($report){
85     ### detailed information about Cs analysed
86     print REPORT "Final Cytosine Methylation Report\n",'='x33,"\n";
87 
88     my $total_number_of_C = $counting{total_meCHG_count}+$counting{total_meCHH_count}+$counting{total_meCpG_count}+$counting{total_unmethylated_CHG_count}+$counting{total_unmethylated_CHH_count}+$counting{total_unmethylated_CpG_count};
89     print REPORT "Total number of C's analysed:\t$total_number_of_C\n\n";
90 
91     print REPORT "Total methylated C's in CpG context:\t$counting{total_meCpG_count}\n";
92     print REPORT "Total methylated C's in CHG context:\t$counting{total_meCHG_count}\n";
93     print REPORT "Total methylated C's in CHH context:\t$counting{total_meCHH_count}\n\n";
94 
95     print REPORT "Total C to T conversions in CpG context:\t$counting{total_unmethylated_CpG_count}\n";
96     print REPORT "Total C to T conversions in CHG context:\t$counting{total_unmethylated_CHG_count}\n";
97     print REPORT "Total C to T conversions in CHH context:\t$counting{total_unmethylated_CHH_count}\n\n";
98 
99     ### calculating methylated CpG percentage if applicable
00     if ($percent_meCpG){
01       print REPORT "C methylated in CpG context:\t${percent_meCpG}%\n";
02     }
03     else{
04       print REPORT "Can't determine percentage of methylated Cs in CpG context if value was 0\n";
05     }
06 
07     ### 2-Context Output
08     if ($merge_non_CpG){
09       if ($percent_non_CpG_methylation){
10 	print REPORT "C methylated in non-CpG context:\t${percent_non_CpG_methylation}%\n\n\n";
11       }
12       else{
13 	print REPORT "Can't determine percentage of methylated Cs in non-CpG context if value was 0\n\n\n";
14       }
15     }
16 
17     ### 3 Context Output
18     else{
19       ### calculating methylated CHG percentage if applicable
20       if ($percent_meCHG){
21 	print REPORT "C methylated in CHG context:\t${percent_meCHG}%\n";
22       }
23       else{
24 	print REPORT "Can't determine percentage of methylated Cs in CHG context if value was 0\n";
25       }
26 
27       ### calculating methylated CHH percentage if applicable
28       if ($percent_meCHH){
29 	print REPORT "C methylated in CHH context:\t${percent_meCHH}%\n\n\n";
30       }
31       else{
32 	print REPORT "Can't determine percentage of methylated Cs in CHH context if value was 0\n\n\n";
33       }
34     }
35   }
36 
37   ### detailed information about Cs analysed for on-screen report
38   print "Final Cytosine Methylation Report\n",'='x33,"\n";
39 
40   my $total_number_of_C = $counting{total_meCHG_count}+$counting{total_meCHH_count}+$counting{total_meCpG_count}+$counting{total_unmethylated_CHG_count}+$counting{total_unmethylated_CHH_count}+$counting{total_unmethylated_CpG_count};
41   print "Total number of C's analysed:\t$total_number_of_C\n\n";
42 
43   print "Total methylated C's in CpG context:\t$counting{total_meCpG_count}\n";
44   print "Total methylated C's in CHG context:\t$counting{total_meCHG_count}\n";
45   print "Total methylated C's in CHH context:\t$counting{total_meCHH_count}\n\n";
46 
47   print "Total C to T conversions in CpG context:\t$counting{total_unmethylated_CpG_count}\n";
48   print "Total C to T conversions in CHG context:\t$counting{total_unmethylated_CHG_count}\n";
49   print "Total C to T conversions in CHH context:\t$counting{total_unmethylated_CHH_count}\n\n";
50 
51   ### printing methylated CpG percentage if applicable
52   if ($percent_meCpG){
53     print "C methylated in CpG context:\t${percent_meCpG}%\n";
54   }
55   else{
56     print "Can't determine percentage of methylated Cs in CpG context if value was 0\n";
57   }
58 
59   ### 2-Context Output
60   if ($merge_non_CpG){
61     if ($percent_non_CpG_methylation){
62       print "C methylated in non-CpG context:\t${percent_non_CpG_methylation}%\n\n\n";
63     }
64     else{
65       print "Can't determine percentage of methylated Cs in non-CpG context if value was 0\n\n\n";
66     }
67   }
68 
69   ### 3-Context Output
70   else{
71     ### printing methylated CHG percentage if applicable
72     if ($percent_meCHG){
73       print "C methylated in CHG context:\t${percent_meCHG}%\n";
74     }
75     else{
76       print "Can't determine percentage of methylated Cs in CHG context if value was 0\n";
77     }
78 
79     ### printing methylated CHH percentage if applicable
80     if ($percent_meCHH){
81       print "C methylated in CHH context:\t${percent_meCHH}%\n\n\n";
82     }
83     else{
84       print "Can't determine percentage of methylated Cs in CHH context if value was 0\n\n\n";
85     }
86   }
87 }
88 
89 
90 
91 sub print_individual_C_methylation_states_paired_end_files{
92 
93   my ($meth_call,$chrom,$start,$id,$strand,$filehandle_index,$no_overlap,$end_read_1,$cigar,$read_identity) = @_;
94 
95   ### we will use the read identity for the M-bias plot to discriminate read 1 and read 2
96   die "Read identity was neither 1 nor 2: $read_identity\n\n" unless ($read_identity == 1 or $read_identity == 2);
97 
98   my @methylation_calls = split(//,$meth_call);
99 
00   #################################################################
01   ### . for bases not involving cytosines                       ###
02   ### X for methylated C in CHG context (was protected)         ###
03   ### x for not methylated C in CHG context (was converted)     ###
04   ### H for methylated C in CHH context (was protected)         ###
05   ### h for not methylated C in CHH context (was converted)     ###
06   ### Z for methylated C in CpG context (was protected)         ###
07   ### z for not methylated C in CpG context (was converted)     ###
08   ### U for methylated C in Unknown context (was protected)     ###
09   ### u for not methylated C in Unknown context (was converted) ###
10   #################################################################
11 
12   my $methyl_CHG_count = 0;
13   my $methyl_CHH_count = 0;
14   my $methyl_CpG_count = 0;
15   my $unmethylated_CHG_count = 0;
16   my $unmethylated_CHH_count = 0;
17   my $unmethylated_CpG_count = 0;
18 
19   my $pos_offset = 0; # this is only relevant for SAM reads with insertions or deletions
20   my $cigar_offset = 0; # again, this is only relevant for SAM reads containing indels
21   my @comp_cigar;
22 
23   ### Checking whether the CIGAR string is a linear genomic match or whether if requires indel processing
24   if ($cigar =~ /^\d+M$/){
25     # this check speeds up the extraction process by up to 60%!!!
26   }
27   else{ # parsing CIGAR string
28     my @len;
29     my @ops;
30     @len = split (/\D+/,$cigar); # storing the length per operation
31     @ops = split (/\d+/,$cigar); # storing the operation
32     shift @ops; # remove the empty first element
33 
34     die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len == scalar @ops);
35 
36     foreach my $index (0..$#len){
37       foreach (1..$len[$index]){
38 	# print  "$ops[$index]";
39 	push @comp_cigar, $ops[$index];
40       }
41     }
42     # warn "\nDetected CIGAR string: $cigar\n";
43     # warn "Length of methylation call: ",length $meth_call,"\n";
44     # warn "number of operations: ",scalar @ops,"\n";
45     # warn "number of length digits: ",scalar @len,"\n\n";
46     # print @comp_cigar,"\n";
47     # print "$meth_call\n\n";
48     # sleep (1);
49   }
50 
51   if ($strand eq '-') {
52 
53     ### the  CIGAR string needs to be reversed, the methylation call has already been reversed above
54     if (@comp_cigar){
55       @comp_cigar  = reverse@comp_cigar; # the CIGAR string needs to be reversed for all reads aligning to the reverse strand, too
56     }
57     #  print "reverse CIGAR string: @comp_cigar\n";
58 
59     ### the start position of paired-end files has already been corrected, see above
60   }
61 
62   ### THIS IS AN OPTIONAL 2-CONTEXT (CpG and non-CpG) SECTION IF --merge_non_CpG was specified
63 
64   if ($merge_non_CpG) {
65     if ($no_overlap) { # this has to be read 2...
66 
67       ### single-file CpG and non-CpG context output
68       if ($full) {
69 	if ($strand eq '+') {
70 	  for my $index (0..$#methylation_calls) {
71 	
72 	    if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
73 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
74 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
75 	      $cigar_offset += $cigar_mod;
76 	      $pos_offset += $pos_mod;
77 	    }
78 	
79 	    ### Returning as soon as the methylation calls start overlapping
80 	    if ($start+$index+$pos_offset >= $end_read_1) {
81 	      return;
82 	    }
83 
84 	    if ($methylation_calls[$index] eq 'X') {
85 	      $counting{total_meCHG_count}++;
86 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
87 	      if ($read_identity == 1){
88 		$mbias_1{CHG}->{$index+1}->{meth}++;
89 	      }
90 	      else{
91 		$mbias_2{CHG}->{$index+1}->{meth}++;
92 	      }
93 	    }
94 	    elsif ($methylation_calls[$index] eq 'x') {
95 	      $counting{total_unmethylated_CHG_count}++;
96 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
97 	      if ($read_identity == 1){
98 		$mbias_1{CHG}->{$index+1}->{un}++;
99 	      }
00 	      else{
01 		$mbias_2{CHG}->{$index+1}->{un}++;
02 	      }
03 	    }
04 	    elsif ($methylation_calls[$index] eq 'Z') {
05 	      $counting{total_meCpG_count}++;
06 	      print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
07 	      if ($read_identity == 1){
08 		$mbias_1{CpG}->{$index+1}->{meth}++;
09 	      }
10 	      else{
11 		$mbias_2{CpG}->{$index+1}->{meth}++;
12 	      }
13 	    }
14 	    elsif ($methylation_calls[$index] eq 'z') {
15 	      $counting{total_unmethylated_CpG_count}++;
16 	      print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
17 	      if ($read_identity == 1){
18 		$mbias_1{CpG}->{$index+1}->{un}++;
19 	      }
20 	      else{
21 		$mbias_2{CpG}->{$index+1}->{un}++;
22 	      }
23 	    }
24 	    elsif ($methylation_calls[$index] eq 'H') {
25 	      $counting{total_meCHH_count}++;
26 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
27 	      if ($read_identity == 1){
28 		$mbias_1{CHH}->{$index+1}->{meth}++;
29 	      }
30 	      else{
31 		$mbias_2{CHH}->{$index+1}->{meth}++;
32 	      }
33 	    }
34 	    elsif ($methylation_calls[$index] eq 'h') {
35 	      $counting{total_unmethylated_CHH_count}++;
36 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
37 	      if ($read_identity == 1){
38 		$mbias_1{CHH}->{$index+1}->{un}++;
39 	      }
40 	      else{
41 		$mbias_2{CHH}->{$index+1}->{un}++;
42 	      }
43 	    }
44 	    elsif ($methylation_calls[$index] eq '.'){}
45 	    elsif (lc$methylation_calls[$index] eq 'u'){}
46 	    else{
47 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n" unless($mbias_only);
48 	    }
49 	  }
50 	}
51 	elsif ($strand eq '-') {
52 	  for my $index (0..$#methylation_calls) {
53 	
54 	    if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
55 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
56 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
57 	      $cigar_offset += $cigar_mod;
58 	      $pos_offset += $pos_mod;
59 	    }
60 	
61 	    ### Returning as soon as the methylation calls start overlapping
62 	    if ($start-$index+$pos_offset <= $end_read_1) {
63 	      return;
64 	    }
65 	
66 	    if ($methylation_calls[$index] eq 'X') {
67 	      $counting{total_meCHG_count}++;
68 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
69 	      if ($read_identity == 1){
70 		$mbias_1{CHG}->{$index+1}->{meth}++;
71 	      }
72 	      else{
73 		$mbias_2{CHG}->{$index+1}->{meth}++;
74 	      }
75 	    }
76 	    elsif ($methylation_calls[$index] eq 'x') {
77 	      $counting{total_unmethylated_CHG_count}++;
78 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
79 	      if ($read_identity == 1){
80 		$mbias_1{CHG}->{$index+1}->{un}++;
81 	      }
82 	      else{
83 		$mbias_2{CHG}->{$index+1}->{un}++;
84 	      }
85 	    }
86 	    elsif ($methylation_calls[$index] eq 'Z') {
87 	      $counting{total_meCpG_count}++;
88 	      print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
89 	      if ($read_identity == 1){
90 		$mbias_1{CpG}->{$index+1}->{meth}++;
91 	      }
92 	      else{
93 		$mbias_2{CpG}->{$index+1}->{meth}++;
94 	      }
95 	    }
96 	    elsif ($methylation_calls[$index] eq 'z') {
97 	      $counting{total_unmethylated_CpG_count}++;
98 	      print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
99 	      if ($read_identity == 1){
00 		$mbias_1{CpG}->{$index+1}->{un}++;
01 	      }
02 	      else{
03 		$mbias_2{CpG}->{$index+1}->{un}++;
04 	      }
05 	    }
06 	    elsif ($methylation_calls[$index] eq 'H') {
07 	      $counting{total_meCHH_count}++;
08 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
09 	      if ($read_identity == 1){
10 		$mbias_1{CHH}->{$index+1}->{meth}++;
11 	      }
12 	      else{
13 		$mbias_2{CHH}->{$index+1}->{meth}++;
14 	      }
15 	    }
16 	    elsif ($methylation_calls[$index] eq 'h') {
17 	      $counting{total_unmethylated_CHH_count}++;
18 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
19 	      if ($read_identity == 1){
20 		$mbias_1{CHH}->{$index+1}->{un}++;
21 	      }
22 	      else{
23 		$mbias_2{CHH}->{$index+1}->{un}++;
24 	      }
25 	    }
26 	    elsif ($methylation_calls[$index] eq '.') {}
27 	    elsif (lc$methylation_calls[$index] eq 'u'){}
28 	    else{
29 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n" unless($mbias_only);
30 	    }
31 	  }
32 	} else {
33 	  die "The read orientation was neither + nor -: '$strand'\n";
34 	}
35       }
36 
37       ### strand-specific methylation output
38       else {
39 	if ($strand eq '+') {
40 	  for my $index (0..$#methylation_calls) {
41 
42 	    if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
43 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
44 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
45 	      $cigar_offset += $cigar_mod;
46 	      $pos_offset += $pos_mod;
47 	    }
48 
49 	    ### Returning as soon as the methylation calls start overlapping
50 	    if ($start+$index+$pos_offset >= $end_read_1) {
51 	      return;
52 	    }
53 
54 	    if ($methylation_calls[$index] eq 'X') {
55 	      $counting{total_meCHG_count}++;
56 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
57 	      if ($read_identity == 1){
58 		$mbias_1{CHG}->{$index+1}->{meth}++;
59 	      }
60 	      else{
61 		$mbias_2{CHG}->{$index+1}->{meth}++;
62 	      }
63 	    }
64 	    elsif ($methylation_calls[$index] eq 'x') {
65 	      $counting{total_unmethylated_CHG_count}++;
66 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
67 	      if ($read_identity == 1){
68 		$mbias_1{CHG}->{$index+1}->{un}++;
69 	      }
70 	      else{
71 		$mbias_2{CHG}->{$index+1}->{un}++;
72 	      }
73 	    }
74 	    elsif ($methylation_calls[$index] eq 'Z') {
75 	      $counting{total_meCpG_count}++;
76 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
77 	      if ($read_identity == 1){
78 		$mbias_1{CpG}->{$index+1}->{meth}++;
79 	      }
80 	      else{
81 		$mbias_2{CpG}->{$index+1}->{meth}++;
82 	      }
83 	    }
84 	    elsif ($methylation_calls[$index] eq 'z') {
85 	      $counting{total_unmethylated_CpG_count}++;
86 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
87 	      if ($read_identity == 1){
88 		$mbias_1{CpG}->{$index+1}->{un}++;
89 	      }
90 	      else{
91 		$mbias_2{CpG}->{$index+1}->{un}++;
92 	      }
93 	    }
94 	    elsif ($methylation_calls[$index] eq 'H') {
95 	      $counting{total_meCHH_count}++;
96 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
97 	      if ($read_identity == 1){
98 		$mbias_1{CHH}->{$index+1}->{meth}++;
99 	      }
00 	      else{
01 		$mbias_2{CHH}->{$index+1}->{meth}++;
02 	      }
03 	    }
04 	    elsif ($methylation_calls[$index] eq 'h') {
05 	      $counting{total_unmethylated_CHH_count}++;
06 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
07 	      if ($read_identity == 1){
08 		$mbias_1{CHH}->{$index+1}->{un}++;
09 	      }
10 	      else{
11 		$mbias_2{CHH}->{$index+1}->{un}++;
12 	      }
13 	    }
14 	    elsif ($methylation_calls[$index] eq '.') {}
15 	    elsif (lc$methylation_calls[$index] eq 'u'){}
16 	    else{
17 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
18 	    }	
19 	  }
20 	} elsif ($strand eq '-') {
21 	  for my $index (0..$#methylation_calls) {
22 
23 	    if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
24 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
25 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
26 	      $cigar_offset += $cigar_mod;
27 	      $pos_offset += $pos_mod;
28 	    }
29 
30 	    ### Returning as soon as the methylation calls start overlapping
31 	    if ($start-$index+$pos_offset <= $end_read_1) {
32 	      return;
33 	    }
34 
35 	    if ($methylation_calls[$index] eq 'X') {
36 	      $counting{total_meCHG_count}++;
37 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
38 	      if ($read_identity == 1){
39 		$mbias_1{CHG}->{$index+1}->{meth}++;
40 	      }
41 	      else{
42 		$mbias_2{CHG}->{$index+1}->{meth}++;
43 	      }
44 	    }
45 	    elsif ($methylation_calls[$index] eq 'x') {
46 	      $counting{total_unmethylated_CHG_count}++;
47 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
48 	      if ($read_identity == 1){
49 		$mbias_1{CHG}->{$index+1}->{un}++;
50 	      }
51 	      else{
52 		$mbias_2{CHG}->{$index+1}->{un}++;
53 	      }
54 	    }
55 	    elsif ($methylation_calls[$index] eq 'Z') {
56 	      $counting{total_meCpG_count}++;
57 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
58 	      if ($read_identity == 1){
59 		$mbias_1{CpG}->{$index+1}->{meth}++;
60 	      }
61 	      else{
62 		$mbias_2{CpG}->{$index+1}->{meth}++;
63 	      }
64 	    }
65 	    elsif ($methylation_calls[$index] eq 'z') {
66 	      $counting{total_unmethylated_CpG_count}++;
67 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
68 	      if ($read_identity == 1){
69 		$mbias_1{CpG}->{$index+1}->{un}++;
70 	      }
71 	      else{
72 		$mbias_2{CpG}->{$index+1}->{un}++;
73 	      }
74 	    }
75 	    elsif ($methylation_calls[$index] eq 'H') {
76 	      $counting{total_meCHH_count}++;
77 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
78 	      if ($read_identity == 1){
79 		$mbias_1{CHH}->{$index+1}->{meth}++;
80 	      }
81 	      else{
82 		$mbias_2{CHH}->{$index+1}->{meth}++;
83 	      }
84 	    }
85 	    elsif ($methylation_calls[$index] eq 'h') {
86 	      $counting{total_unmethylated_CHH_count}++;
87 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
88 	      if ($read_identity == 1){
89 		$mbias_1{CHH}->{$index+1}->{un}++;
90 	      }
91 	      else{
92 		$mbias_2{CHH}->{$index+1}->{un}++;
93 	      }
94 	    }
95 	    elsif ($methylation_calls[$index] eq '.') {}
96 	    elsif (lc$methylation_calls[$index] eq 'u'){}
97 	    else{
98 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
99 	    }
00 	  }
01 	} else {
02 	  die "The strand orientation was neither + nor -: '$strand'/n";
03 	}
04       }
05     }
06 
07     ### this is the default paired-end procedure allowing overlaps and using every single C position
08     ### Still within the 2-CONTEXT ONLY optional section
09     else {
10       ### single-file CpG and non-CpG context output
11       if ($full) {
12 	if ($strand eq '+') {
13 	  for my $index (0..$#methylation_calls) {
14 
15 	    if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
16 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
17 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
18 	      $cigar_offset += $cigar_mod;
19 	      $pos_offset += $pos_mod;
20 	    }
21 
22 	    if ($methylation_calls[$index] eq 'X') {
23 	      $counting{total_meCHG_count}++;
24 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
25 	      if ($read_identity == 1){
26 		$mbias_1{CHG}->{$index+1}->{meth}++;
27 	      }
28 	      else{
29 		$mbias_2{CHG}->{$index+1}->{meth}++;
30 	      }
31 	    }
32 	    elsif ($methylation_calls[$index] eq 'x') {
33 	      $counting{total_unmethylated_CHG_count}++;
34 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
35 	      if ($read_identity == 1){
36 		$mbias_1{CHG}->{$index+1}->{un}++;
37 	      }
38 	      else{
39 		$mbias_2{CHG}->{$index+1}->{un}++;
40 	      }
41 	    }
42 	    elsif ($methylation_calls[$index] eq 'Z') {
43 	      $counting{total_meCpG_count}++;
44 	      print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
45 	      if ($read_identity == 1){
46 		$mbias_1{CpG}->{$index+1}->{meth}++;
47 	      }
48 	      else{
49 		$mbias_2{CpG}->{$index+1}->{meth}++;
50 	      }
51 	    }
52 	    elsif ($methylation_calls[$index] eq 'z') {
53 	      $counting{total_unmethylated_CpG_count}++;
54 	      print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
55 	      if ($read_identity == 1){
56 		$mbias_1{CpG}->{$index+1}->{un}++;
57 	      }
58 	      else{
59 		$mbias_2{CpG}->{$index+1}->{un}++;
60 	      }
61 	    }
62 	    elsif ($methylation_calls[$index] eq 'H') {
63 	      $counting{total_meCHH_count}++;
64 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
65 	      if ($read_identity == 1){
66 		$mbias_1{CHH}->{$index+1}->{meth}++;
67 	      }
68 	      else{
69 		$mbias_2{CHH}->{$index+1}->{meth}++;
70 	      }
71 	    }
72 	    elsif ($methylation_calls[$index] eq 'h') {
73 	      $counting{total_unmethylated_CHH_count}++;
74 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
75 	      if ($read_identity == 1){
76 		$mbias_1{CHH}->{$index+1}->{un}++;
77 	      }
78 	      else{
79 		$mbias_2{CHH}->{$index+1}->{un}++;
80 	      }
81 	    }
82 	    elsif ($methylation_calls[$index] eq '.') {}
83 	    elsif (lc$methylation_calls[$index] eq 'u'){}
84 	    else{
85 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n" unless($mbias_only);
86 	    }
87 	  }
88 	} elsif ($strand eq '-') {
89 	  for my $index (0..$#methylation_calls) {
90 
91 	    if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
92 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
93 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
94 	      $cigar_offset += $cigar_mod;
95 	      $pos_offset += $pos_mod;
96 	    }
97 
98 	    if ($methylation_calls[$index] eq 'X') {
99 	      $counting{total_meCHG_count}++;
00 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
01 	       if ($read_identity == 1){
02 		$mbias_1{CHG}->{$index+1}->{meth}++;
03 	      }
04 	      else{
05 		$mbias_2{CHG}->{$index+1}->{meth}++;
06 	      }
07 	    }
08 	    elsif ($methylation_calls[$index] eq 'x') {
09 	      $counting{total_unmethylated_CHG_count}++;
10 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
11 	      if ($read_identity == 1){
12 		$mbias_1{CHG}->{$index+1}->{un}++;
13 	      }
14 	      else{
15 		$mbias_2{CHG}->{$index+1}->{un}++;
16 	      }
17 	    }
18 	    elsif ($methylation_calls[$index] eq 'Z') {
19 	      $counting{total_meCpG_count}++;
20 	      print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
21 	      if ($read_identity == 1){
22 		$mbias_1{CpG}->{$index+1}->{meth}++;
23 	      }
24 	      else{
25 		$mbias_2{CpG}->{$index+1}->{meth}++;
26 	      }
27 	    }
28 	    elsif ($methylation_calls[$index] eq 'z') {
29 	      $counting{total_unmethylated_CpG_count}++;
30 	      print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
31 	      if ($read_identity == 1){
32 		$mbias_1{CpG}->{$index+1}->{un}++;
33 	      }
34 	      else{
35 		$mbias_2{CpG}->{$index+1}->{un}++;
36 	      }
37 	    }
38 	    elsif ($methylation_calls[$index] eq 'H') {
39 	      $counting{total_meCHH_count}++;
40 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
41 	      if ($read_identity == 1){
42 		$mbias_1{CHH}->{$index+1}->{meth}++;
43 	      }
44 	      else{
45 		$mbias_2{CHH}->{$index+1}->{meth}++;
46 	      }
47 	    }
48 	    elsif ($methylation_calls[$index] eq 'h') {
49 	      $counting{total_unmethylated_CHH_count}++;
50 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
51 	      if ($read_identity == 1){
52 		$mbias_1{CHH}->{$index+1}->{un}++;
53 	      }
54 	      else{
55 		$mbias_2{CHH}->{$index+1}->{un}++;
56 	      }
57 	    }
58 	    elsif ($methylation_calls[$index] eq '.') {}
59 	    elsif (lc$methylation_calls[$index] eq 'u'){}
60 	    else{
61 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n" unless($mbias_only);
62 	    }
63 	  }
64 	} else {
65 	  die "The strand orientation as neither + nor -: '$strand'\n";
66 	}
67       }
68 
69       ### strand-specific methylation output
70       ### still within the 2-CONTEXT optional section
71       else {
72 	if ($strand eq '+') {
73 	  for my $index (0..$#methylation_calls) {
74 
75 	    if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
76 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
77 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
78 	      $cigar_offset += $cigar_mod;
79 	      $pos_offset += $pos_mod;
80 	    }
81 
82 	    if ($methylation_calls[$index] eq 'X') {
83 	      $counting{total_meCHG_count}++;
84 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
85 	       if ($read_identity == 1){
86 		$mbias_1{CHG}->{$index+1}->{meth}++;
87 	      }
88 	      else{
89 		$mbias_2{CHG}->{$index+1}->{meth}++;
90 	      }
91 	    }
92 	    elsif ($methylation_calls[$index] eq 'x') {
93 	      $counting{total_unmethylated_CHG_count}++;
94 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
95 	      if ($read_identity == 1){
96 		$mbias_1{CHG}->{$index+1}->{un}++;
97 	      }
98 	      else{
99 		$mbias_2{CHG}->{$index+1}->{un}++;
00 	      }
01 	    }
02 	    elsif ($methylation_calls[$index] eq 'Z') {
03 	      $counting{total_meCpG_count}++;
04 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
05 	      if ($read_identity == 1){
06 		$mbias_1{CpG}->{$index+1}->{meth}++;
07 	      }
08 	      else{
09 		$mbias_2{CpG}->{$index+1}->{meth}++;
10 	      }
11 	    }
12 	    elsif ($methylation_calls[$index] eq 'z') {
13 	      $counting{total_unmethylated_CpG_count}++;
14 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
15 	      if ($read_identity == 1){
16 		$mbias_1{CpG}->{$index+1}->{un}++;
17 	      }
18 	      else{
19 		$mbias_2{CpG}->{$index+1}->{un}++;
20 	      }
21 	    }
22 	    elsif ($methylation_calls[$index] eq 'H') {
23 	      $counting{total_meCHH_count}++;
24 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
25 	      if ($read_identity == 1){
26 		$mbias_1{CHH}->{$index+1}->{meth}++;
27 	      }
28 	      else{
29 		$mbias_2{CHH}->{$index+1}->{meth}++;
30 	      }
31 	    }
32 	    elsif ($methylation_calls[$index] eq 'h') {
33 	      $counting{total_unmethylated_CHH_count}++;
34 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
35 	      if ($read_identity == 1){
36 		$mbias_1{CHH}->{$index+1}->{un}++;
37 	      }
38 	      else{
39 		$mbias_2{CHH}->{$index+1}->{un}++;
40 	      }
41 	    }
42 	    elsif ($methylation_calls[$index] eq '.') {}
43 	    elsif (lc$methylation_calls[$index] eq 'u'){}
44 	    else{
45 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
46 	    }
47 	  }
48 	} elsif ($strand eq '-') {
49 	  for my $index (0..$#methylation_calls) {
50 
51 	    if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
52 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
53 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
54 	      $cigar_offset += $cigar_mod;
55 	      $pos_offset += $pos_mod;
56 	    }
57 
58 	    if ($methylation_calls[$index] eq 'X') {
59 	      $counting{total_meCHG_count}++;
60 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
61 	      if ($read_identity == 1){
62 		$mbias_1{CHG}->{$index+1}->{meth}++;
63 	      }
64 	      else{
65 		$mbias_2{CHG}->{$index+1}->{meth}++;
66 	      }
67 	    }
68 	    elsif ($methylation_calls[$index] eq 'x') {
69 	      $counting{total_unmethylated_CHG_count}++;
70 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
71 	      if ($read_identity == 1){
72 		$mbias_1{CHG}->{$index+1}->{un}++;
73 	      }
74 	      else{
75 		$mbias_2{CHG}->{$index+1}->{un}++;
76 	      }
77 	    }
78 	    elsif ($methylation_calls[$index] eq 'Z') {
79 	      $counting{total_meCpG_count}++;
80 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
81 	      if ($read_identity == 1){
82 		$mbias_1{CpG}->{$index+1}->{meth}++;
83 	      }
84 	      else{
85 		$mbias_2{CpG}->{$index+1}->{meth}++;
86 	      }
87 	    }
88 	    elsif ($methylation_calls[$index] eq 'z') {
89 	      $counting{total_unmethylated_CpG_count}++;
90 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
91 	      if ($read_identity == 1){
92 		$mbias_1{CpG}->{$index+1}->{un}++;
93 	      }
94 	      else{
95 		$mbias_2{CpG}->{$index+1}->{un}++;
96 	      }
97 	    }
98 	    elsif ($methylation_calls[$index] eq 'H') {
99 	      $counting{total_meCHH_count}++;
00 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
01 	      if ($read_identity == 1){
02 		$mbias_1{CHH}->{$index+1}->{meth}++;
03 	      }
04 	      else{
05 		$mbias_2{CHH}->{$index+1}->{meth}++;
06 	      }
07 	    }
08 	    elsif ($methylation_calls[$index] eq 'h') {
09 	      $counting{total_unmethylated_CHH_count}++;
10 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
11 	      if ($read_identity == 1){
12 		$mbias_1{CHH}->{$index+1}->{un}++;
13 	      }
14 	      else{
15 		$mbias_2{CHH}->{$index+1}->{un}++;
16 	      }
17 	    }
18 	    elsif ($methylation_calls[$index] eq '.') {}
19 	    elsif (lc$methylation_calls[$index] eq 'u'){}
20 	    else{
21 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
22 	    }
23 	  }
24 	} else {
25 	  die "The strand orientation as neither + nor -: '$strand'\n";
26 	}
27       }
28     }
29   }
30 
31   ############################################
32   ### THIS IS THE DEFAULT 3-CONTEXT OUTPUT ###
33   ############################################
34 
35   elsif ($no_overlap) {
36     ### single-file CpG, CHG and CHH context output
37     if ($full) {
38       if ($strand eq '+') {
39 	for my $index (0..$#methylation_calls) {
40 	
41 	  if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
42 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
43 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
44 	    $cigar_offset += $cigar_mod;
45 	    $pos_offset += $pos_mod;
46 	  }
47 
48 	  ### Returning as soon as the methylation calls start overlapping
49 	  if ($start+$index+$pos_offset >= $end_read_1) {
50 	    return;
51 	  }
52 	
53 	  if ($methylation_calls[$index] eq 'X') {
54 	    $counting{total_meCHG_count}++;
55 	    print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
56 	    if ($read_identity == 1){
57 	      $mbias_1{CHG}->{$index+1}->{meth}++;
58 	    }
59 	    else{
60 	      $mbias_2{CHG}->{$index+1}->{meth}++;
61 	    }
62 	  }
63 	  elsif ($methylation_calls[$index] eq 'x') {
64 	    $counting{total_unmethylated_CHG_count}++;
65 	    print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
66 	    if ($read_identity == 1){
67 	      $mbias_1{CHG}->{$index+1}->{un}++;
68 	    }
69 	    else{
70 	      $mbias_2{CHG}->{$index+1}->{un}++;
71 	    }
72 	  }
73 	  elsif ($methylation_calls[$index] eq 'Z') {
74 	    $counting{total_meCpG_count}++;
75 	    print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
76 	    if ($read_identity == 1){
77 	      $mbias_1{CpG}->{$index+1}->{meth}++;
78 	    }
79 	    else{
80 	      $mbias_2{CpG}->{$index+1}->{meth}++;
81 	    }
82 	  }
83 	  elsif ($methylation_calls[$index] eq 'z') {
84 	    $counting{total_unmethylated_CpG_count}++;
85 	    print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
86 	    if ($read_identity == 1){
87 	      $mbias_1{CpG}->{$index+1}->{un}++;
88 	    }
89 	    else{
90 	      $mbias_2{CpG}->{$index+1}->{un}++;
91 	    }
92 	  }
93 	  elsif ($methylation_calls[$index] eq 'H') {
94 	    $counting{total_meCHH_count}++;
95 	    print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
96 	    if ($read_identity == 1){
97 	      $mbias_1{CHH}->{$index+1}->{meth}++;
98 	    }
99 	    else{
00 	      $mbias_2{CHH}->{$index+1}->{meth}++;
01 	    }
02 	  } 
03 	  elsif ($methylation_calls[$index] eq 'h') {
04 	    $counting{total_unmethylated_CHH_count}++;
05 	    print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
06 	    if ($read_identity == 1){
07 	      $mbias_1{CHH}->{$index+1}->{un}++;
08 	    }
09 	    else{
10 	      $mbias_2{CHH}->{$index+1}->{un}++;
11 	    }
12 	  }
13 	  elsif ($methylation_calls[$index] eq '.') {}
14 	  elsif (lc$methylation_calls[$index] eq 'u'){}
15 	  else{
16 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
17 	  }
18 	}
19       } elsif ($strand eq '-') {
20 	for my $index (0..$#methylation_calls) {
21 
22 	  if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
23 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
24 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
25 	    $cigar_offset += $cigar_mod;
26 	    $pos_offset += $pos_mod;
27 	  }
28 
29 	  ### Returning as soon as the methylation calls start overlapping
30 	  if ($start-$index+$pos_offset <= $end_read_1) {
31 	    return;
32 	  }
33 	
34 	  if ($methylation_calls[$index] eq 'X') {
35 	    $counting{total_meCHG_count}++;
36 	    print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
37 	    if ($read_identity == 1){
38 	      $mbias_1{CHG}->{$index+1}->{meth}++;
39 	    }
40 	    else{
41 	      $mbias_2{CHG}->{$index+1}->{meth}++;
42 	    }
43 	  }
44 	  elsif ($methylation_calls[$index] eq 'x') {
45 	    $counting{total_unmethylated_CHG_count}++;
46 	    print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
47 	    if ($read_identity == 1){
48 	      $mbias_1{CHG}->{$index+1}->{un}++;
49 	    }
50 	    else{
51 	      $mbias_2{CHG}->{$index+1}->{un}++;
52 	    }
53 	  }
54 	  elsif ($methylation_calls[$index] eq 'Z') {
55 	    $counting{total_meCpG_count}++;
56 	    print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
57 	    if ($read_identity == 1){
58 	      $mbias_1{CpG}->{$index+1}->{meth}++;
59 	    }
60 	    else{
61 	      $mbias_2{CpG}->{$index+1}->{meth}++;
62 	    }
63 	  }
64 	  elsif ($methylation_calls[$index] eq 'z') {
65 	    $counting{total_unmethylated_CpG_count}++;
66 	    print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
67 	    if ($read_identity == 1){
68 	      $mbias_1{CpG}->{$index+1}->{un}++;
69 	    }
70 	    else{
71 	      $mbias_2{CpG}->{$index+1}->{un}++;
72 	    }
73 	  }
74 	  elsif ($methylation_calls[$index] eq 'H') {
75 	    $counting{total_meCHH_count}++;
76 	    print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
77 	    if ($read_identity == 1){
78 	      $mbias_1{CHH}->{$index+1}->{meth}++;
79 	    }
80 	    else{
81 	      $mbias_2{CHH}->{$index+1}->{meth}++;
82 	    }
83 	  }
84 	  elsif ($methylation_calls[$index] eq 'h') {
85 	    $counting{total_unmethylated_CHH_count}++;
86 	    print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
87 	    if ($read_identity == 1){
88 	      $mbias_1{CHH}->{$index+1}->{un}++;
89 	    }
90 	    else{
91 	      $mbias_2{CHH}->{$index+1}->{un}++;
92 	    }
93 	  }
94 	  elsif ($methylation_calls[$index] eq '.') {}
95 	  elsif (lc$methylation_calls[$index] eq 'u'){}
96 	  else{
97 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
98 	  }
99 	}
00       } else {
01 	die "The strand orientation as neither + nor -: '$strand'\n";
02       }
03     }
04 
05     ### strand-specific methylation output
06     else {
07       if ($strand eq '+') {
08 	for my $index (0..$#methylation_calls) {
09 
10 	  if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
11 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
12 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
13 	    $cigar_offset += $cigar_mod;
14 	    $pos_offset += $pos_mod;
15 	  }
16 
17 	  ### Returning as soon as the methylation calls start overlapping
18 	  if ($start+$index+$pos_offset >= $end_read_1) {
19 	    return;
20 	  }
21 	
22 	  if ($methylation_calls[$index] eq 'X') {
23 	    $counting{total_meCHG_count}++;
24 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
25 	    if ($read_identity == 1){
26 	      $mbias_1{CHG}->{$index+1}->{meth}++;
27 	    }
28 	    else{
29 	      $mbias_2{CHG}->{$index+1}->{meth}++;
30 	    }
31 	  }
32 	  elsif ($methylation_calls[$index] eq 'x') {
33 	    $counting{total_unmethylated_CHG_count}++;
34 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
35 	    if ($read_identity == 1){
36 	      $mbias_1{CHG}->{$index+1}->{un}++;
37 	    }
38 	    else{
39 	      $mbias_2{CHG}->{$index+1}->{un}++;
40 	    }
41 	  }
42 	  elsif ($methylation_calls[$index] eq 'Z') {
43 	    $counting{total_meCpG_count}++;
44 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
45 	    if ($read_identity == 1){
46 	      $mbias_1{CpG}->{$index+1}->{meth}++;
47 	    }
48 	    else{
49 	      $mbias_2{CpG}->{$index+1}->{meth}++;
50 	    }
51 	  }
52 	  elsif ($methylation_calls[$index] eq 'z') {
53 	    $counting{total_unmethylated_CpG_count}++;
54 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
55 	    if ($read_identity == 1){
56 	      $mbias_1{CpG}->{$index+1}->{un}++;
57 	    }
58 	    else{
59 	      $mbias_2{CpG}->{$index+1}->{un}++;
60 	    }
61 	  }
62 	  elsif ($methylation_calls[$index] eq 'H') {
63 	    $counting{total_meCHH_count}++;
64 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
65 	    if ($read_identity == 1){
66 	      $mbias_1{CHH}->{$index+1}->{meth}++;
67 	    }
68 	    else{
69 	      $mbias_2{CHH}->{$index+1}->{meth}++;
70 	    }
71 	  }
72 	  elsif ($methylation_calls[$index] eq 'h') {
73 	    $counting{total_unmethylated_CHH_count}++;
74 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
75 	    if ($read_identity == 1){
76 	      $mbias_1{CHH}->{$index+1}->{un}++;
77 	    }
78 	    else{
79 	      $mbias_2{CHH}->{$index+1}->{un}++;
80 	    }
81 	  }
82 	  elsif ($methylation_calls[$index] eq '.') {}
83 	  elsif (lc$methylation_calls[$index] eq 'u'){}
84 	  else{
85 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
86 	  }	
87 	}
88       } elsif ($strand eq '-') {
89 	for my $index (0..$#methylation_calls) {
90 
91 	  if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
92 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
93 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
94 	    $cigar_offset += $cigar_mod;
95 	    $pos_offset += $pos_mod;
96 	  }
97 
98 	  ### Returning as soon as the methylation calls start overlapping
99 	  if ($start-$index+$pos_offset <= $end_read_1) {
00 	    return;
01 	  }
02 	
03 	  if ($methylation_calls[$index] eq 'X') {
04 	    $counting{total_meCHG_count}++;
05 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
06 	    if ($read_identity == 1){
07 	      $mbias_1{CHG}->{$index+1}->{meth}++;
08 	    }
09 	    else{
10 	      $mbias_2{CHG}->{$index+1}->{meth}++;
11 	    }
12 	  }
13 	  elsif ($methylation_calls[$index] eq 'x') {
14 	    $counting{total_unmethylated_CHG_count}++;
15 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
16 	    if ($read_identity == 1){
17 	      $mbias_1{CHG}->{$index+1}->{un}++;
18 	    }
19 	    else{
20 	      $mbias_2{CHG}->{$index+1}->{un}++;
21 	    }
22 	  }
23 	  elsif ($methylation_calls[$index] eq 'Z') {
24 	    $counting{total_meCpG_count}++;
25 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
26 	    if ($read_identity == 1){
27 	      $mbias_1{CpG}->{$index+1}->{meth}++;
28 	    }
29 	    else{
30 	      $mbias_2{CpG}->{$index+1}->{meth}++;
31 	    }
32 	  }
33 	  elsif ($methylation_calls[$index] eq 'z') {
34 	    $counting{total_unmethylated_CpG_count}++;
35 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
36 	    if ($read_identity == 1){
37 	      $mbias_1{CpG}->{$index+1}->{un}++;
38 	    }
39 	    else{
40 	      $mbias_2{CpG}->{$index+1}->{un}++;
41 	    }
42 	  }
43 	  elsif ($methylation_calls[$index] eq 'H') {
44 	    $counting{total_meCHH_count}++;
45 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
46 	    if ($read_identity == 1){
47 	      $mbias_1{CHH}->{$index+1}->{meth}++;
48 	    }
49 	    else{
50 	      $mbias_2{CHH}->{$index+1}->{meth}++;
51 	    }
52 	  }
53 	  elsif ($methylation_calls[$index] eq 'h') {
54 	    $counting{total_unmethylated_CHH_count}++;
55 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
56 	    if ($read_identity == 1){
57 	      $mbias_1{CHH}->{$index+1}->{un}++;
58 	    }
59 	    else{
60 	      $mbias_2{CHH}->{$index+1}->{un}++;
61 	    }
62 	  }
63 	  elsif ($methylation_calls[$index] eq '.') {}
64 	  elsif (lc$methylation_calls[$index] eq 'u'){}
65 	  else{
66 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
67 	  }
68 	}
69       } else {
70 	die "The strand orientation as neither + nor -: '$strand'\n";
71       }
72     }
73   }
74 
75   ### this is the default paired-end procedure allowing overlaps and using every single C position
76   else {
77     ### single-file CpG, CHG and CHH context output
78     if ($full) {
79       if ($strand eq '+') {
80 	for my $index (0..$#methylation_calls) {
81 
82 	  if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
83 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
84 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
85 	    $cigar_offset += $cigar_mod;
86 	    $pos_offset += $pos_mod;
87 	  }
88 
89 	  if ($methylation_calls[$index] eq 'X') {
90 	    $counting{total_meCHG_count}++;
91 	    print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
92 	    if ($read_identity == 1){
93 	      $mbias_1{CHG}->{$index+1}->{meth}++;
94 	    }
95 	    else{
96 	      $mbias_2{CHG}->{$index+1}->{meth}++;
97 	    }
98 	  }
99 	  elsif ($methylation_calls[$index] eq 'x') {
00 	    $counting{total_unmethylated_CHG_count}++;
01 	    print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
02 	    if ($read_identity == 1){
03 	      $mbias_1{CHG}->{$index+1}->{un}++;
04 	    }
05 	    else{
06 	      $mbias_2{CHG}->{$index+1}->{un}++;
07 	    }
08 	  }
09 	  elsif ($methylation_calls[$index] eq 'Z') {
10 	    $counting{total_meCpG_count}++;
11 	    print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
12 	    if ($read_identity == 1){
13 	      $mbias_1{CpG}->{$index+1}->{meth}++;
14 	    }
15 	    else{
16 	      $mbias_2{CpG}->{$index+1}->{meth}++;
17 	    }
18 	  }
19 	  elsif ($methylation_calls[$index] eq 'z') {
20 	    $counting{total_unmethylated_CpG_count}++;
21 	    print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
22 	    if ($read_identity == 1){
23 	      $mbias_1{CpG}->{$index+1}->{un}++;
24 	    }
25 	    else{
26 	      $mbias_2{CpG}->{$index+1}->{un}++;
27 	    }
28 	  }
29 	  elsif ($methylation_calls[$index] eq 'H') {
30 	    $counting{total_meCHH_count}++;
31 	    print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
32 	    if ($read_identity == 1){
33 	      $mbias_1{CHH}->{$index+1}->{meth}++;
34 	    }
35 	    else{
36 	      $mbias_2{CHH}->{$index+1}->{meth}++;
37 	    }
38 	  }
39 	  elsif ($methylation_calls[$index] eq 'h') {
40 	    $counting{total_unmethylated_CHH_count}++;
41 	    print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
42 	    if ($read_identity == 1){
43 	      $mbias_1{CHH}->{$index+1}->{un}++;
44 	    }
45 	    else{
46 	      $mbias_2{CHH}->{$index+1}->{un}++;
47 	    }
48 	  }
49 	  elsif ($methylation_calls[$index] eq '.') {}
50 	  elsif (lc$methylation_calls[$index] eq 'u'){}
51 	  else{
52 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
53 	  }
54 	}
55       } elsif ($strand eq '-') {
56 	for my $index (0..$#methylation_calls) {
57 
58 	  if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
59 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
60 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
61 	    $cigar_offset += $cigar_mod;
62 	    $pos_offset += $pos_mod;
63 	  }
64 
65 	  if ($methylation_calls[$index] eq 'X') {
66 	    $counting{total_meCHG_count}++;
67 	    print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
68 	    if ($read_identity == 1){
69 	      $mbias_1{CHG}->{$index+1}->{meth}++;
70 	    }
71 	    else{
72 	      $mbias_2{CHG}->{$index+1}->{meth}++;
73 	    }
74 	  }
75 	  elsif ($methylation_calls[$index] eq 'x') {
76 	    $counting{total_unmethylated_CHG_count}++;
77 	    print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
78 	    if ($read_identity == 1){
79 	      $mbias_1{CHG}->{$index+1}->{un}++;
80 	    }
81 	    else{
82 	      $mbias_2{CHG}->{$index+1}->{un}++;
83 	    }
84 	  }
85 	  elsif ($methylation_calls[$index] eq 'Z') {
86 	    $counting{total_meCpG_count}++;
87 	    print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
88 	    if ($read_identity == 1){
89 	      $mbias_1{CpG}->{$index+1}->{meth}++;
90 	    }
91 	    else{
92 	      $mbias_2{CpG}->{$index+1}->{meth}++;
93 	    }
94 	  }
95 	  elsif ($methylation_calls[$index] eq 'z') {
96 	    $counting{total_unmethylated_CpG_count}++;
97 	    print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
98 	    if ($read_identity == 1){
99 	      $mbias_1{CpG}->{$index+1}->{un}++;
00 	    }
01 	    else{
02 	      $mbias_2{CpG}->{$index+1}->{un}++;
03 	    }
04 	  }
05 	  elsif ($methylation_calls[$index] eq 'H') {
06 	    $counting{total_meCHH_count}++;
07 	    print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
08 	    if ($read_identity == 1){
09 	      $mbias_1{CHH}->{$index+1}->{meth}++;
10 	    }
11 	    else{
12 	      $mbias_2{CHH}->{$index+1}->{meth}++;
13 	    }
14 	  }
15 	  elsif ($methylation_calls[$index] eq 'h') {
16 	    $counting{total_unmethylated_CHH_count}++;
17 	    print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
18 	    if ($read_identity == 1){
19 	      $mbias_1{CHH}->{$index+1}->{un}++;
20 	    }
21 	    else{
22 	      $mbias_2{CHH}->{$index+1}->{un}++;
23 	    }
24 	  }
25 	  elsif ($methylation_calls[$index] eq '.') {}
26 	  elsif (lc$methylation_calls[$index] eq 'u'){}
27 	  else{
28 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
29 	  }
30 	}
31       } else {
32 	die "The strand orientation as neither + nor -: '$strand'\n";
33       }
34     }
35 
36     ### strand-specific methylation output
37     else {
38       if ($strand eq '+') {
39 	for my $index (0..$#methylation_calls) {
40 
41 	  if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
42 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
43 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
44 	    $cigar_offset += $cigar_mod;
45 	    $pos_offset += $pos_mod;
46 	  }
47 
48 	  if ($methylation_calls[$index] eq 'X') {
49 	    $counting{total_meCHG_count}++;
50 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
51 	    if ($read_identity == 1){
52 	      $mbias_1{CHG}->{$index+1}->{meth}++;
53 	    }
54 	    else{
55 	      $mbias_2{CHG}->{$index+1}->{meth}++;
56 	    }
57 	  }
58 	  elsif ($methylation_calls[$index] eq 'x') {
59 	    $counting{total_unmethylated_CHG_count}++;
60 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
61 	    if ($read_identity == 1){
62 	      $mbias_1{CHG}->{$index+1}->{un}++;
63 	    }
64 	    else{
65 	      $mbias_2{CHG}->{$index+1}->{un}++;
66 	    }
67 	  }
68 	  elsif ($methylation_calls[$index] eq 'Z') {
69 	    $counting{total_meCpG_count}++;
70 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
71 	    if ($read_identity == 1){
72 	      $mbias_1{CpG}->{$index+1}->{meth}++;
73 	    }
74 	    else{
75 	      $mbias_2{CpG}->{$index+1}->{meth}++;
76 	    }
77 	  }
78 	  elsif ($methylation_calls[$index] eq 'z') {
79 	    $counting{total_unmethylated_CpG_count}++;
80 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
81 	    if ($read_identity == 1){
82 	      $mbias_1{CpG}->{$index+1}->{un}++;
83 	    }
84 	    else{
85 	      $mbias_2{CpG}->{$index+1}->{un}++;
86 	    }
87 	  }
88 	  elsif ($methylation_calls[$index] eq 'H') {
89 	    $counting{total_meCHH_count}++;
90 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
91 	    if ($read_identity == 1){
92 	      $mbias_1{CHH}->{$index+1}->{meth}++;
93 	    }
94 	    else{
95 	      $mbias_2{CHH}->{$index+1}->{meth}++;
96 	    }
97 	  }
98 	  elsif ($methylation_calls[$index] eq 'h') {
99 	    $counting{total_unmethylated_CHH_count}++;
00 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
01 	    if ($read_identity == 1){
02 	      $mbias_1{CHH}->{$index+1}->{un}++;
03 	    }
04 	    else{
05 	      $mbias_2{CHH}->{$index+1}->{un}++;
06 	    }
07 	  }
08 	  elsif ($methylation_calls[$index] eq '.') {}
09 	  elsif (lc$methylation_calls[$index] eq 'u'){}
10 	  else{
11 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
12 	  }
13 	}
14       } elsif ($strand eq '-') {
15 	for my $index (0..$#methylation_calls) {
16 
17 	  if ($cigar and @comp_cigar){ # only needed for SAM reads with InDels
18 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
19 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
20 	    $cigar_offset += $cigar_mod;
21 	    $pos_offset += $pos_mod;
22 	  }
23 
24 	  if ($methylation_calls[$index] eq 'X') {
25 	    $counting{total_meCHG_count}++;
26 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
27 	    if ($read_identity == 1){
28 	      $mbias_1{CHG}->{$index+1}->{meth}++;
29 	    }
30 	    else{
31 	      $mbias_2{CHG}->{$index+1}->{meth}++;
32 	    }
33 	  }
34 	  elsif ($methylation_calls[$index] eq 'x') {
35 	    $counting{total_unmethylated_CHG_count}++;
36 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
37 	    if ($read_identity == 1){
38 	      $mbias_1{CHG}->{$index+1}->{un}++;
39 	    }
40 	    else{
41 	      $mbias_2{CHG}->{$index+1}->{un}++;
42 	    }
43 	  }
44 	  elsif ($methylation_calls[$index] eq 'Z') {
45 	    $counting{total_meCpG_count}++;
46 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
47 	    if ($read_identity == 1){
48 	      $mbias_1{CpG}->{$index+1}->{meth}++;
49 	    }
50 	    else{
51 	      $mbias_2{CpG}->{$index+1}->{meth}++;
52 	    }
53 	  }
54 	  elsif ($methylation_calls[$index] eq 'z') {
55 	    $counting{total_unmethylated_CpG_count}++;
56 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
57 	    if ($read_identity == 1){
58 	      $mbias_1{CpG}->{$index+1}->{un}++;
59 	    }
60 	    else{
61 	      $mbias_2{CpG}->{$index+1}->{un}++;
62 	    }
63 	  }
64 	  elsif ($methylation_calls[$index] eq 'H') {
65 	    $counting{total_meCHH_count}++;
66 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
67 	    if ($read_identity == 1){
68 	      $mbias_1{CHH}->{$index+1}->{meth}++;
69 	    }
70 	    else{
71 	      $mbias_2{CHH}->{$index+1}->{meth}++;
72 	    }
73 	  }
74 	  elsif ($methylation_calls[$index] eq 'h') {
75 	    $counting{total_unmethylated_CHH_count}++;
76 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
77 	    if ($read_identity == 1){
78 	      $mbias_1{CHH}->{$index+1}->{un}++;
79 	    }
80 	    else{
81 	      $mbias_2{CHH}->{$index+1}->{un}++;
82 	    }
83 	  }
84 	  elsif ($methylation_calls[$index] eq '.') {}
85 	  elsif (lc$methylation_calls[$index] eq 'u'){}
86 	  else{
87 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
88 	  }
89 	}
90       } else {
91 	die "The strand orientation as neither + nor -: '$strand'\n";
92       }
93     }
94   }
95 }
96 
97 sub check_cigar_string {
98   my ($index,$cigar_offset,$pos_offset,$strand,$comp_cigar) = @_;
99   # print "$index\t$cigar_offset\t$pos_offset\t$strand\t";
00   my ($new_cigar_offset,$new_pos_offset) = (0,0);
01 
02   if ($strand eq '+') {
03     #  print "### $strand strand @$comp_cigar[$index + $cigar_offset]\t";
04 
05     if (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'M'){ # sequence position matches the genomic position
06       #  warn "position needs no adjustment\n";
07     }
08 
09     elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'I'){ # insertion in the read sequence
10       $new_pos_offset -= 1; # we need to subtract the length of inserted bases from the genomic position
11       # warn "adjusted genomic position by -1 bp (insertion)\n";
12     }
13 
14     elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'D'){ # deletion in the read sequence
15       $new_cigar_offset += 1; # the composite cigar string does no longer match the methylation call index
16       $new_pos_offset += 1; # we need to add the length of deleted bases to get the genomic position
17       # warn "adjusted genomic position by +1 bp (deletion). Now looping through the CIGAR string until we hit another M or I\n";
18 
19       while ( ($index + $cigar_offset + $new_cigar_offset)  < (scalar @$comp_cigar) ){	
20 	if (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'M'){ # sequence position matches the genomic position
21 	  #  warn "position needs no adjustment\n";
22 	  last;
23 	}
24 	elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'I'){
25 	  $new_pos_offset -= 1; # we need to subtract the length of inserted bases from the genomic position
26 	  # warn "adjusted genomic position by another -1 bp (insertion)\n";
27 	  last;
28 	}
29 	elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'D'){ # deletion in the read sequence
30 	  $new_cigar_offset += 1; # the composite cigar string does no longer match the methylation call index
31 	  $new_pos_offset += 1; # we need to add the length of deleted bases to get the genomic position
32 	  # warn "adjusted genomic position by another +1 bp (deletion)\n";
33 	}
34 	else{
35 	  die "The CIGAR string contained undefined operations in addition to 'M', 'I' and 'D': '@$comp_cigar[$index + $cigar_offset + $new_cigar_offset]'\n";
36 	}
37       }
38     }
39     else{
40       die "The CIGAR string contained undefined operations in addition to 'M', 'I' and 'D': '@$comp_cigar[$index + $cigar_offset + $new_cigar_offset]'\n";
41     }
42   }
43 
44   elsif ($strand eq '-') {
45     # print "### $strand strand @$comp_cigar[$index + $cigar_offset]\t";
46 
47     if (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'M'){ # sequence position matches the genomic position
48      # warn "position needs no adjustment\n";
49     }
50 
51     elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'I'){ # insertion in the read sequence
52       $new_pos_offset += 1; # we need to add the length of inserted bases to the genomic position
53       # warn "adjusted genomic position by +1 bp (insertion)\n";
54     }
55 
56     elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'D'){ # deletion in the read sequence
57       $new_cigar_offset += 1; # the composite cigar string does no longer match the methylation call index
58       $new_pos_offset -= 1; # we need to subtract the length of deleted bases to get the genomic position
59       # warn "adjusted genomic position by -1 bp (deletion). Now looping through the CIGAR string until we hit another M or I\n";
60 
61       while ( ($index + $cigar_offset + $new_cigar_offset)  < (scalar @$comp_cigar) ){	
62 	if (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'M'){ # sequence position matches the genomic position
63 	  # warn "Found new 'M' operation; position needs no adjustment\n";
64 	  last;
65 	}
66 	elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'I'){
67 	  $new_pos_offset += 1; # we need to subtract the length of inserted bases from the genomic position
68 	  # warn "Found new 'I' operation; adjusted genomic position by another +1 bp (insertion)\n";
69 	  last;
70 	}
71 	elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'D'){ # deletion in the read sequence
72 	  $new_cigar_offset += 1; # the composite cigar string does no longer match the methylation call index
73 	  $new_pos_offset -= 1; # we need to subtract the length of deleted bases to get the genomic position
74 	  # warn "adjusted genomic position by another -1 bp (deletion)\n";
75 	}
76 	else{
77 	  die "The CIGAR string contained undefined operations in addition to 'M', 'I' and 'D': '@$comp_cigar[$index + $cigar_offset + $new_cigar_offset]'\n";
78 	}
79       }
80     }
81     else{
82       die "The CIGAR string contained undefined operations in addition to 'M', 'I' and 'D': '@$comp_cigar[$index + $cigar_offset + $new_cigar_offset]'\n";
83     }
84   }
85   # print "new cigar offset: $new_cigar_offset\tnew pos offset: $new_pos_offset\n";
86   return ($new_cigar_offset,$new_pos_offset);
87 }
88 
89 sub print_individual_C_methylation_states_single_end{
90 
91   my ($meth_call,$chrom,$start,$id,$strand,$filehandle_index,$cigar) = @_;
92   my @methylation_calls = split(//,$meth_call);
93 
94   #################################################################
95   ### . for bases not involving cytosines                       ###
96   ### X for methylated C in CHG context (was protected)         ###
97   ### x for not methylated C in CHG context (was converted)     ###
98   ### H for methylated C in CHH context (was protected)         ###
99   ### h for not methylated C in CHH context (was converted)     ###
00   ### Z for methylated C in CpG context (was protected)         ###
01   ### z for not methylated C in CpG context (was converted)     ###
02   #################################################################
03 
04   my $methyl_CHG_count = 0;
05   my $methyl_CHH_count = 0;
06   my $methyl_CpG_count = 0;
07   my $unmethylated_CHG_count = 0;
08   my $unmethylated_CHH_count = 0;
09   my $unmethylated_CpG_count = 0;
10 
11   my $pos_offset = 0; # this is only relevant for SAM reads with insertions or deletions
12   my $cigar_offset = 0; # again, this is only relevant for SAM reads containing indels
13 
14   my @comp_cigar;
15 
16   if ($cigar){ # parsing CIGAR string
17 
18     ### Checking whether the CIGAR string is a linear genomic match or whether if requires indel processing
19     if ($cigar =~ /^\d+M$/){
20       #  warn "See!? I told you so! $cigar\n";
21       # sleep(1);
22     }
23     else{
24 
25       my @len;
26       my @ops;
27 
28       @len = split (/\D+/,$cigar); # storing the length per operation
29       @ops = split (/\d+/,$cigar); # storing the operation
30       shift @ops; # remove the empty first element
31       # die "CIGAR string contained a non-matching number of lengths and operations: id: $id\nmeth call: $meth_call\nCIGAR: $cigar\n".join(" ",@len)."\n".join(" ",@ops)."\n" unless (scalar @len == scalar @ops);
32       die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len == scalar @ops);
33 
34       foreach my $index (0..$#len){
35 	foreach (1..$len[$index]){
36 	  # print  "$ops[$index]";
37 	  push @comp_cigar, $ops[$index];
38 	}
39       }
40     }
41     # warn "\nDetected CIGAR string: $cigar\n";
42     # warn "Length of methylation call: ",length $meth_call,"\n";
43     # warn "number of operations: ",scalar @ops,"\n";
44     # warn "number of length digits: ",scalar @len,"\n\n";
45     # print @comp_cigar,"\n";
46     # print "$meth_call\n\n";
47     # sleep (1);
48   }
49 
50   ### adjusting the start position for all reads mapping to the reverse strand
51   if ($strand eq '-') {
52 
53     if (@comp_cigar){ # only needed for SAM reads with InDels
54       @comp_cigar  = reverse@comp_cigar; # the CIGAR string needs to be reversed for all reads aligning to the reverse strand, too
55       # print @comp_cigar,"\n";
56     }
57 
58     unless ($ignore){  ### if --ignore was specified the start position has already been corrected
59 
60       if ($cigar){ ### SAM format
61 	if ($cigar =~ /^(\d+)M$/){ # linear match
62 	  $start += $1 - 1;
63 	}
64 	else{ # InDel read
65 	  my $MD_count = 0;
66 	  foreach (@comp_cigar){
67 	    ++$MD_count if ($_ eq 'M' or $_ eq 'D'); # Matching bases or deletions affect the genomic position of the 3' ends of reads, insertions don't
68 	  }
69 	  $start += $MD_count - 1;
70 	}
71       }
72       else{ ### vanilla format
73 	$start += length($meth_call)-1;
74       }
75     }
76   }
77 
78   ### THIS IS THE CpG and Non-CpG SECTION (OPTIONAL)
79 
80   ### single-file CpG and other-context output
81   if ($full and $merge_non_CpG) {
82     if ($strand eq '+') {
83       for my $index (0..$#methylation_calls) {
84 	
85 	if ($cigar and @comp_cigar){ # only needed for SAM alignments with InDels
86 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
87 	  # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
88 	  $cigar_offset += $cigar_mod;
89 	  $pos_offset += $pos_mod;
90 	}
91 	
92 	### methylated Cs (any context) will receive a forward (+) orientation
93 	### not methylated Cs (any context) will receive a reverse (-) orientation
94 	if ($methylation_calls[$index] eq 'X') {
95 	  $counting{total_meCHG_count}++;
96 	  print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
97 	  $mbias_1{CHG}->{$index+1}->{meth}++;
98 	}
99 	elsif ($methylation_calls[$index] eq 'x') {
00 	  $counting{total_unmethylated_CHG_count}++;
01 	  print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
02 	  $mbias_1{CHG}->{$index+1}->{un}++;
03 	}
04 	elsif ($methylation_calls[$index] eq 'Z') {
05 	  $counting{total_meCpG_count}++;
06 	  print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
07 	  $mbias_1{CpG}->{$index+1}->{meth}++;
08 	}
09 	elsif ($methylation_calls[$index] eq 'z') {
10 	  $counting{total_unmethylated_CpG_count}++;
11 	  print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
12 	  $mbias_1{CpG}->{$index+1}->{un}++;
13 	}
14 	elsif ($methylation_calls[$index] eq 'H') {
15 	  $counting{total_meCHH_count}++;
16 	  print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
17 	  $mbias_1{CHH}->{$index+1}->{meth}++;
18 	}
19 	elsif ($methylation_calls[$index] eq 'h') {
20 	  $counting{total_unmethylated_CHH_count}++;
21 	  print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
22 	  $mbias_1{CHH}->{$index+1}->{un}++;
23 	}
24 	elsif ($methylation_calls[$index] eq '.') {}
25 	elsif (lc$methylation_calls[$index] eq 'u'){}
26 	else{
27 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
28 	}
29       }
30     }
31     elsif ($strand eq '-') {
32 
33       for my $index (0..$#methylation_calls) {
34 	### methylated Cs (any context) will receive a forward (+) orientation
35 	### not methylated Cs (any context) will receive a reverse (-) orientation
36 
37 	if ($cigar and @comp_cigar){ # only needed for SAM entries with InDels
38 	  # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
39 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
40 	  $cigar_offset += $cigar_mod;
41 	  $pos_offset += $pos_mod;
42 	}
43 
44 	if ($methylation_calls[$index] eq 'X') {
45 	  $counting{total_meCHG_count}++;
46 	  print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
47 	  $mbias_1{CHG}->{$index+1}->{meth}++;
48 	}
49 	elsif ($methylation_calls[$index] eq 'x') {
50 	  $counting{total_unmethylated_CHG_count}++;
51 	  print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
52 	  $mbias_1{CHG}->{$index+1}->{un}++;
53 	}
54 	elsif ($methylation_calls[$index] eq 'Z') {
55 	  $counting{total_meCpG_count}++;
56 	  print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
57 	  $mbias_1{CpG}->{$index+1}->{meth}++;
58 	}
59 	elsif ($methylation_calls[$index] eq 'z') {
60 	  $counting{total_unmethylated_CpG_count}++;
61 	  print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
62 	  $mbias_1{CpG}->{$index+1}->{un}++;
63 	}
64 	elsif ($methylation_calls[$index] eq 'H') {
65 	  $counting{total_meCHH_count}++;
66 	  print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
67 	  $mbias_1{CHH}->{$index+1}->{meth}++;
68 	}
69 	elsif ($methylation_calls[$index] eq 'h') {
70 	  $counting{total_unmethylated_CHH_count}++;
71 	  print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
72 	  $mbias_1{CHH}->{$index+1}->{un}++;
73 	}
74 	elsif ($methylation_calls[$index] eq '.'){}
75 	elsif (lc$methylation_calls[$index] eq 'u'){}
76 	else{
77 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
78 	}
79       }
80     }
81     else {
82       die "The strand information was neither + nor -: $strand\n";
83     }
84   }
85 
86   ### strand-specific methylation output
87   elsif ($merge_non_CpG) {
88     if ($strand eq '+') {
89       for my $index (0..$#methylation_calls) {
90 	### methylated Cs (any context) will receive a forward (+) orientation
91 	### not methylated Cs (any context) will receive a reverse (-) orientation
92 
93 	if ($cigar and @comp_cigar){ # only needed for SAM reads with Indels
94 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
95 	  $cigar_offset += $cigar_mod;
96 	  $pos_offset += $pos_mod;
97 	}
98 
99 	if ($methylation_calls[$index] eq 'X') {
00 	  $counting{total_meCHG_count}++;
01 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
02 	  $mbias_1{CHG}->{$index+1}->{meth}++;
03 	}
04 	elsif ($methylation_calls[$index] eq 'x') {
05 	  $counting{total_unmethylated_CHG_count}++;
06 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
07 	  $mbias_1{CHG}->{$index+1}->{un}++;
08 	}
09 	elsif ($methylation_calls[$index] eq 'Z') {
10 	  $counting{total_meCpG_count}++;
11 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
12 	  $mbias_1{CpG}->{$index+1}->{meth}++;
13 	}
14 	elsif ($methylation_calls[$index] eq 'z') {
15 	  $counting{total_unmethylated_CpG_count}++;
16 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
17 	  $mbias_1{CpG}->{$index+1}->{un}++;
18 	}
19 	elsif ($methylation_calls[$index] eq 'H') {
20 	  $counting{total_meCHH_count}++;
21 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
22 	  $mbias_1{CHH}->{$index+1}->{meth}++;
23 	}
24 	elsif ($methylation_calls[$index] eq 'h') {
25 	  $counting{total_unmethylated_CHH_count}++;
26 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
27 	  $mbias_1{CHH}->{$index+1}->{un}++;
28 	}
29 	elsif ($methylation_calls[$index] eq '.') {}
30 	elsif (lc$methylation_calls[$index] eq 'u'){}
31 	else{
32 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
33 	}
34       }
35     }
36     elsif ($strand eq '-') {
37 
38       for my $index (0..$#methylation_calls) {
39 	### methylated Cs (any context) will receive a forward (+) orientation
40 	### not methylated Cs (any context) will receive a reverse (-) orientation
41 
42     	if ($cigar and @comp_cigar){ # only needed for SAM reads with Indels
43 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
44 	  $cigar_offset += $cigar_mod;
45 	  $pos_offset += $pos_mod;
46 	}
47 
48 	if ($methylation_calls[$index] eq 'X') {
49 	  $counting{total_meCHG_count}++;
50 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
51 	  $mbias_1{CHG}->{$index+1}->{meth}++;
52 	}
53 	elsif ($methylation_calls[$index] eq 'x') {
54 	  $counting{total_unmethylated_CHG_count}++;
55 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
56 	  $mbias_1{CHG}->{$index+1}->{un}++;
57 	}
58 	elsif ($methylation_calls[$index] eq 'Z') {
59 	  $counting{total_meCpG_count}++;
60 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
61 	  $mbias_1{CpG}->{$index+1}->{meth}++;
62 	}
63 	elsif ($methylation_calls[$index] eq 'z') {
64 	  $counting{total_unmethylated_CpG_count}++;
65 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
66 	  $mbias_1{CpG}->{$index+1}->{un}++;
67 	}
68 	elsif ($methylation_calls[$index] eq 'H') {
69 	  $counting{total_meCHH_count}++;
70 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
71 	  $mbias_1{CHH}->{$index+1}->{meth}++;
72 	}
73 	elsif ($methylation_calls[$index] eq 'h') {
74 	  $counting{total_unmethylated_CHH_count}++;
75 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
76 	  $mbias_1{CHH}->{$index+1}->{un}++;
77 	}
78 	elsif ($methylation_calls[$index] eq '.') {}
79 	elsif (lc$methylation_calls[$index] eq 'u'){}
80 	else{
81 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
82 	}
83       }
84     }
85     else {
86       die "The strand information was neither + nor -: $strand\n";
87     }
88   }
89 
90   ### THIS IS THE 3-CONTEXT (CpG, CHG and CHH) DEFAULT SECTION
91 
92   elsif ($full) {
93     if ($strand eq '+') {
94       for my $index (0..$#methylation_calls) {
95 	### methylated Cs (any context) will receive a forward (+) orientation
96 	### not methylated Cs (any context) will receive a reverse (-) orientation
97 
98 	if ($cigar and @comp_cigar){ # only needed for SAM reads with Indels
99 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
00 	  $cigar_offset += $cigar_mod;
01 	  $pos_offset += $pos_mod;
02 	}
03 	
04 	if ($methylation_calls[$index] eq 'X') {
05 	  $counting{total_meCHG_count}++;
06 	  print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
07 	  $mbias_1{CHG}->{$index+1}->{meth}++;
08 	}
09 	elsif ($methylation_calls[$index] eq 'x') {
10 	  $counting{total_unmethylated_CHG_count}++;
11 	  print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
12 	  $mbias_1{CHG}->{$index+1}->{un}++;
13 	}
14 	elsif ($methylation_calls[$index] eq 'Z') {
15 	  $counting{total_meCpG_count}++;
16 	  print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
17 	  $mbias_1{CpG}->{$index+1}->{meth}++;
18 	}
19 	elsif ($methylation_calls[$index] eq 'z') {
20 	  $counting{total_unmethylated_CpG_count}++;
21 	  print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
22 	  $mbias_1{CpG}->{$index+1}->{un}++;
23 	}
24 	elsif ($methylation_calls[$index] eq 'H') {
25 	  $counting{total_meCHH_count}++;
26 	  print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
27 	  $mbias_1{CHH}->{$index+1}->{meth}++;
28 	}
29 	elsif ($methylation_calls[$index] eq 'h') {
30 	  $counting{total_unmethylated_CHH_count}++;
31 	  print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
32 	  $mbias_1{CHH}->{$index+1}->{un}++;
33 	}
34 	elsif ($methylation_calls[$index] eq '.') {}
35 	elsif (lc$methylation_calls[$index] eq 'u'){}
36 	else{
37 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n" unless($mbias_only);
38 	}
39       }
40     }
41     elsif ($strand eq '-') {
42 
43       for my $index (0..$#methylation_calls) {
44 	### methylated Cs (any context) will receive a forward (+) orientation
45 	### not methylated Cs (any context) will receive a reverse (-) orientation
46 
47 	if ($cigar and @comp_cigar){ # only needed for SAM reads with Indels
48 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
49 	  $cigar_offset += $cigar_mod;
50 	  $pos_offset += $pos_mod;
51 	}
52 	
53 	if ($methylation_calls[$index] eq 'X') {
54 	  $counting{total_meCHG_count}++;
55 	  print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
56 	  $mbias_1{CHG}->{$index+1}->{meth}++;
57 	}
58 	elsif ($methylation_calls[$index] eq 'x') {
59 	  $counting{total_unmethylated_CHG_count}++;
60 	  print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
61 	  $mbias_1{CHG}->{$index+1}->{un}++;
62 	}
63 	elsif ($methylation_calls[$index] eq 'Z') {
64 	  $counting{total_meCpG_count}++;
65 	  print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
66 	  $mbias_1{CpG}->{$index+1}->{meth}++;
67 	}
68 	elsif ($methylation_calls[$index] eq 'z') {
69 	  $counting{total_unmethylated_CpG_count}++;
70 	  print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
71 	  $mbias_1{CpG}->{$index+1}->{un}++;
72 	}
73 	elsif ($methylation_calls[$index] eq 'H') {
74 	  $counting{total_meCHH_count}++;
75 	  print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
76 	  $mbias_1{CHH}->{$index+1}->{meth}++;
77 	}
78 	elsif ($methylation_calls[$index] eq 'h') {
79 	  $counting{total_unmethylated_CHH_count}++;
80 	  print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
81 	  $mbias_1{CHH}->{$index+1}->{un}++;
82 	}
83 	elsif ($methylation_calls[$index] eq '.') {}
84 	elsif (lc$methylation_calls[$index] eq 'u'){}
85 	else{
86 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
87 	}
88       }
89     }
90     else {
91       die "The read had a strand orientation which was neither + nor -: $strand\n";
92     }
93   }
94 
95   ### strand-specific methylation output
96   else {
97     if ($strand eq '+') {
98       for my $index (0..$#methylation_calls) {
99 	### methylated Cs (any context) will receive a forward (+) orientation
00 	### not methylated Cs (any context) will receive a reverse (-) orientation
01 
02 	if ($cigar and @comp_cigar){ # only needed for SAM reads with Indels
03 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
04 	  $cigar_offset += $cigar_mod;
05 	  $pos_offset += $pos_mod;
06 	}
07 
08 	if ($methylation_calls[$index] eq 'X') {
09 	  $counting{total_meCHG_count}++;
10 	  print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
11 	  $mbias_1{CHG}->{$index+1}->{meth}++;
12 	}
13 	elsif ($methylation_calls[$index] eq 'x') {
14 	  $counting{total_unmethylated_CHG_count}++;
15 	  print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
16 	  $mbias_1{CHG}->{$index+1}->{un}++;
17 	}
18 	elsif ($methylation_calls[$index] eq 'Z') {
19 	  $counting{total_meCpG_count}++;
20 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
21 	  $mbias_1{CpG}->{$index+1}->{meth}++;
22 	}
23 	elsif ($methylation_calls[$index] eq 'z') {
24 	  $counting{total_unmethylated_CpG_count}++;
25 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
26 	  $mbias_1{CpG}->{$index+1}->{un}++;
27 	}
28 	elsif ($methylation_calls[$index] eq 'H') {
29 	  $counting{total_meCHH_count}++;
30 	  print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
31 	  $mbias_1{CHH}->{$index+1}->{meth}++;
32 	}
33 	elsif ($methylation_calls[$index] eq 'h') {
34 	  $counting{total_unmethylated_CHH_count}++;
35 	  print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
36 	  $mbias_1{CHH}->{$index+1}->{un}++;
37 	}
38 	elsif ($methylation_calls[$index] eq '.') {}
39 	elsif (lc$methylation_calls[$index] eq 'u'){}
40 	else{
41 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
42 	}
43       }
44     }
45     elsif ($strand eq '-') {
46 
47       for my $index (0..$#methylation_calls) {
48 	### methylated Cs (any context) will receive a forward (+) orientation
49 	### not methylated Cs (any context) will receive a reverse (-) orientation
50 
51 	if ($cigar and @comp_cigar){ # only needed for SAM reads with Indels
52 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
53 	  $cigar_offset += $cigar_mod;
54 	  $pos_offset += $pos_mod;
55 	}
56 
57 	if ($methylation_calls[$index] eq 'X') {
58 	  $counting{total_meCHG_count}++;
59 	  print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
60 	  $mbias_1{CHG}->{$index+1}->{meth}++;
61 	}
62 	elsif ($methylation_calls[$index] eq 'x') {
63 	  $counting{total_unmethylated_CHG_count}++;
64 	  print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
65 	  $mbias_1{CHG}->{$index+1}->{un}++;
66 	}
67 	elsif ($methylation_calls[$index] eq 'Z') {
68 	  $counting{total_meCpG_count}++;
69 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
70 	  $mbias_1{CpG}->{$index+1}->{meth}++;
71 	}
72 	elsif ($methylation_calls[$index] eq 'z') {
73 	  $counting{total_unmethylated_CpG_count}++;
74 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
75 	  $mbias_1{CpG}->{$index+1}->{un}++;
76 	}
77 	elsif ($methylation_calls[$index] eq 'H') {
78 	  $counting{total_meCHH_count}++;
79 	  print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
80 	  $mbias_1{CHH}->{$index+1}->{meth}++;
81 	}
82 	elsif ($methylation_calls[$index] eq 'h') {
83 	  $counting{total_unmethylated_CHH_count}++;
84 	  print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n" unless($mbias_only);
85 	  $mbias_1{CHH}->{$index+1}->{un}++;
86 	}
87 	elsif ($methylation_calls[$index] eq '.') {}
88 	elsif (lc$methylation_calls[$index] eq 'u'){}
89 	else{
90 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
91 	}
92       }
93     }
94     else {
95       die "The strand information was neither + nor -: $strand\n";
96     }
97   }
98 }
99 
00 
01 
02 sub print_helpfile{
03 
04  print << 'HOW_TO';
05 
06 
07 DESCRIPTION
08 
09 The following is a brief description of all options to control the Bismark
10 methylation extractor. The script reads in a bisulfite read alignment results file 
11 produced by the Bismark bisulfite mapper and extracts the methylation information
12 for individual cytosines. This information is found in the methylation call field
13 which can contain the following characters:
14 
15        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
16        ~~~   X   for methylated C in CHG context                      ~~~
17        ~~~   x   for not methylated C CHG                             ~~~
18        ~~~   H   for methylated C in CHH context                      ~~~
19        ~~~   h   for not methylated C in CHH context                  ~~~
20        ~~~   Z   for methylated C in CpG context                      ~~~
21        ~~~   z   for not methylated C in CpG context                  ~~~
22        ~~~   U   for methylated C in Unknown context (CN or CHN       ~~~
23        ~~~   u   for not methylated C in Unknown context (CN or CHN)  ~~~
24        ~~~   .   for any bases not involving cytosines                ~~~
25        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
26 
27 The methylation extractor outputs result files for cytosines in CpG, CHG and CHH
28 context (this distinction is actually already made in Bismark itself). As the methylation
29 information for every C analysed can produce files which easily have tens or even hundreds of
30 millions of lines, file sizes can become very large and more difficult to handle. The C
31 methylation info additionally splits cytosine methylation calls up into one of the four possible
32 strands a given bisulfite read aligned against:
33 
34              OT      original top strand
35              CTOT    complementary to original top strand
36 
37              OB      original bottom strand
38              CTOB    complementary to original bottom strand
39 
40 Thus, by default twelve individual output files are being generated per input file (unless
41 --comprehensive is specified, see below). The output files can be imported into a genome
42 viewer, such as SeqMonk, and re-combined into a single data group if desired (in fact
43 unless the bisulfite reads were generated preserving directionality it doesn't make any
44 sense to look at the data in a strand-specific manner). Strand-specific output files can
45 optionally be skipped, in which case only three output files for CpG, CHG or CHH context
46 will be generated. For both the strand-specific and comprehensive outputs there is also
47 the option to merge both non-CpG contexts (CHG and CHH) into one single non-CpG context.
48 
49 
50 The output files are in the following format (tab delimited):
51 
52 <sequence_id>     <strand>      <chromosome>     <position>     <methylation call>
53 
54 
55 USAGE: methylation_extractor [options] <filenames>
56 
57 
58 ARGUMENTS:
59 ==========
60 
61 <filenames>              A space-separated list of Bismark result files in SAM format from
62                          which methylation information is extracted for every cytosine in
63                          the reads. For alignment files in the older custom Bismark output
64                          see option '--vanilla'.
65 
66 OPTIONS:
67 
68 -s/--single-end          Input file(s) are Bismark result file(s) generated from single-end
69                          read data. Specifying either --single-end or --paired-end is
70                          mandatory.
71 
72 -p/--paired-end          Input file(s) are Bismark result file(s) generated from paired-end
73                          read data. Specifying either --paired-end or --single-end is
74                          mandatory.
75 
76 --vanilla                The Bismark result input file(s) are in the old custom Bismark format
77                          (up to version 0.5.x) and not in SAM format which is the default as
78                          of Bismark version 0.6.x or higher. Default: OFF.
79 
80 --no_overlap             For paired-end reads it is theoretically possible that read_1 and
81                          read_2 overlap. This option avoids scoring overlapping methylation
82                          calls twice (only methylation calls of read 1 are used for in the process
83                          since read 1 has historically higher quality basecalls than read 2).
84                          Whilst this option removes a bias towards more methylation calls
85                          in the center of sequenced fragments it may de facto remove a sizable
86                          proportion of the data. This option is highly recommended for paired-end
87                          data.
88 
89 --ignore <int>           Ignore the first <int> bp from the 5' end of Read 1 when processing the
90                          methylation call string. This can remove e.g. a restriction enzyme site
91                          at the start of each read or any other source of bias (e.g. PBAT-Seq data).
92 
93 --ignore_r2 <int>        Ignore the first <int> bp from the 5' end of Read 2 of paired-end sequencing
94                          results only. Since the first couple of bases in Read 2 of BS-Seq experiments
95                          show a severe bias towards non-methylation as a result of end-repairing
96                          sonicated fragments with unmethylated cytosines (see M-bias plot), it is
97                          recommended that the first couple of bp of Read 2 are removed before
98                          starting downstream analysis. Please see the section on M-bias plots in the
99                          Bismark User Guide for more details.
00 
01 --comprehensive          Specifying this option will merge all four possible strand-specific
02                          methylation info into context-dependent output files. The default
03 
04                          contexts are:
05                           - CpG context
06                           - CHG context
07                           - CHH context
08 
09 --merge_non_CpG          This will produce two output files (in --comprehensive mode) or eight
10                          strand-specific output files (default) for Cs in
11                           - CpG context
12                           - non-CpG context
13 
14 --report                 Prints out a short methylation summary as well as the paramaters used to run
15                          this script.
16 
17 --no_header              Suppresses the Bismark version header line in all output files for more convenient
18                          batch processing.
19 
20 -o/--output DIR          Allows specification of a different output directory (absolute or relative
21                          path). If not specified explicitely, the output will be written to the current directory.
22 
23 --samtools_path          The path to your Samtools installation, e.g. /home/user/samtools/. Does not need to be specified
24                          explicitly if Samtools is in the PATH already.
25 
26 --gzip                   The methylation extractor files (CpG_OT_..., CpG_OB_... etc) will be written out in
27                          a GZIP compressed form to save disk space. This option does not work on bedGraph and
28                          genome-wide cytosine reports as they are 'tiny' anyway.
29 
30 --version                Displays version information.
31 
32 -h/--help                Displays this help file and exits.
33 
34 --mbias_only             The methylation extractor will read the entire file but only output the M-bias table and plots as 
35                          well as a report (optional) and then quit. Default: OFF.
36 
37 
38 
39 bedGraph specific options:
40 ==========================
41 
42 --bedGraph               After finishing the methylation extraction, the methylation output is written into a
43                          sorted bedGraph file that reports the position of a given cytosine and its methylation 
44                          state (in %, see details below). The methylation extractor output is temporarily split up into
45                          temporary files, one per chromosome (written into the current directory or folder
46                          specified with -o/--output); these temp files are then used for sorting and deleted
47                          afterwards. By default, only cytosines in CpG context will be sorted. The option
48                          '--CX_context' may be used to report all cytosines irrespective of sequence context
49                          (this will take MUCH longer!). The default folder for temporary files during the sorting
50                          process is the output directory. The bedGraph conversion step is performed by the external
51                          module 'bismark2bedGraph'; this script needs to reside in the same folder as the 
52                          bismark_methylation_extractor itself.
53 
54 
55 --cutoff [threshold]     The minimum number of times a methylation state has to be seen for that nucleotide
56                          before its methylation percentage is reported. Default: 1.
57 
58 --remove_spaces          Replaces whitespaces in the sequence ID field with underscores to allow sorting.
59 
60 
61 --CX/--CX_context        The sorted bedGraph output file contains information on every single cytosine that was covered
62                          in the experiment irrespective of its sequence context. This applies to both forward and
63                          reverse strands. Please be aware that this option may generate large temporary and output files
64                          and may take a long time to sort (up to many hours). Default: OFF.
65                          (i.e. Default = CpG context only).
66 
67 --buffer_size <string>   This allows you to specify the main memory sort buffer when sorting the methylation information.
68                          Either specify a percentage of physical memory by appending % (e.g. --buffer_size 50%) or
69 			 a multiple of 1024 bytes, e.g. 'K' multiplies by 1024, 'M' by 1048576 and so on for 'T' etc. 
70                          (e.g. --buffer_size 20G). For more information on sort type 'info sort' on a command line.
71                          Defaults to 2G.
72 
73 --scaffolds/--gazillion  Users working with unfinished genomes sporting tens or even hundreds of thousands of
74                          scaffolds/contigs/chromosomes frequently encountered errors with pre-sorting reads to
75                          individual chromosome files. These errors were caused by the operating system's limit
76                          of the number of filehandle that can be written to at any one time (typically 1024; to
77                          find out this limit on Linux, type: ulimit -a).
78                          To bypass the limitation of open filehandles, the option --scaffolds does not pre-sort
79                          methylation calls into individual chromosome files. Instead, all input files are
80                          temporarily merged into a single file (unless there is only a single file), and this
81                          file will then be sorted by both chromosome AND position using the Unix sort command.
82                          Please be aware that this option might take a looooong time to complete, depending on
83                          the size of the input files, and the memory you allocate to this process (see --buffer_size).
84                          Nevertheless, it seems to be working.
85 
86 --ample_memory           Using this option will not sort chromosomal positions using the UNIX 'sort' command, but will
87                          instead use two arrays to sort methylated and unmethylated calls. This may result in a faster
88                          sorting process of very large files, but this comes at the cost of a larger memory footprint
89                          (two arrays of the length of the largest human chromosome 1 (~250M bp) consume around 16GB
90                          of RAM). Due to overheads in creating and looping through these arrays it seems that it will
91                          actually be *slower* for small files (few million alignments), and we are currently testing at
92                          which point it is advisable to use this option. Note that --ample_memory is not compatible
93                          with options '--scaffolds/--gazillion' (as it requires pre-sorted files to begin with).
94 
95 
96 
97 Genome-wide cytosine methylation report specific options:
98 =========================================================
99 
00 --cytosine_report        After the conversion to bedGraph has completed, the option '--cytosine_report' produces a
01                          genome-wide methylation report for all cytosines in the genome. By default, the output uses 1-based
02                          chromosome coordinates (zero-based cords are optional) and reports CpG context only (all
03                          cytosine context is optional). The output considers all Cs on both forward and reverse strands and
04                          reports their position, strand, trinucleotide content and methylation state (counts are 0 if not
05                          covered). The cytsoine report conversion step is performed by the external module 
06                          'bedGraph2cytosine'; this script needs to reside in the same folder as the bismark_methylation_extractor
07                          itself.
08 
09 --CX/--CX_context        The output file contains information on every single cytosine in the genome irrespective of
10                          its context. This applies to both forward and reverse strands. Please be aware that this will
11                          generate output files with > 1.1 billion lines for a mammalian genome such as human or mouse.
12                          Default: OFF (i.e. Default = CpG context only).
13 
14 --zero_based             Uses zero-based coordinates like used in e.g. bed files instead of 1-based coordinates. Default: OFF.
15 
16 --genome_folder <path>   Enter the genome folder you wish to use to extract sequences from (full path only). Accepted
17                          formats are FastA files ending with '.fa' or '.fasta'. Specifying a genome folder path is mandatory.
18 
19 --split_by_chromosome    Writes the output into individual files for each chromosome instead of a single output file. Files
20                          will be named to include the input filename and the chromosome number.
21 
22 
23 
24 OUTPUT:
25 
26 The bismark_methylation_extractor output is in the form:
27 ========================================================
28 <seq-ID>  <methylation state*>  <chromosome>  <start position (= end position)>  <methylation call>
29 
30 * Methylated cytosines receive a '+' orientation,
31 * Unmethylated cytosines receive a '-' orientation.
32 
33 
34 
35 The bedGraph output (optional) looks like this (tab-delimited; 0-based start coords, 1-based end coords):
36 =========================================================================================================
37 
38 track type=bedGraph (header line)
39 
40 <chromosome>  <start position>  <end position>  <methylation percentage>
41 
42 
43 
44 The coverage output looks like this (tab-delimited, 1-based genomic coords):
45 ============================================================================
46 
47 <chromosome>  <start position>  <end position>  <methylation percentage>  <count methylated>  <count non-methylated>
48 
49 
50 
51 The genome-wide cytosine methylation output file is tab-delimited in the following format:
52 ==========================================================================================
53 <chromosome>  <position>  <strand>  <count methylated>  <count non-methylated>  <C-context>  <trinucleotide context>
54 
55 
56 
57 This script was last modified on 25 November 2013.
58 
59 HOW_TO
60 }