annotate bismark_methylation_extractor @ 0:36d124f44c0a draft

inital commit

1 #!/usr/bin/perl
2 use warnings;
3 use strict;
4 $|++;
5 use Getopt::Long;
6 use Cwd;
7 use Carp;
8 
9 my @filenames; # input files
10 my %counting;
11 my $parent_dir = getcwd();
12 
13 my %fhs;
14 
15 my $version = 'v0.7.7';
16 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) = process_commandline();
17 
18 
19 ### only needed for bedGraph output
20 my @sorting_files; # if files are to be written to bedGraph format, these are the methylation extractor output files
21 my @methylcalls = qw (0 0 0); # [0] = methylated, [1] = unmethylated, [2] = total
22 my @bedfiles;
23 
24 ### only needed for genome-wide cytosine methylation report
25 my %chromosomes;
26 
27 ##############################################################################################
28 ### Summarising Run Parameters
29 ##############################################################################################
30 
31 ### METHYLATION EXTRACTOR
32 
33 warn "Summarising Bismark methylation extractor parameters:\n";
34 warn '='x63,"\n";
35 
36 if ($single){
37   if ($vanilla){
38     warn "Bismark single-end vanilla format specified\n";
39   }
40   else{
41     warn "Bismark single-end SAM format specified (default)\n"; # default
42   }
43 }
44 elsif ($paired){
45   if ($vanilla){
46     warn "Bismark paired-end vanilla format specified\n";
47   }
48   else{
49     warn "Bismark paired-end SAM format specified (default)\n"; # default
50   }
51 }
52 
53 if ($ignore){
54   warn "First $ignore bases will be disregarded when processing the methylation call string\n";
55 }
56 
57 if ($full){
58   warn "Strand-specific outputs will be skipped. Separate output files for cytosines in CpG, CHG and CHH context will be generated\n";
59 }
60 if ($merge_non_CpG){
61   warn "Merge CHG and CHH context to non-CpG context specified\n";
62 }
63 ### output directory
64 if ($output_dir eq ''){
65   warn "Output will be written to the current directory ('$parent_dir')\n";
66 }
67 else{
68   warn "Output path specified as: $output_dir\n";
69 }
70 
71 
72 sleep (1);
73 
74 ### BEDGRAPH
75 
76 if ($bedGraph){
77   warn "\n\nSummarising bedGraph parameters:\n";
78   warn '='x63,"\n";
79 
80   if ($counts){
81     warn "Generating additional output in bedGraph format including methylating counts (output format: <Chromosome> <Start Position> <End Position> <Methylation Percentage> <count methylated> <count non-methylated>)\n";
82   }
83   else{
84     warn "Generating additional sorted output in bedGraph format (output format: <Chromosome> <Start Position> <End Position> <Methylation Percentage>)\n";
85   }
86 
87   warn "Using a cutoff of $coverage_threshold read(s) to report cytosine positions\n";
88 
89   if ($CX_context){
90     warn "Reporting and sorting methylation information for all cytosine context (sorting may take a long time, you have been warned ...)\n";
91   }
92   else{ # default
93     $CpG_only = 1;
94     warn "Reporting and sorting cytosine methylation information in CpG context only (default)\n";
95   }
96 
97   if ($remove){
98     warn "White spaces in read ID names will be removed prior to sorting\n";
99   }
00 
01   sleep (1);
02 
03   if ($cytosine_report){
04     warn "\n\nSummarising genome-wide cytosine methylation report parameters:\n";
05     warn '='x63,"\n";
06     warn "Generating comprehensive genome-wide cytosine report (output format: <Chromosome> <Start Position> <End Position> <Methylation Percentage> )\n";
07 
08 
09     if ($CX_context){
10       warn "Reporting methylation for all cytosine contexts. Be aware that this will generate enormous files\n";
11     }
12     else{ # default
13       $CpG_only = 1;
14       warn "Reporting cytosine methylation in CpG context only (default)\n";
15     }
16 
17     if ($split_by_chromosome){
18       warn "Splitting the cytosine report output up into individual files for each chromosome\n";
19     }
20 
21     ### Zero-based coordinates
22     if ($zero){
23       warn "Using zero-based genomic coordinates (user-defined)\n";
24     }
25     else{ # default, 1-based coords
26       warn "Using 1-based genomic coordinates (default)\n";
27     }
28 
29     ### GENOME folder
30     if ($genome_folder){
31       unless ($genome_folder =~/\/$/){
32 	$genome_folder =~ s/$/\//;
33       }
34       warn "Genome folder was specified as $genome_folder\n";
35     }
36     else{
37       $genome_folder  = '/data/public/Genomes/Mouse/NCBIM37/';
38       warn "Using the default genome folder /data/public/Genomes/Mouse/NCBIM37/\n";
39     }
40     sleep (1);
41   }
42 }
43 
44 warn "\n";
45 sleep (5);
46 
47 ######################################################
48 ### PROCESSING FILES
49 ######################################################
50 
51 foreach my $filename (@filenames){
52   # resetting counters and filehandles
53   %fhs = ();
54   %counting =(
55 	      total_meCHG_count => 0,
56 	      total_meCHH_count => 0,
57 	      total_meCpG_count => 0,
58 	      total_unmethylated_CHG_count => 0,
59 	      total_unmethylated_CHH_count => 0,
60 	      total_unmethylated_CpG_count => 0,
61 	      sequences_count => 0,
62 	     );
63   @sorting_files = ();
64   @bedfiles = ();
65 
66   process_Bismark_results_file($filename);
67 
68   if ($bedGraph){
69     my $out = $filename;
70     $out =~ s/sam$//;
71     $out =~ s/txt$//;
72     $out =~ s/$/bedGraph/;
73 
74     my $bedGraph_output = $out;
75     open (OUT,'>',$output_dir.$out) or die $!;
76     # warn "Writing bedGraph to file: $out\n";
77 
78     process_bedGraph_output();
79     close OUT or die $!;
80 
81     ### genome-wide cytosine methylation report requires bedGraph processing anyway
82     if ($cytosine_report){
83       my $cytosine_out = $out;
84       $cytosine_out =~ s/bedGraph$//;
85 
86       read_genome_into_memory();
87       generate_genome_wide_cytosine_report($bedGraph_output,$cytosine_out);
88     }
89   }
90 }
91 
92 
93 sub process_commandline{
94   my $help;
95   my $single_end;
96   my $paired_end;
97   my $ignore;
98   my $genomic_fasta;
99   my $full;
00   my $report;
01   my $extractor_version;
02   my $no_overlap;
03   my $merge_non_CpG;
04   my $vanilla;
05   my $output_dir;
06   my $no_header;
07   my $bedGraph;
08   my $coverage_threshold = 1; # Minimum number of reads covering before calling methylation status
09   my $remove;
10   my $counts;
11   my $cytosine_report;
12   my $genome_folder;
13   my $zero;
14   my $CpG_only;
15   my $CX_context;
16   my $split_by_chromosome;
17 
18 
19   my $command_line = GetOptions ('help|man' => \$help,
20 				 'p|paired-end' => \$paired_end,
21 				 's|single-end' => \$single_end,
22 				 'fasta' => \$genomic_fasta,
23 				 'ignore=i' => \$ignore,
24 				 'comprehensive' => \$full,
25 				 'report' => \$report,
26 				 'version' => \$extractor_version,
27 				 'no_overlap' => \$no_overlap,
28 				 'merge_non_CpG' => \$merge_non_CpG, 
29 				 'vanilla' => \$vanilla,
30 				 'o|output=s' => \$output_dir,
31 				 'no_header' => \$no_header,
32 				 'bedGraph' => \$bedGraph,
33 				 "cutoff=i" => \$coverage_threshold,
34 				 "remove_spaces" => \$remove,
35 				 "counts" => \$counts,
36 				 "cytosine_report" => \$cytosine_report,
37 				 'g|genome_folder=s' => \$genome_folder,
38 				 "zero_based" => \$zero,	
39 				 "CX|CX_context" => \$CX_context,
40 				 "split_by_chromosome" => \$split_by_chromosome,
41 				);
42 
43   ### EXIT ON ERROR if there were errors with any of the supplied options
44   unless ($command_line){
45     die "Please respecify command line options\n";
46   }
47 
48   ### HELPFILE
49   if ($help){
50     print_helpfile();
51     exit;
52   }
53 
54   if ($extractor_version){
55     print << "VERSION";
56 
57 
58                            Bismark Methylation Extractor
59 
60    Bismark Extractor Version: $version Copyright 2010-12 Felix Krueger, Babraham Bioinformatics
61                 www.bioinformatics.babraham.ac.uk/projects/bismark/
62 
63 
64 VERSION
65     exit;
66   }
67 
68 
69   ### no files provided
70   unless (@ARGV){
71     die "You need to provide one or more Bismark files to create an individual C methylation output. Please respecify!\n";
72   }
73   @filenames = @ARGV;
74 
75   warn "\n *** Bismark methylation extractor version $version ***\n\n";
76 
77   ### IGNORING <INT> bases at the start of the read when processing the methylation call string
78   unless ($ignore){
79     $ignore = 0;
80   }
81   ### PRINT A REPORT
82   unless ($report){
83     $report = 0;
84   }
85 
86   ### OUTPUT DIR PATH
87   if ($output_dir){
88     unless ($output_dir =~ /\/$/){
89       $output_dir =~ s/$/\//;
90     }
91   }
92   else{
93     $output_dir = '';
94   }
95 
96   ### NO HEADER
97   unless ($no_header){
98     $no_header = 0;
99   }
00 
01   ### OLD (VANILLA) OUTPUT FORMAT
02   unless ($vanilla){
03     $vanilla = 0;
04   }
05 
06   if ($single_end){
07     $paired_end = 0;   ### SINGLE END ALIGNMENTS
08   }
09   elsif ($paired_end){
10     $single_end = 0;   ### PAIRED-END ALIGNMENTS
11   }
12   else{
13     die "Please specify whether the supplied file(s) are in Bismark single-end or paired-end format\n\n";
14   }
15 
16   ### NO OVERLAP
17   if ($no_overlap){
18     die "The option '--no_overlap' can only be specified for paired-end input!\n" unless ($paired_end);
19   }
20   else{
21     $no_overlap = 0;
22   }
23 
24   ### COMPREHENSIVE OUTPUT
25   unless ($full){
26     $full = 0;
27   }
28 
29   ### MERGE NON-CpG context
30   unless ($merge_non_CpG){
31     $merge_non_CpG = 0;
32   }
33 
34   ### remove white spaces in read ID (needed for sorting using the sort command
35   unless ($remove){
36     $remove = 0;
37   }
38 
39   ### COVERAGE THRESHOLD FOR gedGraph OUTPUT
40   unless (defined $coverage_threshold){
41     unless ($coverage_threshold > 0){
42       die "Please select a coverage greater than 0 (positive integers only)\n";
43     }
44     $coverage_threshold = 1;
45   }
46 
47   if ($zero){
48     die "Option '--zero' is only available if  '--cytosine_report' is specified as well. Please respecify\n" unless ($cytosine_report);
49   }
50 
51   if ($CX_context){
52     die "Option '--CX_context' is only available if  '--cytosine_report' or '--bedGraph' is specified as well. Please respecify\n" unless ($cytosine_report or $bedGraph);
53   }
54 
55   if ($cytosine_report){
56 
57     ### GENOME folder
58     if ($genome_folder){
59       unless ($genome_folder =~/\/$/){
60 	$genome_folder =~ s/$/\//;
61       }
62     }
63     else{
64       die "Please specify a genome folder to proceed (full path only)\n";
65     }
66 
67     unless ($bedGraph){
68       warn "Setting the option '--bedGraph' since this is required for the genome-wide cytosine report\n";
69       $bedGraph = 1;
70     }
71     unless ($counts){
72       warn "Setting the option '--counts' since this is required for the genome-wide cytosine report\n";
73       $counts = 1;
74     }
75     warn "\n";
76   }
77 
78   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);
79 }
80 
81 
82 sub process_Bismark_results_file{
83   my $filename = shift;
84 
85   warn "\nNow reading in Bismark result file $filename\n\n";
86 
87   if ($filename =~ /\.gz$/) {
88     open (IN,"zcat $filename |") or die "Can't open gzipped file $filename: $!\n";
89   } else {
90     open (IN,$filename) or die "Can't open file $filename: $!\n";
91   }
92 
93   ### Vanilla and SAM output need to read different numbers of header lines
94   if ($vanilla) {
95     my $bismark_version = <IN>; ## discarding the Bismark version info
96     chomp $bismark_version;
97     $bismark_version =~ s/\r//; # replaces \r line feed
98     $bismark_version =~  s/Bismark version: //;
99     if ($bismark_version =~ /^\@/) {
00       warn "Detected \@ as the first character of the version information. Is it possible that the file is in SAM format?\n\n";
01       sleep (2);
02     }
03 
04     unless ($version eq $bismark_version){
05       die "The methylation extractor and Bismark itself need to be of the same version!\n\nVersions used:\nmethylation extractor: '$version'\nBismark: '$bismark_version'\n";
06     }
07   } else {
08     # If the read is in SAM format (default) it can either start with @ header lines or start with alignments directly.
09     # We are reading from it further down
10   }
11 
12   my $output_filename = (split (/\//,$filename))[-1];
13 
14   ### OPENING OUTPUT-FILEHANDLES
15   if ($report) {
16     my $report_filename = $output_filename;
17     $report_filename =~ s/[\.sam|\.txt]$//;
18     $report_filename =~ s/$/_splitting_report.txt/;
19     $report_filename = $output_dir . $report_filename;
20     open (REPORT,'>',$report_filename) or die "Failed to write to file $report_filename $!\n";
21   }
22 
23   if ($report) {
24     print REPORT "$output_filename\n\n";
25     print REPORT "Parameters used to extract methylation information:\n";
26     if ($paired) {
27       if ($vanilla) {
28 	print REPORT "Bismark result file: paired-end (vanilla Bismark format)\n";
29       } else {
30 	print REPORT "Bismark result file: paired-end (SAM format)\n"; # default
31       }
32     }
33 
34     if ($single) {
35       if ($vanilla) {
36 	print REPORT "Bismark result file: single-end (vanilla Bismark format)\n";
37       } else {
38 	print REPORT "Bismark result file: single-end (SAM format)\n"; # default
39       }
40     }
41 
42     if ($ignore) {
43       print REPORT "Ignoring first $ignore bases\n";
44     }
45 
46     if ($full) {
47       print REPORT "Output specified: comprehensive\n";
48     } else {
49       print REPORT "Output specified: strand-specific (default)\n";
50     }
51 
52     if ($no_overlap) {
53       print REPORT "No overlapping methylation calls specified\n";
54     }
55     if ($genomic_fasta) {
56       print REPORT "Genomic equivalent sequences will be printed out in FastA format\n";
57     }
58     if ($merge_non_CpG) {
59       print REPORT "Methylation in CHG and CHH context will be merged into \"non-CpG context\" output\n";
60     }
61 
62     print REPORT "\n";
63   }
64 
65   ### CpG-context and non-CpG context. THIS SECTION IS OPTIONAL
66   ### if --comprehensive AND --merge_non_CpG was specified we are only writing out one CpG-context and one Any-Other-context result file
67   if ($full and $merge_non_CpG) {
68     my $cpg_output = my $other_c_output = $output_filename;
69     ### C in CpG context
70     $cpg_output =~ s/^/CpG_context_/;
71     $cpg_output =~ s/sam$/txt/;
72     $cpg_output =~ s/$/.txt/ unless ($cpg_output =~ /\.txt$/);
73     $cpg_output = $output_dir . $cpg_output;
74     push @sorting_files,$cpg_output;
75     open ($fhs{CpG_context},'>',$cpg_output) or die "Failed to write to $cpg_output $! \n";
76     print "Writing result file containing methylation information for C in CpG context to $cpg_output\n";
77 
78     unless ($no_header) {
79       print {$fhs{CpG_context}} "Bismark methylation extractor version $version\n";
80     }
81 
82     ### C in any other context than CpG
83     $other_c_output =~ s/^/Non_CpG_context_/;
84     $other_c_output =~ s/sam$/txt/;
85     $other_c_output =~ s/$/.txt/ unless ($other_c_output =~ /\.txt$/);
86     $other_c_output = $output_dir . $other_c_output;
87     push @sorting_files,$other_c_output;
88     open ($fhs{other_context},'>',$other_c_output) or die "Failed to write to $other_c_output $!\n";
89     print "Writing result file containing methylation information for C in any other context to $other_c_output\n";
90 
91     unless ($no_header) {
92       print {$fhs{other_context}} "Bismark methylation extractor version $version\n";
93     }
94   }
95 
96   ### 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
97   elsif ($merge_non_CpG) {
98 
99     my $cpg_ot = my $cpg_ctot = my $cpg_ctob = my $cpg_ob = $output_filename;
00 
01     ### For cytosines in CpG context
02     $cpg_ot =~ s/^/CpG_OT_/;
03     $cpg_ot =~ s/sam$/txt/;
04     $cpg_ot =~ s/$/.txt/ unless ($cpg_ot =~ /\.txt$/);
05     $cpg_ot = $output_dir . $cpg_ot;
06     push @sorting_files,$cpg_ot;
07     open ($fhs{0}->{CpG},'>',$cpg_ot) or die "Failed to write to $cpg_ot $!\n";
08     print "Writing result file containing methylation information for C in CpG context from the original top strand to $cpg_ot\n";
09 
10     unless($no_header){
11       print {$fhs{0}->{CpG}} "Bismark methylation extractor version $version\n";
12     }
13 
14     $cpg_ctot =~ s/^/CpG_CTOT_/;
15     $cpg_ctot =~ s/sam$/txt/;
16     $cpg_ctot =~ s/$/.txt/ unless ($cpg_ctot =~ /\.txt$/);
17     $cpg_ctot = $output_dir . $cpg_ctot;
18     push @sorting_files,$cpg_ctot;
19     open ($fhs{1}->{CpG},'>',$cpg_ctot) or die "Failed to write to $cpg_ctot $!\n";
20     print "Writing result file containing methylation information for C in CpG context from the complementary to original top strand to $cpg_ctot\n";
21 
22     unless($no_header){
23       print {$fhs{1}->{CpG}} "Bismark methylation extractor version $version\n";
24     }
25 
26     $cpg_ctob =~ s/^/CpG_CTOB_/;
27     $cpg_ctob =~ s/sam$/txt/;
28     $cpg_ctob =~ s/$/.txt/ unless ($cpg_ctob =~ /\.txt$/);
29     $cpg_ctob = $output_dir . $cpg_ctob;
30     push @sorting_files,$cpg_ctob;
31     open ($fhs{2}->{CpG},'>',$cpg_ctob) or die "Failed to write to $cpg_ctob $!\n";
32     print "Writing result file containing methylation information for C in CpG context from the complementary to original bottom strand to $cpg_ctob\n";
33 
34     unless($no_header){
35       print {$fhs{2}->{CpG}}  "Bismark methylation extractor version $version\n";
36     }
37 
38     $cpg_ob =~ s/^/CpG_OB_/;
39     $cpg_ob =~ s/sam$/txt/;
40     $cpg_ob =~ s/$/.txt/ unless ($cpg_ob =~ /\.txt$/);
41     $cpg_ob = $output_dir . $cpg_ob;
42     push @sorting_files,$cpg_ob;
43     open ($fhs{3}->{CpG},'>',$cpg_ob) or die "Failed to write to $cpg_ob $!\n";
44     print "Writing result file containing methylation information for C in CpG context from the original bottom strand to $cpg_ob\n\n";
45 
46     unless($no_header){
47       print {$fhs{3}->{CpG}}  "Bismark methylation extractor version $version\n";
48     }
49 
50     ### For cytosines in Non-CpG (CC, CT or CA) context
51     my $other_c_ot = my $other_c_ctot = my $other_c_ctob = my $other_c_ob = $output_filename;
52 
53     $other_c_ot =~ s/^/Non_CpG_OT_/;
54     $other_c_ot =~ s/sam$/txt/;
55     $other_c_ot =~ s/$/.txt/ unless ($other_c_ot =~ /\.txt$/);
56     $other_c_ot = $output_dir . $other_c_ot;
57     push @sorting_files,$other_c_ot;
58     open ($fhs{0}->{other_c},'>',$other_c_ot) or die "Failed to write to $other_c_ot $!\n";
59     print "Writing result file containing methylation information for C in any other context from the original top strand to $other_c_ot\n";
60 
61     unless($no_header){
62       print {$fhs{0}->{other_c}} "Bismark methylation extractor version $version\n";
63     }
64 
65     $other_c_ctot =~ s/^/Non_CpG_CTOT_/;
66     $other_c_ctot =~ s/sam$/txt/;
67     $other_c_ctot =~ s/$/.txt/ unless ($other_c_ctot =~ /\.txt$/);
68     $other_c_ctot = $output_dir . $other_c_ctot;
69     push @sorting_files,$other_c_ctot;
70     open ($fhs{1}->{other_c},'>',$other_c_ctot) or die "Failed to write to $other_c_ctot $!\n";
71     print "Writing result file containing methylation information for C in any other context from the complementary to original top strand to $other_c_ctot\n";
72 
73     unless($no_header){
74       print {$fhs{1}->{other_c}} "Bismark methylation extractor version $version\n";
75     }
76 
77     $other_c_ctob =~ s/^/Non_CpG_CTOB_/;
78     $other_c_ctob =~ s/sam$/txt/;
79     $other_c_ctob =~ s/$/.txt/ unless ($other_c_ctob =~ /\.txt$/);
80     $other_c_ctob = $output_dir . $other_c_ctob;
81     push @sorting_files,$other_c_ctob;
82     open ($fhs{2}->{other_c},'>',$other_c_ctob) or die "Failed to write to $other_c_ctob $!\n";
83     print "Writing result file containing methylation information for C in any other context from the complementary to original bottom strand to $other_c_ctob\n";
84 
85     unless($no_header){
86       print {$fhs{2}->{other_c}} "Bismark methylation extractor version $version\n";
87     }
88 
89     $other_c_ob =~ s/^/Non_CpG_OB_/;
90     $other_c_ob =~ s/sam$/txt/;
91     $other_c_ob =~ s/$/.txt/ unless ($other_c_ob =~ /\.txt$/);
92     $other_c_ob = $output_dir . $other_c_ob;
93     push @sorting_files,$other_c_ob;
94     open ($fhs{3}->{other_c},'>',$other_c_ob) or die "Failed to write to $other_c_ob $!\n";
95     print "Writing result file containing methylation information for C in any other context from the original bottom strand to $other_c_ob\n\n";
96 
97     unless($no_header){
98       print {$fhs{3}->{other_c}} "Bismark methylation extractor version $version\n";
99     }
00   }
01   ### THIS SECTION IS THE DEFAULT (CpG, CHG and CHH context)
02 
03   ### if --comprehensive was specified we are only writing one file per context
04   elsif ($full) {
05     my $cpg_output = my $chg_output =  my $chh_output = $output_filename;
06     ### C in CpG context
07     $cpg_output =~ s/^/CpG_context_/;
08     $cpg_output =~ s/sam$/txt/;
09     $cpg_output =~ s/$/.txt/ unless ($cpg_output =~ /\.txt$/);
10     $cpg_output = $output_dir . $cpg_output;
11     push @sorting_files,$cpg_output;
12     open ($fhs{CpG_context},'>',$cpg_output) or die "Failed to write to $cpg_output $! \n";
13     print "Writing result file containing methylation information for C in CpG context to $cpg_output\n";
14 
15     unless($no_header){
16       print {$fhs{CpG_context}} "Bismark methylation extractor version $version\n";
17     }
18 
19     ### C in CHG context
20     $chg_output =~ s/^/CHG_context_/;
21     $chg_output =~ s/sam$/txt/;
22     $chg_output =~ s/$/.txt/ unless ($chg_output =~ /\.txt$/);
23     $chg_output = $output_dir . $chg_output;
24     push @sorting_files,$chg_output;
25     open ($fhs{CHG_context},'>',$chg_output) or die "Failed to write to $chg_output $!\n";
26     print "Writing result file containing methylation information for C in CHG context to $chg_output\n";
27 
28     unless($no_header){
29       print {$fhs{CHG_context}} "Bismark methylation extractor version $version\n";
30     }
31 
32     ### C in CHH context
33     $chh_output =~ s/^/CHH_context_/;
34     $chh_output =~ s/sam$/txt/;
35     $chh_output =~ s/$/.txt/ unless ($chh_output =~ /\.txt$/);
36     $chh_output = $output_dir . $chh_output;
37     push @sorting_files, $chh_output;
38     open ($fhs{CHH_context},'>',$chh_output) or die "Failed to write to $chh_output $!\n";
39     print "Writing result file containing methylation information for C in CHH context to $chh_output\n";
40 
41     unless($no_header){
42       print {$fhs{CHH_context}} "Bismark methylation extractor version $version\n";
43     }
44   }
45   ### else we will write out 12 different output files, depending on where the (first) unique best alignment was found
46   else {
47     my $cpg_ot = my $cpg_ctot = my $cpg_ctob = my $cpg_ob = $output_filename;
48 
49     ### For cytosines in CpG context
50     $cpg_ot =~ s/^/CpG_OT_/;
51     $cpg_ot =~ s/sam$/txt/;
52     $cpg_ot =~ s/$/.txt/ unless ($cpg_ot =~ /\.txt$/);
53     $cpg_ot = $output_dir . $cpg_ot;
54     push @sorting_files,$cpg_ot;
55     open ($fhs{0}->{CpG},'>',$cpg_ot) or die "Failed to write to $cpg_ot $!\n";
56     print "Writing result file containing methylation information for C in CpG context from the original top strand to $cpg_ot\n";
57 
58     unless($no_header){
59       print {$fhs{0}->{CpG}} "Bismark methylation extractor version $version\n";
60     }
61 
62     $cpg_ctot =~ s/^/CpG_CTOT_/;
63     $cpg_ctot =~ s/sam$/txt/;
64     $cpg_ctot =~ s/$/.txt/ unless ($cpg_ctot =~ /\.txt$/);
65     $cpg_ctot = $output_dir . $cpg_ctot;
66     push @sorting_files,$cpg_ctot;
67     open ($fhs{1}->{CpG},'>',$cpg_ctot) or die "Failed to write to $cpg_ctot $!\n";
68     print "Writing result file containing methylation information for C in CpG context from the complementary to original top strand to $cpg_ctot\n";
69 
70     unless($no_header){
71       print {$fhs{1}->{CpG}} "Bismark methylation extractor version $version\n";
72     }
73 
74     $cpg_ctob =~ s/^/CpG_CTOB_/;
75     $cpg_ctob =~ s/sam$/txt/;
76     $cpg_ctob =~ s/$/.txt/ unless ($cpg_ctob =~ /\.txt$/);
77     $cpg_ctob = $output_dir . $cpg_ctob;
78     push @sorting_files,$cpg_ctob;
79     open ($fhs{2}->{CpG},'>',$cpg_ctob) or die "Failed to write to $cpg_ctob $!\n";
80     print "Writing result file containing methylation information for C in CpG context from the complementary to original bottom strand to $cpg_ctob\n";
81 
82     unless($no_header){
83       print {$fhs{2}->{CpG}}  "Bismark methylation extractor version $version\n";
84     }
85 
86     $cpg_ob =~ s/^/CpG_OB_/;
87     $cpg_ob =~ s/sam$/txt/;
88     $cpg_ob =~ s/$/.txt/ unless ($cpg_ob =~ /\.txt$/);
89     $cpg_ob = $output_dir . $cpg_ob;
90     push @sorting_files,$cpg_ob;
91     open ($fhs{3}->{CpG},'>',$cpg_ob) or die "Failed to write to $cpg_ob $!\n";
92     print "Writing result file containing methylation information for C in CpG context from the original bottom strand to $cpg_ob\n\n";
93 
94     unless($no_header){
95       print {$fhs{3}->{CpG}}  "Bismark methylation extractor version $version\n";
96     }
97 
98     ### For cytosines in CHG context
99     my $chg_ot = my $chg_ctot = my $chg_ctob = my $chg_ob = $output_filename;
00 
01     $chg_ot =~ s/^/CHG_OT_/;
02     $chg_ot =~ s/sam$/txt/;
03     $chg_ot =~ s/$/.txt/ unless ($chg_ot =~ /\.txt$/);
04     $chg_ot = $output_dir . $chg_ot;
05     push @sorting_files,$chg_ot;
06     open ($fhs{0}->{CHG},'>',$chg_ot) or die "Failed to write to $chg_ot $!\n";
07     print "Writing result file containing methylation information for C in CHG context from the original top strand to $chg_ot\n";
08 
09     unless($no_header){
10       print {$fhs{0}->{CHG}} "Bismark methylation extractor version $version\n";
11     }
12 
13     $chg_ctot =~ s/^/CHG_CTOT_/;
14     $chg_ctot =~ s/sam$/txt/;
15     $chg_ctot =~ s/$/.txt/ unless ($chg_ctot =~ /\.txt$/);
16     $chg_ctot = $output_dir . $chg_ctot;
17     push @sorting_files,$chg_ctot;
18     open ($fhs{1}->{CHG},'>',$chg_ctot) or die "Failed to write to $chg_ctot $!\n";
19     print "Writing result file containing methylation information for C in CHG context from the complementary to original top strand to $chg_ctot\n";
20 
21     unless($no_header){
22       print {$fhs{1}->{CHG}} "Bismark methylation extractor version $version\n";
23     }
24 
25     $chg_ctob =~ s/^/CHG_CTOB_/;
26     $chg_ctob =~ s/sam$/txt/;
27     $chg_ctob =~ s/$/.txt/ unless ($chg_ctob =~ /\.txt$/);
28     $chg_ctob = $output_dir . $chg_ctob;
29     push @sorting_files,$chg_ctob;
30     open ($fhs{2}->{CHG},'>',$chg_ctob) or die "Failed to write to $chg_ctob $!\n";
31     print "Writing result file containing methylation information for C in CHG context from the complementary to original bottom strand to $chg_ctob\n";
32 
33     unless($no_header){
34       print {$fhs{2}->{CHG}} "Bismark methylation extractor version $version\n";
35     }
36 
37     $chg_ob =~ s/^/CHG_OB_/;
38     $chg_ob =~ s/sam$/txt/;
39     $chg_ob =~ s/$/.txt/ unless ($chg_ob =~ /\.txt$/);
40     $chg_ob = $output_dir . $chg_ob;
41     push @sorting_files,$chg_ob;
42     open ($fhs{3}->{CHG},'>',$chg_ob) or die "Failed to write to $chg_ob $!\n";
43     print "Writing result file containing methylation information for C in CHG context from the original bottom strand to $chg_ob\n\n";
44 
45     unless($no_header){
46       print {$fhs{3}->{CHG}} "Bismark methylation extractor version $version\n";
47     }
48 
49     ### For cytosines in CHH context
50     my $chh_ot = my $chh_ctot = my $chh_ctob = my $chh_ob = $output_filename;
51 
52     $chh_ot =~ s/^/CHH_OT_/;
53     $chh_ot =~ s/sam$/txt/;
54     $chh_ot =~ s/$/.txt/ unless ($chh_ot =~ /\.txt$/);
55     $chh_ot = $output_dir . $chh_ot;
56     push @sorting_files,$chh_ot;
57     open ($fhs{0}->{CHH},'>',$chh_ot) or die "Failed to write to $chh_ot $!\n";
58     print "Writing result file containing methylation information for C in CHH context from the original top strand to $chh_ot\n";
59 
60     unless($no_header){
61       print {$fhs{0}->{CHH}} "Bismark methylation extractor version $version\n";
62     }
63 
64     $chh_ctot =~ s/^/CHH_CTOT_/;
65     $chh_ctot =~ s/sam$/txt/;
66     $chh_ctot =~ s/$/.txt/ unless ($chh_ctot =~ /\.txt$/);
67     $chh_ctot = $output_dir . $chh_ctot;
68     push @sorting_files,$chh_ctot;
69     open ($fhs{1}->{CHH},'>',$chh_ctot) or die "Failed to write to $chh_ctot $!\n";
70     print "Writing result file containing methylation information for C in CHH context from the complementary to original top strand to $chh_ctot\n";
71 
72     unless($no_header){
73       print {$fhs{1}->{CHH}} "Bismark methylation extractor version $version\n";
74     }
75 
76     $chh_ctob =~ s/^/CHH_CTOB_/;
77     $chh_ctob =~ s/sam$/txt/;
78     $chh_ctob =~ s/$/.txt/ unless ($chh_ctob =~ /\.txt$/);
79     $chh_ctob = $output_dir . $chh_ctob;
80     push @sorting_files,$chh_ctob;
81     open ($fhs{2}->{CHH},'>',$chh_ctob) or die "Failed to write to $chh_ctob $!\n";
82     print "Writing result file containing methylation information for C in CHH context from the complementary to original bottom strand to $chh_ctob\n";
83 
84     unless($no_header){
85       print {$fhs{2}->{CHH}} "Bismark methylation extractor version $version\n";
86     }
87 
88     $chh_ob =~ s/^/CHH_OB_/;
89     $chh_ob =~ s/sam$/txt/;
90     $chh_ob =~ s/$/.txt/ unless ($chh_ob =~ /\.txt$/);
91     $chh_ob = $output_dir . $chh_ob;
92     push @sorting_files,$chh_ob;
93     open ($fhs{3}->{CHH},'>',$chh_ob) or die "Failed to write to $chh_ob $!\n";
94     print "Writing result file containing methylation information for C in CHH context from the original bottom strand to $chh_ob\n\n";
95 
96     unless($no_header){
97       print {$fhs{3}->{CHH}} "Bismark methylation extractor version $version\n";
98     }
99   }
00 
01   my $methylation_call_strings_processed = 0;
02   my $line_count = 0;
03 
04   ### proceeding differently now for single-end or paired-end Bismark files
05 
06   ### PROCESSING SINGLE-END RESULT FILES
07   if ($single) {
08 
09     ### also proceeding differently now for SAM format or vanilla Bismark format files
10     if ($vanilla) {		# old vanilla Bismark output format
11       while (<IN>) {
12 	++$line_count;
13 	warn "Processed lines: $line_count\n" if ($line_count%500000==0);
14 	
15 	### $seq here is the chromosomal sequence (to use for the repeat analysis for example)
16 	my ($id,$strand,$chrom,$start,$seq,$meth_call,$read_conversion,$genome_conversion) = (split("\t"))[0,1,2,3,6,7,8,9];
17 	
18 	### 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
19 	### last position
20 	chomp $genome_conversion;
21 
22 	my $index;
23 	if ($meth_call) {
24 
25 	  if ($read_conversion eq 'CT' and $genome_conversion eq 'CT') { ## original top strand
26 	    $index = 0;
27 	  } elsif ($read_conversion eq 'GA' and $genome_conversion eq 'CT') { ## complementary to original top strand
28 	    $index = 1;
29 	  } elsif ($read_conversion eq 'CT' and $genome_conversion eq 'GA') { ## original bottom strand
30 	    $index = 3;
31 	  } elsif ($read_conversion eq 'GA' and $genome_conversion eq 'GA') { ## complementary to original bottom strand
32 	    $index = 2;
33 	  } else {
34 	    die "Unexpected combination of read and genome conversion: '$read_conversion' / '$genome_conversion'\n";
35 	  }
36 	
37 	  ### Clipping off the first <int> number of bases from the methylation call string as specified with --ignore <int>
38 	  if ($ignore) {
39 	    $meth_call = substr($meth_call,$ignore,length($meth_call)-$ignore);	
40 	
41 	    ### If we are clipping off some bases at the start we need to adjust the start position of the alignments accordingly!
42 	    if ($strand eq '+') {
43 	      $start += $ignore;
44 	    } elsif ($strand eq '-') {
45 	      $start += length($meth_call)-1; ## $meth_call is already shortened!
46 	    } else {
47 	      die "Alignment did not have proper strand information: $strand\n";
48 	    }
49 	  }
50 	  ### printing out the methylation state of every C in the read
51 	  print_individual_C_methylation_states_single_end($meth_call,$chrom,$start,$id,$strand,$index);
52 	
53 	  ++$methylation_call_strings_processed; # 1 per single-end result
54 	}
55       }
56     } else {		  # processing single-end SAM format (default)
57       while (<IN>) {
58 	### SAM format can either start with header lines (starting with @) or start with alignments directly
59 	if (/^\@/) {	     # skipping header lines (starting with @)
60 	  warn "skipping SAM header line:\t$_";
61 	  next;
62 	}
63 
64 	++$line_count;
65 	warn "Processed lines: $line_count\n" if ($line_count%500000==0);
66 	
67 	# example read in SAM format
68 	# 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
69 	###
70 
71 	# < 0.7.6 my ($id,$chrom,$start,$meth_call,$read_conversion,$genome_conversion) = (split("\t"))[0,2,3,13,14,15];
72 	# < 0.7.6 $meth_call =~ s/^XM:Z://;
73 	# < 0.7.6 $read_conversion =~ s/^XR:Z://;
74 	# < 0.7.6 $genome_conversion =~ s/^XG:Z://;	
75 
76 	my ($id,$chrom,$start,$cigar) = (split("\t"))[0,2,3,5];
77 
78 	### detecting the following SAM flags in case the SAM entry was shuffled by CRAM or Goby compression/decompression
79 	my $meth_call;	  ### Thanks to Zachary Zeno for this solution
80 	my $read_conversion;
81 	my $genome_conversion;
82 	
83 	while ( /(XM|XR|XG):Z:([^\t]+)/g ) {
84 	  my $tag = $1;
85 	  my $value = $2;
86 
87 	  if ($tag eq "XM") {
88 	    $meth_call = $value;
89 	    $meth_call =~ s/\r//;
90 	  } elsif ($tag eq "XR") {
91 	    $read_conversion = $value;
92 	    $read_conversion =~ s/\r//;
93 	  } elsif ($tag eq "XG") {
94 	    $genome_conversion = $value;
95 	    $genome_conversion =~ s/\r//;
96 	  }
97 	}
98 
99 	my $strand;
00 	chomp $genome_conversion;
01 	# print "$meth_call\n$read_conversion\n$genome_conversion\n";
02 	
03 	my $index;
04 	if ($meth_call) {
05 	  if ($read_conversion eq 'CT' and $genome_conversion eq 'CT') { ## original top strand
06 	    $index = 0;
07 	    $strand = '+';
08 	  } elsif ($read_conversion eq 'GA' and $genome_conversion eq 'CT') { ## complementary to original top strand
09 	    $index = 1;
10 	    $strand = '-';
11 	  } elsif ($read_conversion eq 'GA' and $genome_conversion eq 'GA') { ## complementary to original bottom strand
12 	    $index = 2;
13 	    $strand = '+';
14 	  } elsif ($read_conversion eq 'CT' and $genome_conversion eq 'GA') { ## original bottom strand
15 	    $index = 3;
16 	    $strand = '-';
17 	  } else {
18 	    die "Unexpected combination of read and genome conversion: '$read_conversion' / '$genome_conversion'\n";
19 	  }
20 	
21 	  ### If the read is in SAM format we need to reverse the methylation call if the read has been reverse-complemented for the output
22 	  if ($strand eq '-') {
23 	    $meth_call = reverse $meth_call;
24 	  }
25 	
26 	  ### Clipping off the first <int> number of bases from the methylation call string as specified with --ignore <int>
27 	  if ($ignore) {
28 	    # print "\n\n$meth_call\n";
29 	    $meth_call = substr($meth_call,$ignore,length($meth_call)-$ignore);	
30 	    # print "$meth_call\n";
31 	    ### If we are ignoring a part of the sequence we also need to adjust the cigar string accordingly
32 
33 	    my @len = split (/\D+/,$cigar); # storing the length per operation
34 	    my @ops = split (/\d+/,$cigar); # storing the operation
35 	    shift @ops;		# remove the empty first element
36 	    die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len == scalar @ops);
37 		
38 	    my @comp_cigar; # building an array with all CIGAR operations
39 	    foreach my $index (0..$#len) {
40 	      foreach (1..$len[$index]) {
41 		# print  "$ops[$index]";
42 		push @comp_cigar, $ops[$index];
43 	      }
44 	    }
45 	    # print "original CIGAR: $cigar\n";
46 	    # print "original CIGAR: @comp_cigar\n";
47 
48 	    ### If we are clipping off some bases at the start we need to adjust the start position of the alignments accordingly!
49 	    if ($strand eq '+') {
50 	
51 	      my $D_count = 0; # counting all deletions that affect the ignored genomic position, i.e. Deletions and insertions
52 	      my $I_count = 0;
53 
54 	      for (1..$ignore) {
55 		my $op = shift @comp_cigar; # adjusting composite CIGAR string by removing $ignore operations from the start
56 		# print "$_ deleted $op\n";
57 
58 		while ($op eq 'D') { # repeating this for deletions (D)
59 		  $D_count++;
60 		  $op = shift @comp_cigar;
61 		  # print "$_ deleted $op\n";
62 		}
63 		if ($op eq 'I') { # adjusting the genomic position for insertions (I)
64 		  $I_count++;
65 		}
66 	      }
67 	      $start += $ignore + $D_count - $I_count;
68 	      # print "start $start\t ignore: $ignore\t D count: $D_count I_count: $I_count\n";
69 	    } elsif ($strand eq '-') {
70 
71 	      for (1..$ignore) {
72 		my $op = pop @comp_cigar; # adjusting composite CIGAR string by removing $ignore operations, here the last value of the array
73 		while ($op eq 'D') { # repeating this for deletions (D)
74 		  $op = pop @comp_cigar;
75 		}
76 	      }
77 
78 	      ### For reverse strand alignments we need to determine the number of matching bases (M) or deletions (D) in the read from the CIGAR
79 	      ### string to be able to work out the starting position of the read which is on the 3' end of the sequence
80 	      my $MD_count = 0;	# counting all operations that affect the genomic position, i.e. M and D. Insertions do not affect the start position
81 	      foreach (@comp_cigar) {
82 		++$MD_count if ($_ eq 'M' or $_ eq 'D');
83 	      }
84 	      $start += $MD_count - 1;
85 	    }
86 	
87 	    ### reconstituting shortened CIGAR string
88 	    my $new_cigar;
89 	    my $count = 0;
90 	    my $last_op;
91 	    # print "ignore adjusted: @comp_cigar\n";
92 	    foreach my $op (@comp_cigar) {
93 	      unless (defined $last_op){
94 		$last_op = $op;
95 		++$count;
96 		next;
97 	      }
98 	      if ($last_op eq $op) {
99 		++$count;
00 	      } else {
01 		$new_cigar .= "$count$last_op";
02 		$last_op = $op;
03 		$count = 1;
04 	      }
05 	    }
06 	    $new_cigar .= "$count$last_op"; # appending the last operation and count
07 	    $cigar = $new_cigar;
08 	    # print "ignore adjusted scalar: $cigar\n";
09 	  }
10 	}
11 	### printing out the methylation state of every C in the read
12 	print_individual_C_methylation_states_single_end($meth_call,$chrom,$start,$id,$strand,$index,$cigar);
13 	
14 	++$methylation_call_strings_processed; # 1 per single-end result
15       }
16     }
17   }
18 
19   ### PROCESSING PAIRED-END RESULT FILES
20   elsif ($paired) {
21 
22     ### proceeding differently now for SAM format or vanilla Bismark format files
23     if ($vanilla) {	# old vanilla Bismark paired-end output format
24       while (<IN>) {
25 	++$line_count;
26 	warn "processed line: $line_count\n" if ($line_count%500000==0);
27 
28 	### $seq here is the chromosomal sequence (to use for the repeat analysis for example)
29 	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];
30 
31 	my $index;
32 	chomp $genome_conversion;
33 	
34 	if ($first_read_conversion eq 'CT' and $genome_conversion eq 'CT') {
35 	  $index = 0;		## this is OT
36 	} elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'GA') {
37 	  $index = 2;		## this is CTOB!!!
38 	} elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'CT') {
39 	  $index = 1;		## this is CTOT!!!
40 	} elsif ($first_read_conversion eq 'CT' and $genome_conversion eq 'GA') {
41 	  $index = 3;		## this is OB
42 	} else {
43 	  die "Unexpected combination of read and genome conversion: $first_read_conversion / $genome_conversion\n";
44 	}
45 	
46 	if ($meth_call_1 and $meth_call_2) {
47 	  ### Clipping off the first <int> number of bases from the methylation call strings as specified with '--ignore <int>'
48 	  if ($ignore) {
49 	    $meth_call_1 = substr($meth_call_1,$ignore,length($meth_call_1)-$ignore);
50 	    $meth_call_2 = substr($meth_call_2,$ignore,length($meth_call_2)-$ignore);
51 
52 	    ### we also need to adjust the start and end positions of the alignments accordingly if '--ignore' was specified
53 	    $start_read_1 += $ignore;
54 	    $end_read_2   -= $ignore;
55 	  }
56 	  my $end_read_1;
57 	  my $start_read_2;
58 
59 	  if ($strand eq '+') {
60 
61 	    $end_read_1 = $start_read_1+length($meth_call_1)-1;
62 	    $start_read_2 = $end_read_2-length($meth_call_2)+1;
63 		
64 	    ## we first pass the first read which is in + orientation on the forward strand
65 	    print_individual_C_methylation_states_paired_end_files($meth_call_1,$chrom,$start_read_1,$id,'+',$index,0,0);
66 	
67 	    # we next pass the second read which is in - orientation on the reverse strand
68 	    ### 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
69 	    print_individual_C_methylation_states_paired_end_files($meth_call_2,$chrom,$end_read_2,$id,'-',$index,$no_overlap,$end_read_1);
70 	  } else {
71 	
72 	    $end_read_1 = $start_read_1+length($meth_call_2)-1;	# read 1 is the second reported read!
73 	    $start_read_2 = $end_read_2-length($meth_call_1)+1;	# read 2 is the first reported read!
74 
75 	    ## we first pass the first read which is in - orientation on the reverse strand
76 	    print_individual_C_methylation_states_paired_end_files($meth_call_1,$chrom,$end_read_2,$id,'-',$index,0,0);
77 
78 	    # we next pass the second read which is in + orientation on the forward strand
79 	    ### 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
80 	    print_individual_C_methylation_states_paired_end_files($meth_call_2,$chrom,$start_read_1,$id,'+',$index,$no_overlap,$start_read_2);
81 	  }
82 	
83 	  $methylation_call_strings_processed += 2; # paired-end = 2 methylation calls
84 	}	
85       }
86     } else {	      # Bismark paired-end SAM output format (default)
87       while (<IN>) {
88 	### SAM format can either start with header lines (starting with @) or start with alignments directly
89 	if (/^\@/) {	     # skipping header lines (starting with @)
90 	  warn "skipping SAM header line:\t$_";
91 	  next;
92 	}
93 	
94 	++$line_count;
95 	warn "Processed lines: $line_count\n" if ($line_count%500000==0);
96 	
97 	# example paired-end reads in SAM format (2 consecutive lines)
98 	# 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
99 	# 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
00 	
01 	#  < 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];
02 
03 	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
04 	my $meth_call_1;
05 	my $first_read_conversion;
06 	my $genome_conversion;
07 	
08 	while ( /(XM|XR|XG):Z:([^\t]+)/g ) {
09 	  my $tag = $1;
10 	  my $value = $2;
11 
12 	  if ($tag eq "XM") {
13 	    $meth_call_1 = $value;
14 	    $meth_call_1 =~ s/\r//;
15 	  } elsif ($tag eq "XR") {
16 	    $first_read_conversion = $value;
17 	    $first_read_conversion =~ s/\r//;
18 	  } elsif ($tag eq "XG") {
19 	    $genome_conversion = $value;
20 	    $genome_conversion =~ s/\r//;
21 	  }
22 	}
23 
24 	$_ = <IN>;		# reading in the paired read
25 
26 	# < version 0.7.6 my ($id_2,$start_read_2,$meth_call_2,$second_read_conversion) = (split("\t"))[0,3,13,14];
27 	# < version 0.7.6 $meth_call_1 =~ s/^XM:Z://;
28 	# < version 0.7.6 $meth_call_2 =~ s/^XM:Z://;
29 	# < version 0.7.6 $first_read_conversion =~ s/^XR:Z://;
30 	# < version 0.7.6 $second_read_conversion =~ s/^XR:Z://;
31 
32 	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
33 
34 	my $meth_call_2;
35 	my $second_read_conversion;
36 	
37 	while ( /(XM|XR):Z:([^\t]+)/g ) {
38 	  my $tag = $1;
39 	  my $value = $2;
40 
41 	  if ($tag eq "XM") {
42 	    $meth_call_2 = $value;
43 	    $meth_call_2 =~ s/\r//;
44 	  } elsif ($tag eq "XR") {
45 	    $second_read_conversion = $value;
46 	    $second_read_conversion = s/\r//;
47 	  }
48 	}
49 	
50 	# < version 0.7.6 $genome_conversion =~ s/^XG:Z://;	
51 	chomp $genome_conversion; # in case it captured a new line character	
52 
53 	# print join ("\t",$meth_call_1,$meth_call_2,$first_read_conversion,$second_read_conversion,$genome_conversion),"\n";
54 
55 	my $index;
56 	my $strand;
57 
58 	if ($first_read_conversion eq 'CT' and $genome_conversion eq 'CT') {
59 	  $index = 0;		## this is OT
60 	  $strand = '+';
61 	} elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'CT') {
62 	  $index = 1;		## this is CTOT
63 	  $strand = '-';
64 	} elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'GA') {
65 	  $index = 2;		## this is CTOB
66 	  $strand = '+';
67 	} elsif ($first_read_conversion eq 'CT' and $genome_conversion eq 'GA') {
68 	  $index = 3;		## this is OB
69 	  $strand = '-';
70 	} else {
71 	  die "Unexpected combination of read and genome conversion: $first_read_conversion / $genome_conversion\n";
72 	}
73 
74 	### reversing the methylation call of the read that was reverse-complemented
75 	if ($strand eq '+') {
76 	  $meth_call_2 = reverse $meth_call_2;
77 	} else {
78 	  $meth_call_1 = reverse $meth_call_1;
79 	}
80 
81 	if ($meth_call_1 and $meth_call_2) {
82 
83 	  my $end_read_1;
84 	
85 	  ### READ 1
86 	  my @len_1 = split (/\D+/,$cigar_1); # storing the length per operation
87 	  my @ops_1 = split (/\d+/,$cigar_1); # storing the operation
88 	  shift @ops_1;		# remove the empty first element
89 	  die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len_1 == scalar @ops_1);
90 	
91 	  my @comp_cigar_1; # building an array with all CIGAR operations
92 	  foreach my $index (0..$#len_1) {
93 	    foreach (1..$len_1[$index]) {
94 	      # print  "$ops_1[$index]";
95 	      push @comp_cigar_1, $ops_1[$index];
96 	    }
97 	  }
98 	  # print "original CIGAR read 1: $cigar_1\n";
99 	  # print "original CIGAR read 1: @comp_cigar_1\n";
00 
01 	  ### READ 2
02 	  my @len_2 = split (/\D+/,$cigar_2); # storing the length per operation
03 	  my @ops_2 = split (/\d+/,$cigar_2); # storing the operation
04 	  shift @ops_2;		# remove the empty first element
05 	  die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len_2 == scalar @ops_2);
06 	  my @comp_cigar_2; # building an array with all CIGAR operations for read 2
07 	  foreach my $index (0..$#len_2) {
08 	    foreach (1..$len_2[$index]) {
09 	      # print  "$ops_2[$index]";
10 	      push @comp_cigar_2, $ops_2[$index];
11 	    }
12 	  }
13 	  # print "original CIGAR read 2: $cigar_2\n";
14 	  # print "original CIGAR read 2: @comp_cigar_2\n";
15 	
16 	  if ($ignore) {
17 	    ### Clipping off the first <int> number of bases from the methylation call strings as specified with '--ignore <int>'	
18 	    ### the methylation calls have already been reversed where necessary
19 	    $meth_call_1 = substr($meth_call_1,$ignore,length($meth_call_1)-$ignore);
20 	    $meth_call_2 = substr($meth_call_2,$ignore,length($meth_call_2)-$ignore);
21 
22 	    ### If we are ignoring a part of the sequence we also need to adjust the cigar string accordingly
23 
24 	    if ($strand eq '+') {
25 		
26 	      ### if the (read 1) strand information is '+', read 1 needs to be trimmed from the start
27 	      my $D_count_1 = 0; # counting all deletions that affect the ignored genomic position for read 1, i.e. Deletions and insertions
28 	      my $I_count_1 = 0;
29 
30 	      for (1..$ignore) {
31 		my $op = shift @comp_cigar_1; # adjusting composite CIGAR string of read 1 by removing $ignore operations from the start
32 		# print "$_ deleted $op\n";
33 		
34 		while ($op eq 'D') { # repeating this for deletions (D)
35 		  $D_count_1++;
36 		  $op = shift @comp_cigar_1;
37 		  # print "$_ deleted $op\n";
38 		}
39 		if ($op eq 'I') { # adjusting the genomic position for insertions (I)
40 		  $I_count_1++;
41 		}
42 	      }
43 		
44 	      $start_read_1 += $ignore + $D_count_1 - $I_count_1;
45 	      # print "start read 1 $start_read_1\t ignore: $ignore\t D count 1: $D_count_1\tI_count 1: $I_count_1\n";
46 	
47 	      ### if the (read 1) strand information is '+', read 2 needs to be trimmed from the back
48 
49 	      for (1..$ignore) {
50 		my $op = pop @comp_cigar_2; # adjusting composite CIGAR string by removing $ignore operations, here the last value of the array
51 		while ($op eq 'D') { # repeating this for deletions (D)
52 		  $op = pop @comp_cigar_2;
53 		}
54 	      }
55 	      # the start position of reads mapping to the reverse strand is being adjusted further below
56 	    } elsif ($strand eq '-') {
57 	
58 	      ### if the (read 1) strand information is '-', read 1 needs to be trimmed from the back
59 	      for (1..$ignore) {
60 		my $op = pop @comp_cigar_1; # adjusting composite CIGAR string by removing $ignore operations, here the last value of the array
61 		while ($op eq 'D') { # repeating this for deletions (D)
62 		  $op = pop @comp_cigar_1;
63 		}
64 	      }
65 	      # the start position of reads mapping to the reverse strand is being adjusted further below
66 
67 	      ### if the (read 1) strand information is '-', read 2 needs to be trimmed from the start
68 	      my $D_count_2 = 0; # counting all deletions that affect the ignored genomic position for read 2, i.e. Deletions and insertions
69 	      my $I_count_2 = 0;
70 		
71 	      for (1..$ignore) {
72 		my $op = shift @comp_cigar_2; # adjusting composite CIGAR string of read 2 by removing $ignore operations from the start
73 		# print "$_ deleted $op\n";
74 		
75 		while ($op eq 'D') { # repeating this for deletions (D)
76 		  $D_count_2++;
77 		  $op = shift @comp_cigar_2;
78 		  # print "$_ deleted $op\n";
79 		}
80 		if ($op eq 'I') { # adjusting the genomic position for insertions (I)
81 		  $I_count_2++;
82 		}
83 	      }
84 		
85 	      $start_read_2 += $ignore + $D_count_2 - $I_count_2;
86 	      # print "start read 2 $start_read_2\t ignore: $ignore\t D count 2: $D_count_2\tI_count 2: $I_count_2\n";
87 	
88 	    }
89 	
90 	    ### reconstituting shortened CIGAR string 1
91 	    my $new_cigar_1;
92 	    my $count_1 = 0;
93 	    my $last_op_1;
94 	    # print "ignore adjusted CIGAR 1: @comp_cigar_1\n";
95 	    foreach my $op (@comp_cigar_1) {
96 	      unless (defined $last_op_1){
97 		$last_op_1 = $op;
98 		++$count_1;
99 		next;
00 	      }
01 	      if ($last_op_1 eq $op) {
02 		++$count_1;
03 	      } else {
04 		$new_cigar_1 .= "$count_1$last_op_1";
05 		$last_op_1 = $op;
06 		$count_1 = 1;
07 	      }
08 	    }
09 	    $new_cigar_1 .= "$count_1$last_op_1"; # appending the last operation and count
10 	    $cigar_1 = $new_cigar_1;
11 	    # print "ignore adjusted CIGAR 1 scalar: $cigar_1\n";
12 
13 	    ### reconstituting shortened CIGAR string 2
14 	    my $new_cigar_2;
15 	    my $count_2 = 0;
16 	    my $last_op_2;
17 	    # print "ignore adjusted CIGAR 2: @comp_cigar_2\n";
18 	    foreach my $op (@comp_cigar_2) {
19 	      unless (defined $last_op_2){
20 		$last_op_2 = $op;
21 		++$count_2;
22 		next;
23 	      }
24 	      if ($last_op_2 eq $op) {
25 		++$count_2;
26 	      } else {
27 		$new_cigar_2 .= "$count_2$last_op_2";
28 		$last_op_2 = $op;
29 		$count_2 = 1;
30 	      }
31 	    }
32 	    $new_cigar_2 .= "$count_2$last_op_2"; # appending the last operation and count
33 	    $cigar_2 = $new_cigar_2;
34 	    # print "ignore adjusted CIGAR 2 scalar: $cigar_2\n";
35 
36 	  }
37 	
38 	  if ($strand eq '+') {
39 	    ### adjusting the start position for all reads mapping to the reverse strand, in this case read 2
40 	    @comp_cigar_2  = reverse@comp_cigar_2; # the CIGAR string needs to be reversed for all reads aligning to the reverse strand, too
41 	    # print "reverse: @comp_cigar_2\n";
42 	
43 	    my $MD_count_1 = 0;
44 	    foreach (@comp_cigar_1) {
45 	      ++$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
46 	    }
47 
48 	    my $MD_count_2 = 0;
49 	    foreach (@comp_cigar_2) {
50 	      ++$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
51 	    }
52 
53 	    $end_read_1 = $start_read_1 + $MD_count_1 - 1;
54 	    $start_read_2 += $MD_count_2 - 1; ## Passing on the start position on the reverse strand
55 	  } else {
56 	    ### adjusting the start position for all reads mapping to the reverse strand, in this case read 1
57 	
58 	    @comp_cigar_1  = reverse@comp_cigar_1; # the CIGAR string needs to be reversed for all reads aligning to the reverse strand, too
59 	    # print "reverse: @comp_cigar_1\n";
60 
61 	    my $MD_count_1 = 0;
62 	    foreach (@comp_cigar_1) {
63 	      ++$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
64 	    }
65 
66 	    $end_read_1 = $start_read_1;	
67 	    $start_read_1 +=  $MD_count_1 - 1; ### Passing on the start position on the reverse strand
68 
69 	  }
70 
71 	  if ($strand eq '+') {
72 	    ## we first pass the first read which is in + orientation on the forward strand
73 	    print_individual_C_methylation_states_paired_end_files($meth_call_1,$chrom,$start_read_1,$id_1,'+',$index,0,0,$cigar_1);
74 	
75 	    # we next pass the second read which is in - orientation on the reverse strand
76 	    ### 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
77 	    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);
78 	  } else {
79 	    ## we first pass the first read which is in - orientation on the reverse strand
80 	    print_individual_C_methylation_states_paired_end_files($meth_call_1,$chrom,$start_read_1,$id_1,'-',$index,0,0,$cigar_1);
81 	
82 	    # we next pass the second read which is in + orientation on the forward strand
83 	    ### 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
84 	    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);
85 	  }
86 	
87 	  $methylation_call_strings_processed += 2; # paired-end = 2 methylation calls
88 	}	
89       }
90     }
91   } else {
92     die "Single-end or paired-end reads not specified properly\n";
93   }
94 
95   print "\n\nProcessed $line_count lines from $filename in total\n";
96   print "Total number of methylation call strings processed: $methylation_call_strings_processed\n\n";
97   if ($report) {
98     print REPORT "Total number of methylation call strings processed: $methylation_call_strings_processed\n\n";
99   }
00   print_splitting_report ();
01 }
02 
03 
04 
05 sub print_splitting_report{
06 
07   ### Calculating methylation percentages if applicable
08 
09   my $percent_meCpG;
10   if (($counting{total_meCpG_count}+$counting{total_unmethylated_CpG_count}) > 0){
11     $percent_meCpG = sprintf("%.1f",100*$counting{total_meCpG_count}/($counting{total_meCpG_count}+$counting{total_unmethylated_CpG_count}));
12   }
13 
14   my $percent_meCHG;
15   if (($counting{total_meCHG_count}+$counting{total_unmethylated_CHG_count}) > 0){
16     $percent_meCHG = sprintf("%.1f",100*$counting{total_meCHG_count}/($counting{total_meCHG_count}+$counting{total_unmethylated_CHG_count}));
17   }
18 
19   my $percent_meCHH;
20   if (($counting{total_meCHH_count}+$counting{total_unmethylated_CHH_count}) > 0){
21     $percent_meCHH = sprintf("%.1f",100*$counting{total_meCHH_count}/($counting{total_meCHH_count}+$counting{total_unmethylated_CHH_count}));
22   }
23 
24   my $percent_non_CpG_methylation;
25   if ($merge_non_CpG){
26     if ( ($counting{total_meCHH_count}+$counting{total_unmethylated_CHH_count}+$counting{total_meCHG_count}+$counting{total_unmethylated_CHG_count}) > 0){
27       $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} ) );
28     }
29   }
30 
31   if ($report){
32     ### detailed information about Cs analysed
33     print REPORT "Final Cytosine Methylation Report\n",'='x33,"\n";
34 
35     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};
36     print REPORT "Total number of C's analysed:\t$total_number_of_C\n\n";
37 
38     print REPORT "Total methylated C's in CpG context:\t$counting{total_meCpG_count}\n";
39     print REPORT "Total methylated C's in CHG context:\t$counting{total_meCHG_count}\n";
40     print REPORT "Total methylated C's in CHH context:\t$counting{total_meCHH_count}\n\n";
41 
42     print REPORT "Total C to T conversions in CpG context:\t$counting{total_unmethylated_CpG_count}\n";
43     print REPORT "Total C to T conversions in CHG context:\t$counting{total_unmethylated_CHG_count}\n";
44     print REPORT "Total C to T conversions in CHH context:\t$counting{total_unmethylated_CHH_count}\n\n";
45 
46     ### calculating methylated CpG percentage if applicable
47     if ($percent_meCpG){
48       print REPORT "C methylated in CpG context:\t${percent_meCpG}%\n";
49     }
50     else{
51       print REPORT "Can't determine percentage of methylated Cs in CpG context if value was 0\n";
52     }
53 
54     ### 2-Context Output
55     if ($merge_non_CpG){
56       if ($percent_non_CpG_methylation){
57 	print REPORT "C methylated in non-CpG context:\t${percent_non_CpG_methylation}%\n\n\n";
58       }
59       else{
60 	print REPORT "Can't determine percentage of methylated Cs in non-CpG context if value was 0\n\n\n";
61       }
62     }
63 
64     ### 3 Context Output
65     else{
66       ### calculating methylated CHG percentage if applicable
67       if ($percent_meCHG){
68 	print REPORT "C methylated in CHG context:\t${percent_meCHG}%\n";
69       }
70       else{
71 	print REPORT "Can't determine percentage of methylated Cs in CHG context if value was 0\n";
72       }
73 
74       ### calculating methylated CHH percentage if applicable
75       if ($percent_meCHH){
76 	print REPORT "C methylated in CHH context:\t${percent_meCHH}%\n\n\n";
77       }
78       else{
79 	print REPORT "Can't determine percentage of methylated Cs in CHH context if value was 0\n\n\n";
80       }
81     }
82   }
83 
84   ### detailed information about Cs analysed for on-screen report
85   print "Final Cytosine Methylation Report\n",'='x33,"\n";
86 
87   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};
88   print "Total number of C's analysed:\t$total_number_of_C\n\n";
89 
90   print "Total methylated C's in CpG context:\t$counting{total_meCpG_count}\n";
91   print "Total methylated C's in CHG context:\t$counting{total_meCHG_count}\n";
92   print "Total methylated C's in CHH context:\t$counting{total_meCHH_count}\n\n";
93 
94   print "Total C to T conversions in CpG context:\t$counting{total_unmethylated_CpG_count}\n";
95   print "Total C to T conversions in CHG context:\t$counting{total_unmethylated_CHG_count}\n";
96   print "Total C to T conversions in CHH context:\t$counting{total_unmethylated_CHH_count}\n\n";
97 
98   ### printing methylated CpG percentage if applicable
99   if ($percent_meCpG){
00     print "C methylated in CpG context:\t${percent_meCpG}%\n";
01   }
02   else{
03     print "Can't determine percentage of methylated Cs in CpG context if value was 0\n";
04   }
05 
06   ### 2-Context Output
07   if ($merge_non_CpG){
08     if ($percent_non_CpG_methylation){
09       print "C methylated in non-CpG context:\t${percent_non_CpG_methylation}%\n\n\n";
10     }
11     else{
12       print "Can't determine percentage of methylated Cs in non-CpG context if value was 0\n\n\n";
13     }
14   }
15 
16   ### 3-Context Output
17   else{
18     ### printing methylated CHG percentage if applicable
19     if ($percent_meCHG){
20       print "C methylated in CHG context:\t${percent_meCHG}%\n";
21     }
22     else{
23       print "Can't determine percentage of methylated Cs in CHG context if value was 0\n";
24     }
25 
26     ### printing methylated CHH percentage if applicable
27     if ($percent_meCHH){
28       print "C methylated in CHH context:\t${percent_meCHH}%\n\n\n";
29     }
30     else{
31       print "Can't determine percentage of methylated Cs in CHH context if value was 0\n\n\n";
32     }
33   }
34 }
35 
36 
37 
38 
39 
40 sub print_individual_C_methylation_states_paired_end_files{
41 
42   my ($meth_call,$chrom,$start,$id,$strand,$filehandle_index,$no_overlap,$end_read_1,$cigar) = @_;
43   my @methylation_calls = split(//,$meth_call);
44 
45   #################################################################
46   ### . for bases not involving cytosines                       ###
47   ### X for methylated C in CHG context (was protected)         ###
48   ### x for not methylated C in CHG context (was converted)     ###
49   ### H for methylated C in CHH context (was protected)         ###
50   ### h for not methylated C in CHH context (was converted)     ###
51   ### Z for methylated C in CpG context (was protected)         ###
52   ### z for not methylated C in CpG context (was converted)     ###
53   #################################################################
54 
55   my $methyl_CHG_count = 0;
56   my $methyl_CHH_count = 0;
57   my $methyl_CpG_count = 0;
58   my $unmethylated_CHG_count = 0;
59   my $unmethylated_CHH_count = 0;
60   my $unmethylated_CpG_count = 0;
61 
62   my @len;
63   my @ops;
64   my $pos_offset = 0; # this is only relevant for SAM reads with insertions or deletions
65   my $cigar_offset = 0; # again, this is only relevant for SAM reads containing indels
66   my @comp_cigar;
67 
68   if ($cigar){ # parsing CIGAR string
69     @len = split (/\D+/,$cigar); # storing the length per operation
70     @ops = split (/\d+/,$cigar); # storing the operation
71     shift @ops; # remove the empty first element
72     die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len == scalar @ops);
73 
74     foreach my $index (0..$#len){
75       foreach (1..$len[$index]){
76 	# print  "$ops[$index]";
77 	push @comp_cigar, $ops[$index];
78       }
79     }
80     # warn "\nDetected CIGAR string: $cigar\n";
81     # warn "Length of methylation call: ",length $meth_call,"\n";
82     # warn "number of operations: ",scalar @ops,"\n";
83     # warn "number of length digits: ",scalar @len,"\n\n";
84     # print @comp_cigar,"\n";
85     # print "$meth_call\n\n";
86     # sleep (1);
87   }
88 
89 
90   if ($strand eq '-') {
91 
92     ### the  CIGAR string needs to be reversed, the methylation call has already been reversed above
93     @comp_cigar  = reverse@comp_cigar; # the CIGAR string needs to be reversed for all reads aligning to the reverse strand, too
94     #  print "reverse CIGAR string: @comp_cigar\n";
95 
96     ### the start position of paired-end files has already been corrected, see above
97   }
98 
99   ### THIS IS AN OPTIONAL 2-CONTEXT (CpG and non-CpG) SECTION IF --merge_non_CpG was specified
00 
01   if ($merge_non_CpG) {
02 
03     if ($no_overlap) {
04 
05       ### single-file CpG and non-CpG context output
06       if ($full) {
07 	if ($strand eq '+') {
08 	  for my $index (0..$#methylation_calls) {
09 	
10 	    if ($cigar){ # only needed for SAM files
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{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
25 	    } elsif ($methylation_calls[$index] eq 'x') {
26 	      $counting{total_unmethylated_CHG_count}++;
27 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
28 	    } elsif ($methylation_calls[$index] eq 'Z') {
29 	      $counting{total_meCpG_count}++;
30 	      print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
31 	    } elsif ($methylation_calls[$index] eq 'z') {
32 	      $counting{total_unmethylated_CpG_count}++;
33 	      print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
34 	    } elsif ($methylation_calls[$index] eq 'H') {
35 	      $counting{total_meCHH_count}++;
36 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
37 	    } elsif ($methylation_calls[$index] eq 'h') {
38 	      $counting{total_unmethylated_CHH_count}++;
39 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
40 	    }
41 	    elsif ($methylation_calls[$index] eq '.'){}
42 	    else{
43 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
44 	    }
45 	  }
46 	} elsif ($strand eq '-') {
47 	  for my $index (0..$#methylation_calls) {
48 	
49 	    if ($cigar){ # only needed for SAM files
50 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
51 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
52 	      $cigar_offset += $cigar_mod;
53 	      $pos_offset += $pos_mod;
54 	    }
55 	
56 	    ### Returning as soon as the methylation calls start overlapping
57 	    if ($start-$index+$pos_offset <= $end_read_1) {
58 	      return;
59 	    }
60 
61 	    if ($methylation_calls[$index] eq 'X') {
62 	      $counting{total_meCHG_count}++;
63 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
64 	    } elsif ($methylation_calls[$index] eq 'x') {
65 	      $counting{total_unmethylated_CHG_count}++;
66 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
67 	    } elsif ($methylation_calls[$index] eq 'Z') {
68 	      $counting{total_meCpG_count}++;
69 	      print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
70 	    } elsif ($methylation_calls[$index] eq 'z') {
71 	      $counting{total_unmethylated_CpG_count}++;
72 	      print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
73 	    } elsif ($methylation_calls[$index] eq 'H') {
74 	      $counting{total_meCHH_count}++;
75 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
76 	    } elsif ($methylation_calls[$index] eq 'h') {
77 	      $counting{total_unmethylated_CHH_count}++;
78 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
79 	    }
80 	    elsif ($methylation_calls[$index] eq '.') {}
81 	    else{
82 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
83 	    }
84 	  }
85 	} else {
86 	  die "The read orientation was neither + nor -: '$strand'\n";
87 	}
88       }
89 
90       ### strand-specific methylation output
91       else {
92 	if ($strand eq '+') {
93 	  for my $index (0..$#methylation_calls) {
94 
95 	    if ($cigar){ # only needed for SAM files
96 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
97 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
98 	      $cigar_offset += $cigar_mod;
99 	      $pos_offset += $pos_mod;
00 	    }
01 
02 	    ### Returning as soon as the methylation calls start overlapping
03 	    if ($start+$index+$pos_offset >= $end_read_1) {
04 	      return;
05 	    }
06 
07 	    if ($methylation_calls[$index] eq 'X') {
08 	      $counting{total_meCHG_count}++;
09 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
10 	    } elsif ($methylation_calls[$index] eq 'x') {
11 	      $counting{total_unmethylated_CHG_count}++;
12 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
13 	    } elsif ($methylation_calls[$index] eq 'Z') {
14 	      $counting{total_meCpG_count}++;
15 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
16 	    } elsif ($methylation_calls[$index] eq 'z') {
17 	      $counting{total_unmethylated_CpG_count}++;
18 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
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";
22 	    } elsif ($methylation_calls[$index] eq 'h') {
23 	      $counting{total_unmethylated_CHH_count}++;
24 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
25 	    }
26 	    elsif ($methylation_calls[$index] eq '.') {}
27 	    else{
28 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
29 	    }	
30 	  }
31 	} elsif ($strand eq '-') {
32 	  for my $index (0..$#methylation_calls) {
33 
34 	    if ($cigar){ # only needed for SAM files
35 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
36 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
37 	      $cigar_offset += $cigar_mod;
38 	      $pos_offset += $pos_mod;
39 	    }
40 
41 	    ### Returning as soon as the methylation calls start overlapping
42 	    if ($start-$index+$pos_offset <= $end_read_1) {
43 	      return;
44 	    }
45 
46 	    if ($methylation_calls[$index] eq 'X') {
47 	      $counting{total_meCHG_count}++;
48 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
49 	    } elsif ($methylation_calls[$index] eq 'x') {
50 	      $counting{total_unmethylated_CHG_count}++;
51 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
52 	    } elsif ($methylation_calls[$index] eq 'Z') {
53 	      $counting{total_meCpG_count}++;
54 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
55 	    } elsif ($methylation_calls[$index] eq 'z') {
56 	      $counting{total_unmethylated_CpG_count}++;
57 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
58 	    } elsif ($methylation_calls[$index] eq 'H') {
59 	      $counting{total_meCHH_count}++;
60 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
61 	    } elsif ($methylation_calls[$index] eq 'h') {
62 	      $counting{total_unmethylated_CHH_count}++;
63 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
64 	    }
65 	    elsif ($methylation_calls[$index] eq '.') {}
66 	    else{
67 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
68 	    }
69 	  }
70 	} else {
71 	  die "The strand orientation was neither + nor -: '$strand'/n";
72 	}
73       }
74     }
75 
76     ### this is the default paired-end procedure allowing overlaps and using every single C position
77     ### Still within the 2-CONTEXT ONLY optional section
78     else {
79       ### single-file CpG and non-CpG context output
80       if ($full) {
81 	if ($strand eq '+') {
82 	  for my $index (0..$#methylation_calls) {
83 
84 	    if ($cigar){ # only needed for SAM files
85 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
86 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
87 	      $cigar_offset += $cigar_mod;
88 	      $pos_offset += $pos_mod;
89 	    }
90 
91 	    if ($methylation_calls[$index] eq 'X') {
92 	      $counting{total_meCHG_count}++;
93 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
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";
97 	    } elsif ($methylation_calls[$index] eq 'Z') {
98 	      $counting{total_meCpG_count}++;
99 	      print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
00 	    } elsif ($methylation_calls[$index] eq 'z') {
01 	      $counting{total_unmethylated_CpG_count}++;
02 	      print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
03 	    } elsif ($methylation_calls[$index] eq 'H') {
04 	      $counting{total_meCHH_count}++;
05 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
06 	    } elsif ($methylation_calls[$index] eq 'h') {
07 	      $counting{total_unmethylated_CHH_count}++;
08 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
09 	    }
10 	    elsif ($methylation_calls[$index] eq '.') {}
11 	    else{
12 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
13 	    }
14 	  }
15 	} elsif ($strand eq '-') {
16 	  for my $index (0..$#methylation_calls) {
17 
18 	    if ($cigar){ # only needed for SAM files
19 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
20 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
21 	      $cigar_offset += $cigar_mod;
22 	      $pos_offset += $pos_mod;
23 	    }
24 
25 	    if ($methylation_calls[$index] eq 'X') {
26 	      $counting{total_meCHG_count}++;
27 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
28 	    } elsif ($methylation_calls[$index] eq 'x') {
29 	      $counting{total_unmethylated_CHG_count}++;
30 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
31 	    } elsif ($methylation_calls[$index] eq 'Z') {
32 	      $counting{total_meCpG_count}++;
33 	      print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
34 	    } elsif ($methylation_calls[$index] eq 'z') {
35 	      $counting{total_unmethylated_CpG_count}++;
36 	      print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
37 	    } elsif ($methylation_calls[$index] eq 'H') {
38 	      $counting{total_meCHH_count}++;
39 	      print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
40 	    } elsif ($methylation_calls[$index] eq 'h') {
41 	      $counting{total_unmethylated_CHH_count}++;
42 	      print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
43 	    }
44 	    elsif ($methylation_calls[$index] eq '.') {}
45 	    else{
46 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
47 	    }
48 	  }
49 	} else {
50 	  die "The strand orientation as neither + nor -: '$strand'\n";
51 	}
52       }
53 
54       ### strand-specific methylation output
55       ### still within the 2-CONTEXT optional section
56       else {
57 	if ($strand eq '+') {
58 	  for my $index (0..$#methylation_calls) {
59 
60 	    if ($cigar){ # only needed for SAM files
61 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
62 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
63 	      $cigar_offset += $cigar_mod;
64 	      $pos_offset += $pos_mod;
65 	    }
66 
67 	    if ($methylation_calls[$index] eq 'X') {
68 	      $counting{total_meCHG_count}++;
69 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
70 	    } elsif ($methylation_calls[$index] eq 'x') {
71 	      $counting{total_unmethylated_CHG_count}++;
72 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
73 	    } elsif ($methylation_calls[$index] eq 'Z') {
74 	      $counting{total_meCpG_count}++;
75 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
76 	    } elsif ($methylation_calls[$index] eq 'z') {
77 	      $counting{total_unmethylated_CpG_count}++;
78 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
79 	    } elsif ($methylation_calls[$index] eq 'H') {
80 	      $counting{total_meCHH_count}++;
81 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
82 	    } elsif ($methylation_calls[$index] eq 'h') {
83 	      $counting{total_unmethylated_CHH_count}++;
84 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
85 	    }
86 	    elsif ($methylation_calls[$index] eq '.') {}
87 	    else{
88 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
89 	    }
90 	  }
91 	} elsif ($strand eq '-') {
92 	  for my $index (0..$#methylation_calls) {
93 
94 	    if ($cigar){ # only needed for SAM files
95 	      # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
96 	      my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
97 	      $cigar_offset += $cigar_mod;
98 	      $pos_offset += $pos_mod;
99 	    }
00 
01 	    if ($methylation_calls[$index] eq 'X') {
02 	      $counting{total_meCHG_count}++;
03 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
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";
07 	    } elsif ($methylation_calls[$index] eq 'Z') {
08 	      $counting{total_meCpG_count}++;
09 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
10 	    } elsif ($methylation_calls[$index] eq 'z') {
11 	      $counting{total_unmethylated_CpG_count}++;
12 	      print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
13 	    } elsif ($methylation_calls[$index] eq 'H') {
14 	      $counting{total_meCHH_count}++;
15 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
16 	    } elsif ($methylation_calls[$index] eq 'h') {
17 	      $counting{total_unmethylated_CHH_count}++;
18 	      print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
19 	    }
20 	    elsif ($methylation_calls[$index] eq '.') {}
21 	    else{
22 	      die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
23 	    }
24 	  }
25 	} else {
26 	  die "The strand orientation as neither + nor -: '$strand'\n";
27 	}
28       }
29     }
30   }
31 
32   ############################################
33   ### THIS IS THE DEFAULT 3-CONTEXT OUTPUT ###
34   ############################################
35 
36   elsif ($no_overlap) {
37     ### single-file CpG, CHG and CHH context output
38     if ($full) {
39       if ($strand eq '+') {
40 	for my $index (0..$#methylation_calls) {
41 	
42 	  if ($cigar){ # only needed for SAM files
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{CHG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
57 	  } elsif ($methylation_calls[$index] eq 'x') {
58 	    $counting{total_unmethylated_CHG_count}++;
59 	    print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
60 	  } elsif ($methylation_calls[$index] eq 'Z') {
61 	    $counting{total_meCpG_count}++;
62 	    print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
63 	  } elsif ($methylation_calls[$index] eq 'z') {
64 	    $counting{total_unmethylated_CpG_count}++;
65 	    print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
66 	  } elsif ($methylation_calls[$index] eq 'H') {
67 	    $counting{total_meCHH_count}++;
68 	    print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
69 	  } elsif ($methylation_calls[$index] eq 'h') {
70 	    $counting{total_unmethylated_CHH_count}++;
71 	    print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
72 	  }
73 	  elsif ($methylation_calls[$index] eq '.') {}
74 	  else{
75 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
76 	  }
77 	}
78       } elsif ($strand eq '-') {
79 	for my $index (0..$#methylation_calls) {
80 
81 	  if ($cigar){ # only needed for SAM files
82 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
83 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
84 	    $cigar_offset += $cigar_mod;
85 	    $pos_offset += $pos_mod;
86 	  }
87 
88 	  ### Returning as soon as the methylation calls start overlapping
89 	  if ($start-$index+$pos_offset <= $end_read_1) {
90 	    return;
91 	  }
92 	
93 	  if ($methylation_calls[$index] eq 'X') {
94 	    $counting{total_meCHG_count}++;
95 	    print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
96 	  } elsif ($methylation_calls[$index] eq 'x') {
97 	    $counting{total_unmethylated_CHG_count}++;
98 	    print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
99 	  } elsif ($methylation_calls[$index] eq 'Z') {
00 	    $counting{total_meCpG_count}++;
01 	    print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
02 	  } elsif ($methylation_calls[$index] eq 'z') {
03 	    $counting{total_unmethylated_CpG_count}++;
04 	    print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
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";
08 	  } elsif ($methylation_calls[$index] eq 'h') {
09 	    $counting{total_unmethylated_CHH_count}++;
10 	    print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
11 	  }
12 	  elsif ($methylation_calls[$index] eq '.') {}
13 	  else{
14 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
15 	  }
16 	}
17       } else {
18 	die "The strand orientation as neither + nor -: '$strand'\n";
19       }
20     }
21 
22     ### strand-specific methylation output
23     else {
24       if ($strand eq '+') {
25 	for my $index (0..$#methylation_calls) {
26 
27 	  if ($cigar){ # only needed for SAM files
28 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
29 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
30 	    $cigar_offset += $cigar_mod;
31 	    $pos_offset += $pos_mod;
32 	  }
33 
34 	  ### Returning as soon as the methylation calls start overlapping
35 	  if ($start+$index+$pos_offset >= $end_read_1) {
36 	    return;
37 	  }
38 	
39 	  if ($methylation_calls[$index] eq 'X') {
40 	    $counting{total_meCHG_count}++;
41 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
42 	  } elsif ($methylation_calls[$index] eq 'x') {
43 	    $counting{total_unmethylated_CHG_count}++;
44 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
45 	  } elsif ($methylation_calls[$index] eq 'Z') {
46 	    $counting{total_meCpG_count}++;
47 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
48 	  } elsif ($methylation_calls[$index] eq 'z') {
49 	    $counting{total_unmethylated_CpG_count}++;
50 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
51 	  } elsif ($methylation_calls[$index] eq 'H') {
52 	    $counting{total_meCHH_count}++;
53 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
54 	  } elsif ($methylation_calls[$index] eq 'h') {
55 	    $counting{total_unmethylated_CHH_count}++;
56 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
57 	  }
58 	  elsif ($methylation_calls[$index] eq '.') {}
59 	  else{
60 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
61 	  }	
62 	}
63       } elsif ($strand eq '-') {
64 	for my $index (0..$#methylation_calls) {
65 
66 	  if ($cigar){ # only needed for SAM files
67 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
68 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
69 	    $cigar_offset += $cigar_mod;
70 	    $pos_offset += $pos_mod;
71 	  }
72 
73 	  ### Returning as soon as the methylation calls start overlapping
74 	  if ($start-$index+$pos_offset <= $end_read_1) {
75 	    return;
76 	  }
77 	
78 	  if ($methylation_calls[$index] eq 'X') {
79 	    $counting{total_meCHG_count}++;
80 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
81 	  } elsif ($methylation_calls[$index] eq 'x') {
82 	    $counting{total_unmethylated_CHG_count}++;
83 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
84 	  } elsif ($methylation_calls[$index] eq 'Z') {
85 	    $counting{total_meCpG_count}++;
86 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
87 	  } elsif ($methylation_calls[$index] eq 'z') {
88 	    $counting{total_unmethylated_CpG_count}++;
89 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
90 	  } elsif ($methylation_calls[$index] eq 'H') {
91 	    $counting{total_meCHH_count}++;
92 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
93 	  } elsif ($methylation_calls[$index] eq 'h') {
94 	    $counting{total_unmethylated_CHH_count}++;
95 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
96 	  }
97 	  elsif ($methylation_calls[$index] eq '.') {}
98 	  else{
99 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
00 	  }
01 	}
02       } else {
03 	die "The strand orientation as neither + nor -: '$strand'\n";
04       }
05     }
06   }
07 
08   ### this is the default paired-end procedure allowing overlaps and using every single C position
09   else {
10     ### single-file CpG, CHG and CHH context output
11     if ($full) {
12       if ($strand eq '+') {
13 	for my $index (0..$#methylation_calls) {
14 
15 	  if ($cigar){ # only needed for SAM files
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{CHG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
25 	  } elsif ($methylation_calls[$index] eq 'x') {
26 	    $counting{total_unmethylated_CHG_count}++;
27 	    print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
28 	  } elsif ($methylation_calls[$index] eq 'Z') {
29 	    $counting{total_meCpG_count}++;
30 	    print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
31 	  } elsif ($methylation_calls[$index] eq 'z') {
32 	    $counting{total_unmethylated_CpG_count}++;
33 	    print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
34 	  } elsif ($methylation_calls[$index] eq 'H') {
35 	    $counting{total_meCHH_count}++;
36 	    print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
37 	  } elsif ($methylation_calls[$index] eq 'h') {
38 	    $counting{total_unmethylated_CHH_count}++;
39 	    print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
40 	  }
41 	  elsif ($methylation_calls[$index] eq '.') {}
42 	  else{
43 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
44 	  }
45 	}
46       } elsif ($strand eq '-') {
47 	for my $index (0..$#methylation_calls) {
48 
49 	  if ($cigar){ # only needed for SAM files
50 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
51 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
52 	    $cigar_offset += $cigar_mod;
53 	    $pos_offset += $pos_mod;
54 	  }
55 
56 	  if ($methylation_calls[$index] eq 'X') {
57 	    $counting{total_meCHG_count}++;
58 	    print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
59 	  } elsif ($methylation_calls[$index] eq 'x') {
60 	    $counting{total_unmethylated_CHG_count}++;
61 	    print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
62 	  } elsif ($methylation_calls[$index] eq 'Z') {
63 	    $counting{total_meCpG_count}++;
64 	    print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
65 	  } elsif ($methylation_calls[$index] eq 'z') {
66 	    $counting{total_unmethylated_CpG_count}++;
67 	    print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
68 	  } elsif ($methylation_calls[$index] eq 'H') {
69 	    $counting{total_meCHH_count}++;
70 	    print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
71 	  } elsif ($methylation_calls[$index] eq 'h') {
72 	    $counting{total_unmethylated_CHH_count}++;
73 	    print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
74 	  }
75 	  elsif ($methylation_calls[$index] eq '.') {}
76 	  else{
77 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
78 	  }
79 	}
80       } else {
81 	die "The strand orientation as neither + nor -: '$strand'\n";
82       }
83     }
84 
85     ### strand-specific methylation output
86     else {
87       if ($strand eq '+') {
88 	for my $index (0..$#methylation_calls) {
89 
90 	  if ($cigar){ # only needed for SAM files
91 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
92 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
93 	    $cigar_offset += $cigar_mod;
94 	    $pos_offset += $pos_mod;
95 	  }
96 
97 	  if ($methylation_calls[$index] eq 'X') {
98 	    $counting{total_meCHG_count}++;
99 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
00 	  } elsif ($methylation_calls[$index] eq 'x') {
01 	    $counting{total_unmethylated_CHG_count}++;
02 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
03 	  } elsif ($methylation_calls[$index] eq 'Z') {
04 	    $counting{total_meCpG_count}++;
05 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
06 	  } elsif ($methylation_calls[$index] eq 'z') {
07 	    $counting{total_unmethylated_CpG_count}++;
08 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
09 	  } elsif ($methylation_calls[$index] eq 'H') {
10 	    $counting{total_meCHH_count}++;
11 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
12 	  } elsif ($methylation_calls[$index] eq 'h') {
13 	    $counting{total_unmethylated_CHH_count}++;
14 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
15 	  }
16 	  elsif ($methylation_calls[$index] eq '.') {}
17 	  else{
18 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
19 	  }
20 	}
21       } elsif ($strand eq '-') {
22 	for my $index (0..$#methylation_calls) {
23 
24 	  if ($cigar){ # only needed for SAM files
25 	    # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
26 	    my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
27 	    $cigar_offset += $cigar_mod;
28 	    $pos_offset += $pos_mod;
29 	  }
30 
31 	  if ($methylation_calls[$index] eq 'X') {
32 	    $counting{total_meCHG_count}++;
33 	    print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
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";
37 	  } elsif ($methylation_calls[$index] eq 'Z') {
38 	    $counting{total_meCpG_count}++;
39 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
40 	  } elsif ($methylation_calls[$index] eq 'z') {
41 	    $counting{total_unmethylated_CpG_count}++;
42 	    print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
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";
46 	  } elsif ($methylation_calls[$index] eq 'h') {
47 	    $counting{total_unmethylated_CHH_count}++;
48 	    print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
49 	  }
50 	  elsif ($methylation_calls[$index] eq '.') {}
51 	  else{
52 	    die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
53 	  }
54 	}
55       } else {
56 	die "The strand orientation as neither + nor -: '$strand'\n";
57       }
58     }
59   }
60 }
61 
62 sub check_cigar_string {
63   my ($index,$cigar_offset,$pos_offset,$strand,$comp_cigar) = @_;
64   # print "$index\t$cigar_offset\t$pos_offset\t$strand\t";
65   my ($new_cigar_offset,$new_pos_offset) = (0,0);
66 
67   if ($strand eq '+') {
68     #  print "### $strand strand @$comp_cigar[$index + $cigar_offset]\t";
69 
70     if (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'M'){ # sequence position matches the genomic position
71       #  warn "position needs no adjustment\n";
72     }
73 
74     elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'I'){ # insertion in the read sequence
75       $new_pos_offset -= 1; # we need to subtract the length of inserted bases from the genomic position
76       # warn "adjusted genomic position by -1 bp (insertion)\n";
77     }
78 
79     elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'D'){ # deletion in the read sequence
80       $new_cigar_offset += 1; # the composite cigar string does no longer match the methylation call index
81       $new_pos_offset += 1; # we need to add the length of deleted bases to get the genomic position
82       # warn "adjusted genomic position by +1 bp (deletion). Now looping through the CIGAR string until we hit another M or I\n";
83 
84       while ( ($index + $cigar_offset + $new_cigar_offset)  < (scalar @$comp_cigar) ){	
85 	if (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'M'){ # sequence position matches the genomic position
86 	  #  warn "position needs no adjustment\n";
87 	  last;
88 	}
89 	elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'I'){
90 	  $new_pos_offset -= 1; # we need to subtract the length of inserted bases from the genomic position
91 	  # warn "adjusted genomic position by another -1 bp (insertion)\n";
92 	  last;
93 	}
94 	elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'D'){ # deletion in the read sequence
95 	  $new_cigar_offset += 1; # the composite cigar string does no longer match the methylation call index
96 	  $new_pos_offset += 1; # we need to add the length of deleted bases to get the genomic position
97 	  # warn "adjusted genomic position by another +1 bp (deletion)\n";
98 	}
99 	else{
00 	  die "The CIGAR string contained undefined operations in addition to 'M', 'I' and 'D': '@$comp_cigar[$index + $cigar_offset + $new_cigar_offset]'\n";
01 	}
02       }
03     }
04     else{
05       die "The CIGAR string contained undefined operations in addition to 'M', 'I' and 'D': '@$comp_cigar[$index + $cigar_offset + $new_cigar_offset]'\n";
06     }
07   }
08 
09   elsif ($strand eq '-') {
10     # print "### $strand strand @$comp_cigar[$index + $cigar_offset]\t";
11 
12     if (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'M'){ # sequence position matches the genomic position
13      # warn "position needs no adjustment\n";
14     }
15 
16     elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'I'){ # insertion in the read sequence
17       $new_pos_offset += 1; # we need to add the length of inserted bases to the genomic position
18       # warn "adjusted genomic position by +1 bp (insertion)\n";
19     }
20 
21     elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'D'){ # deletion in the read sequence
22       $new_cigar_offset += 1; # the composite cigar string does no longer match the methylation call index
23       $new_pos_offset -= 1; # we need to subtract the length of deleted bases to get the genomic position
24       # warn "adjusted genomic position by -1 bp (deletion). Now looping through the CIGAR string until we hit another M or I\n";
25 
26       while ( ($index + $cigar_offset + $new_cigar_offset)  < (scalar @$comp_cigar) ){	
27 	if (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'M'){ # sequence position matches the genomic position
28 	  # warn "Found new 'M' operation; position needs no adjustment\n";
29 	  last;
30 	}
31 	elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'I'){
32 	  $new_pos_offset += 1; # we need to subtract the length of inserted bases from the genomic position
33 	  # warn "Found new 'I' operation; adjusted genomic position by another +1 bp (insertion)\n";
34 	  last;
35 	}
36 	elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'D'){ # deletion in the read sequence
37 	  $new_cigar_offset += 1; # the composite cigar string does no longer match the methylation call index
38 	  $new_pos_offset -= 1; # we need to subtract the length of deleted bases to get the genomic position
39 	  # warn "adjusted genomic position by another -1 bp (deletion)\n";
40 	}
41 	else{
42 	  die "The CIGAR string contained undefined operations in addition to 'M', 'I' and 'D': '@$comp_cigar[$index + $cigar_offset + $new_cigar_offset]'\n";
43 	}
44       }
45     }
46     else{
47       die "The CIGAR string contained undefined operations in addition to 'M', 'I' and 'D': '@$comp_cigar[$index + $cigar_offset + $new_cigar_offset]'\n";
48     }
49   }
50   # print "new cigar offset: $new_cigar_offset\tnew pos offset: $new_pos_offset\n";
51   return ($new_cigar_offset,$new_pos_offset);
52 }
53 
54 sub print_individual_C_methylation_states_single_end{
55 
56   my ($meth_call,$chrom,$start,$id,$strand,$filehandle_index,$cigar) = @_;
57   my @methylation_calls = split(//,$meth_call);
58 
59   #################################################################
60   ### . for bases not involving cytosines                       ###
61   ### X for methylated C in CHG context (was protected)         ###
62   ### x for not methylated C in CHG context (was converted)     ###
63   ### H for methylated C in CHH context (was protected)         ###
64   ### h for not methylated C in CHH context (was converted)     ###
65   ### Z for methylated C in CpG context (was protected)         ###
66   ### z for not methylated C in CpG context (was converted)     ###
67   #################################################################
68 
69   my $methyl_CHG_count = 0;
70   my $methyl_CHH_count = 0;
71   my $methyl_CpG_count = 0;
72   my $unmethylated_CHG_count = 0;
73   my $unmethylated_CHH_count = 0;
74   my $unmethylated_CpG_count = 0;
75 
76   my @len;
77   my @ops;
78   my $pos_offset = 0; # this is only relevant for SAM reads with insertions or deletions
79   my $cigar_offset = 0; # again, this is only relevant for SAM reads containing indels
80 
81   my @comp_cigar;
82 
83   if ($cigar){ # parsing CIGAR string
84     @len = split (/\D+/,$cigar); # storing the length per operation
85     @ops = split (/\d+/,$cigar); # storing the operation
86     shift @ops; # remove the empty first element
87     die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len == scalar @ops);
88 
89     foreach my $index (0..$#len){
90       foreach (1..$len[$index]){
91 	# print  "$ops[$index]";
92 	push @comp_cigar, $ops[$index];
93       }
94     }
95     # warn "\nDetected CIGAR string: $cigar\n";
96     # warn "Length of methylation call: ",length $meth_call,"\n";
97     # warn "number of operations: ",scalar @ops,"\n";
98     # warn "number of length digits: ",scalar @len,"\n\n";
99     # print @comp_cigar,"\n";
00     # print "$meth_call\n\n";
01     # sleep (1);
02   }
03 
04   ### adjusting the start position for all reads mapping to the reverse strand
05   if ($strand eq '-') {
06 
07     @comp_cigar  = reverse@comp_cigar; # the CIGAR string needs to be reversed for all reads aligning to the reverse strand, too
08     # print @comp_cigar,"\n";
09 
10     unless ($ignore){  ### if --ignore was specified the start position has already been corrected
11 
12       if ($cigar){ ### SAM format
13 	my $MD_count = 0;
14 	foreach (@comp_cigar){
15 	  ++$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
16 	}
17 	$start += $MD_count - 1;
18       }
19       else{ ### vanilla format
20 	$start += length($meth_call)-1;
21       }
22     }
23   }
24 
25   ### THIS IS THE CpG and Non-CpG SECTION (OPTIONAL)
26 
27   ### single-file CpG and other-context output
28   if ($full and $merge_non_CpG) {
29     if ($strand eq '+') {
30       for my $index (0..$#methylation_calls) {
31 
32 	if ($cigar){ # only needed for SAM files
33 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
34 	  # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
35 	  $cigar_offset += $cigar_mod;
36 	  $pos_offset += $pos_mod;
37 	}
38 
39 	### methylated Cs (any context) will receive a forward (+) orientation
40 	### not methylated Cs (any context) will receive a reverse (-) orientation
41 	if ($methylation_calls[$index] eq 'X') {
42 	  $counting{total_meCHG_count}++;
43 	  print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
44 	}
45 	elsif ($methylation_calls[$index] eq 'x') {
46 	  $counting{total_unmethylated_CHG_count}++;
47 	  print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
48 	}
49 	elsif ($methylation_calls[$index] eq 'Z') {
50 	  $counting{total_meCpG_count}++;
51 	  print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
52 	}
53 	elsif ($methylation_calls[$index] eq 'z') {
54 	  $counting{total_unmethylated_CpG_count}++;
55 	  print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
56 	}
57 	elsif ($methylation_calls[$index] eq 'H') {
58 	  $counting{total_meCHH_count}++;
59 	  print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
60 	}
61 	elsif ($methylation_calls[$index] eq 'h') {
62 	  $counting{total_unmethylated_CHH_count}++;
63 	  print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
64 	}
65 	elsif ($methylation_calls[$index] eq '.') {
66 	}
67 	else{
68 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
69 	}
70       }
71     }
72     elsif ($strand eq '-') {
73 
74       for my $index (0..$#methylation_calls) {
75 	### methylated Cs (any context) will receive a forward (+) orientation
76 	### not methylated Cs (any context) will receive a reverse (-) orientation
77 
78 	if ($cigar){ # only needed for SAM files
79 	  # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
80 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
81 	  $cigar_offset += $cigar_mod;
82 	  $pos_offset += $pos_mod;
83 	}
84 
85 	if ($methylation_calls[$index] eq 'X') {
86 	  $counting{total_meCHG_count}++;
87 	  print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
88 	}
89 	elsif ($methylation_calls[$index] eq 'x') {
90 	  $counting{total_unmethylated_CHG_count}++;
91 	  print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
92 	}
93 	elsif ($methylation_calls[$index] eq 'Z') {
94 	  $counting{total_meCpG_count}++;
95 	  print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
96 	}
97 	elsif ($methylation_calls[$index] eq 'z') {
98 	  $counting{total_unmethylated_CpG_count}++;
99 	  print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
00 	}
01 	elsif ($methylation_calls[$index] eq 'H') {
02 	  $counting{total_meCHH_count}++;
03 	  print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
04 	}
05 	elsif ($methylation_calls[$index] eq 'h') {
06 	  $counting{total_unmethylated_CHH_count}++;
07 	  print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
08 	}
09 	elsif ($methylation_calls[$index] eq '.'){
10 	}
11 	else{
12 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
13 	}
14       }
15     }
16     else {
17       die "The strand information was neither + nor -: $strand\n";
18     }
19   }
20 
21   ### strand-specific methylation output
22   elsif ($merge_non_CpG) {
23     if ($strand eq '+') {
24       for my $index (0..$#methylation_calls) {
25 	### methylated Cs (any context) will receive a forward (+) orientation
26 	### not methylated Cs (any context) will receive a reverse (-) orientation
27 
28 	if ($cigar){ # only needed for SAM files
29 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
30 	  $cigar_offset += $cigar_mod;
31 	  $pos_offset += $pos_mod;
32 	}
33 
34 	if ($methylation_calls[$index] eq 'X') {
35 	  $counting{total_meCHG_count}++;
36 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
37 	}
38 	elsif ($methylation_calls[$index] eq 'x') {
39 	  $counting{total_unmethylated_CHG_count}++;
40 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
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";
45 	}
46 	elsif ($methylation_calls[$index] eq 'z') {
47 	  $counting{total_unmethylated_CpG_count}++;
48 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
49 	}
50 	elsif ($methylation_calls[$index] eq 'H') {
51 	  $counting{total_meCHH_count}++;
52 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
53 	}
54 	elsif ($methylation_calls[$index] eq 'h') {
55 	  $counting{total_unmethylated_CHH_count}++;
56 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
57 	}
58 	elsif ($methylation_calls[$index] eq '.') {
59 	}
60 	else{
61 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
62 	}
63       }
64     }
65     elsif ($strand eq '-') {
66 
67       for my $index (0..$#methylation_calls) {
68 	### methylated Cs (any context) will receive a forward (+) orientation
69 	### not methylated Cs (any context) will receive a reverse (-) orientation
70 
71     	if ($cigar){ # only needed for SAM files
72 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
73 	  $cigar_offset += $cigar_mod;
74 	  $pos_offset += $pos_mod;
75 	}
76 
77 	if ($methylation_calls[$index] eq 'X') {
78 	  $counting{total_meCHG_count}++;
79 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
80 	}
81 	elsif ($methylation_calls[$index] eq 'x') {
82 	  $counting{total_unmethylated_CHG_count}++;
83 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
84 	}
85 	elsif ($methylation_calls[$index] eq 'Z') {
86 	  $counting{total_meCpG_count}++;
87 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
88 	}
89 	elsif ($methylation_calls[$index] eq 'z') {
90 	  $counting{total_unmethylated_CpG_count}++;
91 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
92 	}
93 	elsif ($methylation_calls[$index] eq 'H') {
94 	  $counting{total_meCHH_count}++;
95 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
96 	}
97 	elsif ($methylation_calls[$index] eq 'h') {
98 	  $counting{total_unmethylated_CHH_count}++;
99 	  print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
00 	}
01 	elsif ($methylation_calls[$index] eq '.') {
02 	}
03 	else{
04 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
05 	}
06       }
07     }
08     else {
09       die "The strand information was neither + nor -: $strand\n";
10     }
11   }
12 
13   ### THIS IS THE 3-CONTEXT (CpG, CHG and CHH) DEFAULT SECTION
14 
15   elsif ($full) {
16     if ($strand eq '+') {
17       for my $index (0..$#methylation_calls) {
18 	### methylated Cs (any context) will receive a forward (+) orientation
19 	### not methylated Cs (any context) will receive a reverse (-) orientation
20 
21 	if ($cigar){ # only needed for SAM files
22 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
23 	  $cigar_offset += $cigar_mod;
24 	  $pos_offset += $pos_mod;
25 	}
26 
27 	if ($methylation_calls[$index] eq 'X') {
28 	  $counting{total_meCHG_count}++;
29 	  print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
30 	} elsif ($methylation_calls[$index] eq 'x') {
31 	  $counting{total_unmethylated_CHG_count}++;
32 	  print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
33 	} elsif ($methylation_calls[$index] eq 'Z') {
34 	  $counting{total_meCpG_count}++;
35 	  print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
36 	} elsif ($methylation_calls[$index] eq 'z') {
37 	  $counting{total_unmethylated_CpG_count}++;
38 	  print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
39 	} elsif ($methylation_calls[$index] eq 'H') {
40 	  $counting{total_meCHH_count}++;
41 	  print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
42 	} elsif ($methylation_calls[$index] eq 'h') {
43 	  $counting{total_unmethylated_CHH_count}++;
44 	  print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
45 	}
46 	elsif ($methylation_calls[$index] eq '.') {}
47 	else{
48 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
49 	}
50       }
51     }
52     elsif ($strand eq '-') {
53 
54       for my $index (0..$#methylation_calls) {
55 	### methylated Cs (any context) will receive a forward (+) orientation
56 	### not methylated Cs (any context) will receive a reverse (-) orientation
57 
58 	if ($cigar){ # only needed for SAM files
59 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
60 	  $cigar_offset += $cigar_mod;
61 	  $pos_offset += $pos_mod;
62 	}
63 	
64 	if ($methylation_calls[$index] eq 'X') {
65 	  $counting{total_meCHG_count}++;
66 	  print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
67 	} elsif ($methylation_calls[$index] eq 'x') {
68 	  $counting{total_unmethylated_CHG_count}++;
69 	  print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
70 	} elsif ($methylation_calls[$index] eq 'Z') {
71 	  $counting{total_meCpG_count}++;
72 	  print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
73 	} elsif ($methylation_calls[$index] eq 'z') {
74 	  $counting{total_unmethylated_CpG_count}++;
75 	  print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
76 	} elsif ($methylation_calls[$index] eq 'H') {
77 	  $counting{total_meCHH_count}++;
78 	  print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
79 	} elsif ($methylation_calls[$index] eq 'h') {
80 	  $counting{total_unmethylated_CHH_count}++;
81 	  print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
82 	}
83 	elsif ($methylation_calls[$index] eq '.') {}
84 	else{
85 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
86 	}
87       }
88     }
89     else {
90       die "The read had a strand orientation which was neither + nor -: $strand\n";
91     }
92   }
93 
94   ### strand-specific methylation output
95   else {
96     if ($strand eq '+') {
97       for my $index (0..$#methylation_calls) {
98 	### methylated Cs (any context) will receive a forward (+) orientation
99 	### not methylated Cs (any context) will receive a reverse (-) orientation
00 
01 	if ($cigar){ # only needed for SAM files
02 	  my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);	
03 	  $cigar_offset += $cigar_mod;
04 	  $pos_offset += $pos_mod;
05 	}
06 
07 	if ($methylation_calls[$index] eq 'X') {
08 	  $counting{total_meCHG_count}++;
09 	  print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
10 	} elsif ($methylation_calls[$index] eq 'x') {
11 	  $counting{total_unmethylated_CHG_count}++;
12 	  print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
13 	} elsif ($methylation_calls[$index] eq 'Z') {
14 	  $counting{total_meCpG_count}++;
15 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
16 	} elsif ($methylation_calls[$index] eq 'z') {
17 	  $counting{total_unmethylated_CpG_count}++;
18 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
19 	} elsif ($methylation_calls[$index] eq 'H') {
20 	  $counting{total_meCHH_count}++;
21 	  print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
22 	} elsif ($methylation_calls[$index] eq 'h') {
23 	  $counting{total_unmethylated_CHH_count}++;
24 	  print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
25 	}
26 	elsif ($methylation_calls[$index] eq '.') {}
27 	else{
28 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
29 	}
30       }
31     }
32     elsif ($strand eq '-') {
33 
34       for my $index (0..$#methylation_calls) {
35 	### methylated Cs (any context) will receive a forward (+) orientation
36 	### not methylated Cs (any context) will receive a reverse (-) orientation
37 
38 	if ($cigar){ # only needed for SAM files
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{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
47 	} elsif ($methylation_calls[$index] eq 'x') {
48 	  $counting{total_unmethylated_CHG_count}++;
49 	  print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
50 	} elsif ($methylation_calls[$index] eq 'Z') {
51 	  $counting{total_meCpG_count}++;
52 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
53 	} elsif ($methylation_calls[$index] eq 'z') {
54 	  $counting{total_unmethylated_CpG_count}++;
55 	  print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
56 	} elsif ($methylation_calls[$index] eq 'H') {
57 	  $counting{total_meCHH_count}++;
58 	  print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
59 	} elsif ($methylation_calls[$index] eq 'h') {
60 	  $counting{total_unmethylated_CHH_count}++;
61 	  print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
62 	}
63 	elsif ($methylation_calls[$index] eq '.') {}
64 	else{
65 	  die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
66 	}
67       }
68     }
69     else {
70       die "The strand information was neither + nor -: $strand\n";
71     }
72   }
73 }
74 
75 
76 
77 #######################################################################################################################################
78 ### bismark2bedGaph section - START
79 #######################################################################################################################################
80 
81 sub process_bedGraph_output{
82   warn  "="x64,"\n";
83   warn "Methylation information will now be written into a bedGraph file\n";
84   warn  "="x64,"\n\n";
85   sleep (2);
86 
87   ### Closing all filehandles so that the Bismark methtylation extractor output doesn't get truncated due to buffering issues
88   foreach my $fh (keys %fhs) {
89     if ($fh =~ /^[1230]$/) {
90       foreach my $context (keys %{$fhs{$fh}}) {
91 	close $fhs{$fh}->{$context} or die $!;
92       }
93     } else {
94       close $fhs{$fh} or die $!;
95     }
96   }
97 
98   ### deciding which files to use for bedGraph conversion
99   foreach my $filename (@sorting_files){
00     # warn "$filename\n";
01     if ($filename =~ /\//){ # if files are in a different output folder we extract the filename again
02       $filename =~ s/.*\///; # replacing everything up to the last slash in the filename
03       # warn "$filename\n";
04     }
05 
06     if ($CX_context){
07       push @bedfiles,$filename;
08     }
09     else{ ## CpG context only (default)
10       if ($filename =~ /^CpG_/){
11 	push @bedfiles,$filename;
12       }
13       else{
14 	# skipping CHH or CHG files
15       }
16     }
17   }
18 
19   warn "Using the following files as Input:\n";
20   print join ("\t",@bedfiles),"\n\n";
21   sleep (2);
22 
23   my %temp_fhs;
24   my @temp_files; # writing all context files (default CpG only) to these files prior to sorting
25 
26   ### changing to the output directory
27   unless ($output_dir eq ''){ # default
28     chdir $output_dir or die "Failed to change directory to $output_dir\n";
29     warn "Changed directory to $output_dir\n";
30   }
31 
32   foreach my $infile (@bedfiles) {
33 
34     if ($remove) {
35       warn "Now replacing whitespaces in the sequence ID field of the Bismark methylation extractor output $infile prior to bedGraph conversion\n\n";
36 
37       open (READ,$infile) or die $!;
38 
39       my $removed_spaces_outfile = $infile;
40       $removed_spaces_outfile =~ s/$/.spaces_removed.txt/;
41 
42       open (REM,'>',$output_dir.$removed_spaces_outfile) or die "Couldn't write to file $removed_spaces_outfile: $!\n";
43 
44       unless ($no_header){
45 	$_ = <READ>;		### Bismark version header
46 	print REM $_;		### Bismark version header
47       }
48 
49       while (<READ>) {
50 	chomp;
51 	my ($id,$strand,$chr,$pos,$context) = (split (/\t/));
52 	$id =~ s/\s+/_/g;
53 	print REM join ("\t",$id,$strand,$chr,$pos,$context),"\n";
54       }
55 	
56       close READ or die $!;
57       close REM or die $!;
58 
59       ### changing the infile name to the new file without spaces
60       $infile = $removed_spaces_outfile;
61     }
62 
63     warn "Now writing methylation information for file $infile to individual files for each chromosome\n";
64     open (IN,$infile) or die $!;
65 
66     ## always ignoring the version header
67     unless ($no_header){
68       $_ = <IN>;		### Bismark version header
69     }
70 	
71     while (<IN>) {
72       chomp;
73       my ($chr) = (split (/\t/))[2];
74 	
75       unless (exists $temp_fhs{$chr}) {
76 	open ($temp_fhs{$chr},'>','chr'.$chr.'.meth_extractor.temp') or die "Failed to open filehandle: $!";
77       }
78       print {$temp_fhs{$chr}} "$_\n";
79     }
80 
81     warn "Finished writing out individual chromosome files for $infile\n";
82   }
83   warn "\n";
84 
85   @temp_files = <*.meth_extractor.temp>;
86 
87   warn "Collecting temporary chromosome file information...\n";
88   sleep (1);
89   warn "processing the following input file(s):\n";
90   warn join ("\n",@temp_files),"\n\n";
91   sleep (1);
92 
93   foreach my $in (@temp_files) {
94     warn "Sorting input file $in by positions\n";
95     open my $ifh, "sort -k3,3 -k4,4n $in |" or die "Input file could not be sorted. $!";
96     # print "Chromosome\tStart Position\tEnd Position\tMethylation Percentage\n";
97 	
98     ############################################# m.a.bentley - moved the variables out of the while loop to hold the current line data {
99 
00     my $name;
01     my $meth_state;
02     my $chr = "";
03     my $pos = 0;
04     my $meth_state2;
05 
06     my $last_pos;
07     my $last_chr;
08 	
09     #############################################  }
10 	
11     while (my $line = <$ifh>) {
12       next if $line =~ /^Bismark/;
13       chomp $line;
14 
15       ########################################### m.a.bentley - (1) set the last_chr and last_pos variables early in the while loop, before the line split (2) removed unnecessary setting of same variables in if statement {
16 
17       $last_chr = $chr;
18       $last_pos = $pos;
19       ($name, $meth_state, $chr, $pos, $meth_state2) = split "\t", $line;
20 
21       if (($last_pos ne $pos) || ($last_chr ne $chr)) {
22 	generate_output($last_chr,$last_pos) if $methylcalls[2] > 0;
23 	@methylcalls = qw (0 0 0);
24       }
25 
26       #############################################  }
27 	
28       my $validated = validate_methylation_call($meth_state, $meth_state2);
29       unless($validated){
30 	warn "Methylation state of sequence ($name) in file ($in) on line $. is inconsistent (meth_state is $meth_state, meth_state2 = $meth_state2)\n";
31 	next;
32       }
33       if ($meth_state eq "+") {
34 	$methylcalls[0]++;
35 	$methylcalls[2]++;
36       } else {
37 	$methylcalls[1]++;
38 	$methylcalls[2]++;
39       }
40     }
41 
42     ############################################# m.a.bentley - set the last_chr and last_pos variables for the last line in the file (outside the while loop's scope using the method i've implemented) {
43 	
44     $last_chr = $chr;
45     $last_pos = $pos;
46     if ($methylcalls[2] > 0) {
47       generate_output($last_chr,$last_pos) if $methylcalls[2] > 0;
48     }
49     #############################################  }
50 	
51     close $ifh or die $!;
52 	
53     @methylcalls = qw (0 0 0); # resetting @methylcalls
54 
55     ### deleting temporary files
56     my $delete = unlink $in;
57     if ($delete) {
58       warn "Successfully deleted the temporary input file $in\n\n";
59     }
60     else {
61       warn "The temporary inputfile $in could not be deleted $!\n\n";
62     }
63   }
64 }
65 
66 sub generate_output{
67   my $methcount = $methylcalls[0];
68   my $nonmethcount = $methylcalls[1];
69   my $totalcount = $methylcalls[2];
70   my $last_chr = shift;
71   my $last_pos = shift;
72   croak "Should not be generating output if there's no reads to this region" unless $totalcount > 0;
73   croak "Total counts ($totalcount) is not the sum of the methylated ($methcount) and unmethylated ($nonmethcount) counts" if $totalcount != ($methcount + $nonmethcount);
74 
75   ############################################# m.a.bentley - declare a new variable 'bed_pos' to distinguish from bismark positions (-1) - previous scripts modified the last_pos variable earlier in the script leading to problems in meth % calculation {
76 
77   my $bed_pos = $last_pos -1; ### Bismark coordinates are 1 based whereas bedGraph coordinates are 0 based.
78   my $meth_percentage;
79   ($totalcount >= $coverage_threshold) ? ($meth_percentage = ($methcount/$totalcount) * 100) : ($meth_percentage = undef);
80   # $meth_percentage =~ s/(\.\d\d).+$/$1/ unless $meth_percentage =~ /^Below/;
81   if (defined $meth_percentage){
82     if ($counts){
83       print OUT "$last_chr\t$bed_pos\t$bed_pos\t$meth_percentage\t$methcount\t$nonmethcount\n";
84     }
85     else{
86       print OUT "$last_chr\t$bed_pos\t$bed_pos\t$meth_percentage\n";
87     }
88   }
89   #############################################  }
90 }
91 
92 sub validate_methylation_call{
93   my $meth_state = shift;
94   croak "Missing (+/-) methylation call" unless defined $meth_state;
95   my $meth_state2 = shift;
96   croak "Missing alphabetical methylation call" unless defined $meth_state2;
97   my $is_consistent;
98   ($meth_state2 =~ /^z/i) ? ($is_consistent = check_CpG_methylation_call($meth_state, $meth_state2)) 
99                           : ($is_consistent = check_nonCpG_methylation_call($meth_state,$meth_state2));
00   return 1 if $is_consistent;
01   return 0;
02 }
03 
04 sub check_CpG_methylation_call{
05   my $meth1 = shift;
06   my $meth2 = shift;
07   return 1 if($meth1 eq "+" && $meth2 eq "Z");
08   return 1 if($meth1 eq "-" && $meth2 eq "z");
09   return 0;
10 }
11 
12 sub check_nonCpG_methylation_call{
13   my $meth1 = shift;
14   my $meth2 = shift;
15   return 1 if($meth1 eq "+" && $meth2 eq "C");
16   return 1 if($meth1 eq "+" && $meth2 eq "X");
17   return 1 if($meth1 eq "+" && $meth2 eq "H");
18   return 1 if($meth1 eq "-" && $meth2 eq "c");
19   return 1 if($meth1 eq "-" && $meth2 eq "x");
20   return 1 if($meth1 eq "-" && $meth2 eq "h");
21   return 0;
22 }
23 
24 #######################################################################################################################################
25 ### bismark2bedGaph section - END
26 #######################################################################################################################################
27 
28 
29 
30 
31 
32 
33 #######################################################################################################################################
34 ### genome-wide cytosine methylation report - START
35 #######################################################################################################################################
36 
37 sub generate_genome_wide_cytosine_report {
38 
39   warn  "="x78,"\n";
40   warn "Methylation information will now be written into a genome-wide cytosine report\n";
41   warn  "="x78,"\n\n";
42   sleep (2);
43 
44   ### changing to the output directory again
45   unless ($output_dir eq ''){ # default
46     chdir $output_dir or die "Failed to change directory to $output_dir\n";
47     # warn "Changed directory to $output_dir\n";
48   }
49 
50   my $in = shift;
51   open (IN,$in) or die $!;
52 
53   my $cytosine_out = shift;
54 
55   if ($CX_context){
56     $cytosine_out =~ s/$/genome-wide_CX_report.txt/;
57   }
58   else{
59     $cytosine_out =~ s/$/genome_wide_CpG_report.txt/;
60   }
61 
62   ### note: we are still in the folder: $output_dir, so we do not have to include this into the open commands
63   unless ($split_by_chromosome){ ### writing all output to a single file (default)
64     open (CYT,'>',$cytosine_out) or die $!;
65     warn "Writing genome-wide cytosine report to: $cytosine_out\n";
66     sleep (3);
67   }
68 
69   my $last_chr;
70   my %chr; # storing reads for one chromosome at a time
71 
72   my $count = 0;
73   while (<IN>){
74     chomp;
75     ++$count;
76     my ($chr,$start,$end,undef,$meth,$nonmeth) = (split /\t/);
77 
78     # defining the first chromosome
79     unless (defined $last_chr){
80       $last_chr = $chr;
81       # warn "Storing all covered cytosine positions for chromosome: $chr\n";
82     }
83 
84     if ($chr eq $last_chr){
85       $chr{$chr}->{$start}->{meth} = $meth;
86       $chr{$chr}->{$start}->{nonmeth} = $nonmeth;
87     }
88     else{
89       warn "Writing cytosine reports for chromosome $last_chr (stored ",scalar keys %{$chr{$last_chr}}," different covered positions)\n";
90 
91       if ($split_by_chromosome){ ## writing output to 1 file per chromosome
92 	my $chromosome_out = $cytosine_out;
93 	$chromosome_out =~ s/txt$/chr${last_chr}.txt/;
94       open (CYT,'>',$chromosome_out) or die $!;
95     }
96 
97       while ( $chromosomes{$last_chr} =~ /([CG])/g){
98 	
99 	my $tri_nt = '';
00 	my $context = '';
01 	my $pos = pos$chromosomes{$last_chr};
02 	
03 	my $strand;
04 	my $meth = 0;
05 	my $nonmeth = 0;
06 	
07 	if ($1 eq 'C'){    # C on forward strand
08 	  $tri_nt = substr ($chromosomes{$last_chr},($pos-1),3);   # positions are 0-based!
09 	  $strand = '+';
10 	}
11 	elsif ($1 eq 'G'){ # C on reverse strand
12 	  $tri_nt = substr ($chromosomes{$last_chr},($pos-3),3);   # positions are 0-based!
13 	  $tri_nt = reverse $tri_nt;
14 	  $tri_nt =~ tr/ACTG/TGAC/;
15 	  $strand = '-';
16 	}
17 	next if (length$tri_nt < 3); # trinucleotide sequence could not be extracted
18 
19 	if (exists $chr{$last_chr}->{($pos-1)}){ # stored positions are 0-based!
20 	  $meth =  $chr{$last_chr}->{$pos-1}->{meth};
21 	  $nonmeth = $chr{$last_chr}->{$pos-1}->{nonmeth};
22 	}
23 
24 	### determining cytosine context	
25 	if ($tri_nt =~ /^CG/){
26 	  $context = 'CG';
27 	}
28 	elsif ($tri_nt =~ /^C.{1}G$/){
29 	  $context = 'CHG';
30 	}
31 	elsif ($tri_nt =~ /^C.{2}$/){
32 	  $context = 'CHH';
33 	}
34 	else{ # if the context can't be determined the positions will not be printed (it will equally not have been reported by Bismark)
35 	  warn "The sequence context could not be determined (found: '$tri_nt'). Skipping.\n";
36 	  next;
37 	}
38 
39 	if ($CpG_only){
40 	  if ($tri_nt =~ /^CG/){ # CpG context is the default
41 	    if ($zero){ # zero based coordinates
42 	      $pos -= 1;
43 	      print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
44 	    }
45 	    else{ # default
46 	      print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
47 	    }
48 	  }
49 	}
50 	else{ ## all cytosines, specified with --CX
51 	  if ($zero){ # zero based coordinates
52 	    $pos -= 1;
53 	    print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
54 	  }
55 	  else{ # default
56 	    print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
57 	  }
58 	}
59       }
60 
61       %chr = (); # resetting the hash
62 
63       # new first entry
64       $last_chr = $chr;
65       $chr{$chr}->{$start}->{meth} = $meth;
66       $chr{$chr}->{$start}->{nonmeth} = $nonmeth;
67     }
68   }
69 
70   # Last found chromosome
71 warn "Writing cytosine reports for chromosome $last_chr (stored ",scalar keys %{$chr{$last_chr}}," different covered positions)\n";
72 
73 if ($split_by_chromosome){ ## writing output to 1 file per chromosome
74   my $chromosome_out = $cytosine_out;
75   $chromosome_out =~ s/txt$/chr${last_chr}.txt/;
76   open (CYT,'>',$chromosome_out) or die $!;
77 }
78 
79   while ( $chromosomes{$last_chr} =~ /([CG])/g){
80 	
81     my $tri_nt;
82     my $context;
83     my $pos = pos$chromosomes{$last_chr};
84 
85     my $strand;
86     my $meth = 0;
87     my $nonmeth = 0;
88 
89     if ($1 eq 'C'){    # C on forward strand
90       $tri_nt = substr ($chromosomes{$last_chr},($pos-1),3);   # positions are 0-based!
91       $strand = '+';
92     }
93     elsif ($1 eq 'G'){ # C on reverse strand
94       $tri_nt = substr ($chromosomes{$last_chr},($pos-3),3);   # positions are 0-based!
95       $tri_nt = reverse $tri_nt;
96       $tri_nt =~ tr/ACTG/TGAC/;
97       $strand = '-';
98     }
99 
00     if (exists $chr{$last_chr}->{($pos-1)}){ # stored positions are 0-based!
01       $meth =  $chr{$last_chr}->{$pos-1}->{meth};
02       $nonmeth = $chr{$last_chr}->{$pos-1}->{nonmeth};
03     }
04 
05     next if (length$tri_nt < 3); # trinucleotide sequence could not be extracted
06 
07     ### determining cytosine context	
08     if ($tri_nt =~ /^CG/){
09       $context = 'CG';
10     }
11     elsif ($tri_nt =~ /^C.{1}G$/){
12       $context = 'CHG';
13     }
14     elsif ($tri_nt =~ /^C.{2}$/){
15       $context = 'CHH';
16     }
17     else{ # if the context can't be determined the positions will not be printed (it will equally not have been reported by Bismark)
18       warn "The cytosine context could not be determined (found: '$tri_nt'). Skipping.\n";
19       next;
20     }
21 	
22     if ($CpG_only){
23       if ($tri_nt =~ /^CG/){ # CpG context is the default
24 	if ($zero){ # zero-based coordinates
25 	  $pos -= 1;
26 	  print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
27 	}
28 	else{ # default
29 	  print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
30 	}
31       }
32     }
33     else{ ## all cytosines, specified with --CX
34       if ($zero){ # zero based coordinates
35 	$pos -= 1;
36 	print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
37       }
38       else{ # default
39 	print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
40       }
41     }
42   }
43   close CYT or die $!;
44 }
45 
46 
47 sub read_genome_into_memory{
48 
49   ## reading in and storing the specified genome in the %chromosomes hash
50   chdir ($genome_folder) or die "Can't move to $genome_folder: $!";
51   warn "Now reading in and storing sequence information of the genome specified in: $genome_folder\n\n";
52 
53   my @chromosome_filenames =  <*.fa>;
54 
55   ### if there aren't any genomic files with the extension .fa we will look for files with the extension .fasta
56   unless (@chromosome_filenames){
57     @chromosome_filenames =  <*.fasta>;
58   }
59   unless (@chromosome_filenames){
60     die "The specified genome folder $genome_folder does not contain any sequence files in FastA format (with .fa or .fasta file extensions)\n";
61   }
62 
63   foreach my $chromosome_filename (@chromosome_filenames){
64 
65     # skipping the tophat entire mouse genome fasta file
66     next if ($chromosome_filename eq 'Mus_musculus.NCBIM37.fa');
67 
68     open (CHR_IN,$chromosome_filename) or die "Failed to read from sequence file $chromosome_filename $!\n";
69     ### first line needs to be a fastA header
70     my $first_line = <CHR_IN>;
71     chomp $first_line;
72     $first_line =~ s/\r//; # removing /r carriage returns
73 
74     ### Extracting chromosome name from the FastA header
75     my $chromosome_name = extract_chromosome_name($first_line);
76 	
77     my $sequence;
78     while (<CHR_IN>){
79       chomp;
80       $_ =~ s/\r//; # removing /r carriage returns
81 
82       if ($_ =~ /^>/){
83 	### storing the previous chromosome in the %chromosomes hash, only relevant for Multi-Fasta-Files (MFA)
84 	if (exists $chromosomes{$chromosome_name}){
85 	  warn "chr $chromosome_name (",length $sequence ," bp)\n";
86 	  die "Exiting because chromosome name already exists. Please make sure all chromosomes have a unique name!\n";
87 	}
88 	else {
89 	  if (length($sequence) == 0){
90 	    warn "Chromosome $chromosome_name in the multi-fasta file $chromosome_filename did not contain any sequence information!\n";
91 	  }
92 	  warn "chr $chromosome_name (",length $sequence ," bp)\n";
93 	  $chromosomes{$chromosome_name} = $sequence;
94 	}
95 	### resetting the sequence variable
96 	$sequence = '';
97 	### setting new chromosome name
98 	$chromosome_name = extract_chromosome_name($_);
99       }
00       else{
01 	$sequence .= uc$_;
02       }
03     }
04 
05     if (exists $chromosomes{$chromosome_name}){
06       warn "chr $chromosome_name (",length $sequence ," bp)\t";
07       die "Exiting because chromosome name already exists. Please make sure all chromosomes have a unique name.\n";
08     }
09     else{
10       if (length($sequence) == 0){
11 	warn "Chromosome $chromosome_name in the file $chromosome_filename did not contain any sequence information!\n";
12       }
13       warn "chr $chromosome_name (",length $sequence ," bp)\n";
14       $chromosomes{$chromosome_name} = $sequence;
15     }
16   }
17   warn "\n";
18   chdir $parent_dir or die "Failed to move to directory $parent_dir\n";
19 }
20 
21 sub extract_chromosome_name {
22   ## Bowtie extracts the first string after the inition > in the FASTA file, so we are doing this as well
23   my $fasta_header = shift;
24   if ($fasta_header =~ s/^>//){
25     my ($chromosome_name) = split (/\s+/,$fasta_header);
26     return $chromosome_name;
27   }
28   else{
29     die "The specified chromosome ($fasta_header) file doesn't seem to be in FASTA format as required!\n";
30   }
31 }
32 
33 #######################################################################################################################################
34 ### genome-wide cytosine methylation report - END
35 #######################################################################################################################################
36 
37 
38 
39 
40 sub print_helpfile{
41 
42  print << 'HOW_TO';
43 
44 
45 DESCRIPTION
46 
47 The following is a brief description of all options to control the Bismark
48 methylation extractor. The script reads in a bisulfite read alignment results file 
49 produced by the Bismark bisulfite mapper and extracts the methylation information
50 for individual cytosines. This information is found in the methylation call field
51 which can contain the following characters:
52 
53        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
54        ~~~   X   for methylated C in CHG context (was protected)     ~~~
55        ~~~   x   for not methylated C CHG (was converted)            ~~~
56        ~~~   H   for methylated C in CHH context (was protected)     ~~~
57        ~~~   h   for not methylated C in CHH context (was converted) ~~~
58        ~~~   Z   for methylated C in CpG context (was protected)     ~~~
59        ~~~   z   for not methylated C in CpG context (was converted) ~~~
60        ~~~   .   for any bases not involving cytosines               ~~~
61        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
62 
63 The methylation extractor outputs result files for cytosines in CpG, CHG and CHH
64 context (this distinction is actually already made in Bismark itself). As the methylation
65 information for every C analysed can produce files which easily have tens or even hundreds of
66 millions of lines, file sizes can become very large and more difficult to handle. The C
67 methylation info additionally splits cytosine methylation calls up into one of the four possible
68 strands a given bisulfite read aligned against:
69 
70              OT      original top strand
71              CTOT    complementary to original top strand
72 
73              OB      original bottom strand
74              CTOB    complementary to original bottom strand
75 
76 Thus, by default twelve individual output files are being generated per input file (unless
77 --comprehensive is specified, see below). The output files can be imported into a genome
78 viewer, such as SeqMonk, and re-combined into a single data group if desired (in fact
79 unless the bisulfite reads were generated preserving directionality it doesn't make any
80 sense to look at the data in a strand-specific manner). Strand-specific output files can
81 optionally be skipped, in which case only three output files for CpG, CHG or CHH context
82 will be generated. For both the strand-specific and comprehensive outputs there is also
83 the option to merge both non-CpG contexts (CHG and CHH) into one single non-CpG context.
84 
85 
86 The output files are in the following format (tab delimited):
87 
88 <sequence_id>     <strand>      <chromosome>     <position>     <methylation call>
89 
90 
91 USAGE: methylation_extractor [options] <filenames>
92 
93 
94 ARGUMENTS:
95 
96 <filenames>              A space-separated list of Bismark result files in SAM format from
97                          which methylation information is extracted for every cytosine in
98                          the reads. For alignment files in the older custom Bismark output
99                          see option '--vanilla'.
00 
01 OPTIONS:
02 
03 -s/--single-end          Input file(s) are Bismark result file(s) generated from single-end
04                          read data. Specifying either --single-end or --paired-end is
05                          mandatory.
06 
07 -p/--paired-end          Input file(s) are Bismark result file(s) generated from paired-end
08                          read data. Specifying either --paired-end or --single-end is
09                          mandatory.
10 
11 --vanilla                The Bismark result input file(s) are in the old custom Bismark format
12                          (up to version 0.5.x) and not in SAM format which is the default as
13                          of Bismark version 0.6.x or higher. Default: OFF.
14 
15 --no_overlap             For paired-end reads it is theoretically possible that read_1 and
16                          read_2 overlap. This option avoids scoring overlapping methylation
17                          calls twice (only methylation calls of read 1 are used for in the process
18                          since read 1 has historically higher quality basecalls than read 2).
19                          Whilst this option removes a bias towards more methylation calls
20                          in the center of sequenced fragments it may de facto remove a sizable
21                          proportion of the data. This option is highly recommended for paired-end
22                          data.
23 
24 --ignore <int>           Ignore the first <int> bp at the 5' end of each read when processing the
25                          methylation call string. This can remove e.g. a restriction enzyme site
26                          at the start of each read.
27 
28 --comprehensive          Specifying this option will merge all four possible strand-specific 
29                          methylation info into context-dependent output files. The default 
30                          contexts are:
31                           - CpG context
32                           - CHG context
33                           - CHH context
34 
35 --merge_non_CpG          This will produce two output files (in --comprehensive mode) or eight
36                          strand-specific output files (default) for Cs in
37                           - CpG context
38                           - non-CpG context
39 
40 --report                 Prints out a short methylation summary as well as the paramaters used to run
41                          this script.
42 
43 --no_header              Suppresses the Bismark version header line in all output files for more convenient
44                          batch processing.
45 
46 -o/--output DIR          Allows specification of a different output directory (absolute or relative
47                          path). If not specified explicitely, the output will be written to the current directory.
48 
49 --version                Displays version information.
50 
51 -h/--help                Displays this help file and exits.
52 
53 
54 
55 bedGraph specific options:
56 
57 --bedGraph               After finishing the methylation extraction, the methylation output is written into a
58                          sorted bedGraph file that reports the position of a given cytosine and its methylation 
59                          state (in %, seem details below). The methylation extractor output is temporarily split up into
60                          temporary files, one per chromosome (written into the current directory or folder
61                          specified with -o/--output); these temp files are then used for sorting and deleted
62                          afterwards. By default, only cytosines in CpG context will be sorted. The option
63                          '--CX_context' may be used to report all cyosines irrespective of sequence context
64                          (this will take MUCH longer!).
65 
66 
67 --cutoff [threshold]     The minimum number of times a methylation state has to be seen for that nucleotide
68                          before its methylation percentage is reported. Default: 1.
69 
70 --remove_spaces          Replaces whitespaces in the sequence ID field with underscores to allow sorting.
71 
72 
73 --counts                 Adds two additional columns to the output file to enable further calculations:
74                              col 5: number of methylated calls
75                              col 6: number of unmethylated calls
76                          This option is required if '--cytosine_report' is specified (and will be set automatically if
77                          necessary).
78 
79 --CX/--CX_context        The sorted bedGraph output file contains information on every single cytosine that was covered
80                          in the experiment irrespective of its sequence context. This applies to both forward and
81                          reverse strands. Please be aware that this option may generate large temporary and output files
82                          and may take a long time to sort (up to many hours). Default: OFF.
83                          (i.e. Default = CpG context only).
84 
85 
86 
87 Genome-wide cytosine methylation report specific options:
88 
89 --cytosine_report        After the conversion to bedGraph has completed, the option '--cytosine_report' produces a
90                          genome-wide methylation report for all cytosines in the genome. By default, the output uses 1-based
91                          chromosome coordinates (zero-based cords are optional) and reports CpG context only (all
92                          cytosine context is optional). The output considers all Cs on both forward and reverse strands and
93                          reports their position, strand, trinucleotide content and methylation state (counts are 0 if not
94                          covered).
95 
96 --CX/--CX_context        The output file contains information on every single cytosine in the genome irrespective of
97                          its context. This applies to both forward and reverse strands. Please be aware that this will
98                          generate output files with > 1.1 billion lines for a mammalian genome such as human or mouse.
99                          Default: OFF (i.e. Default = CpG context only).
00 
01 --zero_based             Uses zero-based coordinates like used in e.g. bed files instead of 1-based coordinates. Default: OFF.
02 
03 --genome_folder <path>   Enter the genome folder you wish to use to extract sequences from (full path only). Accepted
04                          formats are FastA files ending with '.fa' or '.fasta'. Specifying a genome folder path is mandatory.
05 
06 --split_by_chromosome    Writes the output into individual files for each chromosome instead of a single output file. Files
07                          will be named to include the input filename and the chromosome number.
08 
09 
10 
11 OUTPUT:
12 
13 The bismark_methylation_extractor output is in the form:
14 ========================================================
15 <seq-ID>  <methylation state*>  <chromosome>  <start position (= end position)>  <methylation call>
16 
17 * Methylated cytosines receive a '+' orientation,
18 * Unmethylated cytosines receive a '-' orientation.
19 
20 
21 
22 The bedGraph output (optional) looks like this (tab-delimited):
23 ===============================================================
24 <chromosome>  <start position>  <end position>  <methylation percentage>
25 
26 
27 
28 The genome-wide cytosine methylation output file is tab-delimited in the following format:
29 ==========================================================================================
30 <chromosome>  <position>  <strand>  <count methylated>  <count non-methylated>  <C-context>  <trinucleotide context>
31 
32 
33 
34 This script was last modified on 02 Oct 2012.
35 
36 HOW_TO
37 }