Mercurial > repos > bjoern-gruening > bismark

comparison bismark_methylation_extractor @ 0:36d124f44c0a draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
inital commit
author bjoern-gruening
date Tue, 25 Dec 2012 05:45:46 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:36d124f44c0a
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 }
100
101 sleep (1);
102
103 if ($cytosine_report){
104 warn "\n\nSummarising genome-wide cytosine methylation report parameters:\n";
105 warn '='x63,"\n";
106 warn "Generating comprehensive genome-wide cytosine report (output format: <Chromosome> <Start Position> <End Position> <Methylation Percentage> )\n";
107
108
109 if ($CX_context){
110 warn "Reporting methylation for all cytosine contexts. Be aware that this will generate enormous files\n";
111 }
112 else{ # default
113 $CpG_only = 1;
114 warn "Reporting cytosine methylation in CpG context only (default)\n";
115 }
116
117 if ($split_by_chromosome){
118 warn "Splitting the cytosine report output up into individual files for each chromosome\n";
119 }
120
121 ### Zero-based coordinates
122 if ($zero){
123 warn "Using zero-based genomic coordinates (user-defined)\n";
124 }
125 else{ # default, 1-based coords
126 warn "Using 1-based genomic coordinates (default)\n";
127 }
128
129 ### GENOME folder
130 if ($genome_folder){
131 unless ($genome_folder =~/\/$/){
132 $genome_folder =~ s/$/\//;
133 }
134 warn "Genome folder was specified as $genome_folder\n";
135 }
136 else{
137 $genome_folder = '/data/public/Genomes/Mouse/NCBIM37/';
138 warn "Using the default genome folder /data/public/Genomes/Mouse/NCBIM37/\n";
139 }
140 sleep (1);
141 }
142 }
143
144 warn "\n";
145 sleep (5);
146
147 ######################################################
148 ### PROCESSING FILES
149 ######################################################
150
151 foreach my $filename (@filenames){
152 # resetting counters and filehandles
153 %fhs = ();
154 %counting =(
155 total_meCHG_count => 0,
156 total_meCHH_count => 0,
157 total_meCpG_count => 0,
158 total_unmethylated_CHG_count => 0,
159 total_unmethylated_CHH_count => 0,
160 total_unmethylated_CpG_count => 0,
161 sequences_count => 0,
162 );
163 @sorting_files = ();
164 @bedfiles = ();
165
166 process_Bismark_results_file($filename);
167
168 if ($bedGraph){
169 my $out = $filename;
170 $out =~ s/sam$//;
171 $out =~ s/txt$//;
172 $out =~ s/$/bedGraph/;
173
174 my $bedGraph_output = $out;
175 open (OUT,'>',$output_dir.$out) or die $!;
176 # warn "Writing bedGraph to file: $out\n";
177
178 process_bedGraph_output();
179 close OUT or die $!;
180
181 ### genome-wide cytosine methylation report requires bedGraph processing anyway
182 if ($cytosine_report){
183 my $cytosine_out = $out;
184 $cytosine_out =~ s/bedGraph$//;
185
186 read_genome_into_memory();
187 generate_genome_wide_cytosine_report($bedGraph_output,$cytosine_out);
188 }
189 }
190 }
191
192
193 sub process_commandline{
194 my $help;
195 my $single_end;
196 my $paired_end;
197 my $ignore;
198 my $genomic_fasta;
199 my $full;
200 my $report;
201 my $extractor_version;
202 my $no_overlap;
203 my $merge_non_CpG;
204 my $vanilla;
205 my $output_dir;
206 my $no_header;
207 my $bedGraph;
208 my $coverage_threshold = 1; # Minimum number of reads covering before calling methylation status
209 my $remove;
210 my $counts;
211 my $cytosine_report;
212 my $genome_folder;
213 my $zero;
214 my $CpG_only;
215 my $CX_context;
216 my $split_by_chromosome;
217
218
219 my $command_line = GetOptions ('help|man' => \$help,
220 'p|paired-end' => \$paired_end,
221 's|single-end' => \$single_end,
222 'fasta' => \$genomic_fasta,
223 'ignore=i' => \$ignore,
224 'comprehensive' => \$full,
225 'report' => \$report,
226 'version' => \$extractor_version,
227 'no_overlap' => \$no_overlap,
228 'merge_non_CpG' => \$merge_non_CpG,
229 'vanilla' => \$vanilla,
230 'o|output=s' => \$output_dir,
231 'no_header' => \$no_header,
232 'bedGraph' => \$bedGraph,
233 "cutoff=i" => \$coverage_threshold,
234 "remove_spaces" => \$remove,
235 "counts" => \$counts,
236 "cytosine_report" => \$cytosine_report,
237 'g|genome_folder=s' => \$genome_folder,
238 "zero_based" => \$zero,
239 "CX|CX_context" => \$CX_context,
240 "split_by_chromosome" => \$split_by_chromosome,
241 );
242
243 ### EXIT ON ERROR if there were errors with any of the supplied options
244 unless ($command_line){
245 die "Please respecify command line options\n";
246 }
247
248 ### HELPFILE
249 if ($help){
250 print_helpfile();
251 exit;
252 }
253
254 if ($extractor_version){
255 print << "VERSION";
256
257
258 Bismark Methylation Extractor
259
260 Bismark Extractor Version: $version Copyright 2010-12 Felix Krueger, Babraham Bioinformatics
261 www.bioinformatics.babraham.ac.uk/projects/bismark/
262
263
264 VERSION
265 exit;
266 }
267
268
269 ### no files provided
270 unless (@ARGV){
271 die "You need to provide one or more Bismark files to create an individual C methylation output. Please respecify!\n";
272 }
273 @filenames = @ARGV;
274
275 warn "\n *** Bismark methylation extractor version $version ***\n\n";
276
277 ### IGNORING <INT> bases at the start of the read when processing the methylation call string
278 unless ($ignore){
279 $ignore = 0;
280 }
281 ### PRINT A REPORT
282 unless ($report){
283 $report = 0;
284 }
285
286 ### OUTPUT DIR PATH
287 if ($output_dir){
288 unless ($output_dir =~ /\/$/){
289 $output_dir =~ s/$/\//;
290 }
291 }
292 else{
293 $output_dir = '';
294 }
295
296 ### NO HEADER
297 unless ($no_header){
298 $no_header = 0;
299 }
300
301 ### OLD (VANILLA) OUTPUT FORMAT
302 unless ($vanilla){
303 $vanilla = 0;
304 }
305
306 if ($single_end){
307 $paired_end = 0; ### SINGLE END ALIGNMENTS
308 }
309 elsif ($paired_end){
310 $single_end = 0; ### PAIRED-END ALIGNMENTS
311 }
312 else{
313 die "Please specify whether the supplied file(s) are in Bismark single-end or paired-end format\n\n";
314 }
315
316 ### NO OVERLAP
317 if ($no_overlap){
318 die "The option '--no_overlap' can only be specified for paired-end input!\n" unless ($paired_end);
319 }
320 else{
321 $no_overlap = 0;
322 }
323
324 ### COMPREHENSIVE OUTPUT
325 unless ($full){
326 $full = 0;
327 }
328
329 ### MERGE NON-CpG context
330 unless ($merge_non_CpG){
331 $merge_non_CpG = 0;
332 }
333
334 ### remove white spaces in read ID (needed for sorting using the sort command
335 unless ($remove){
336 $remove = 0;
337 }
338
339 ### COVERAGE THRESHOLD FOR gedGraph OUTPUT
340 unless (defined $coverage_threshold){
341 unless ($coverage_threshold > 0){
342 die "Please select a coverage greater than 0 (positive integers only)\n";
343 }
344 $coverage_threshold = 1;
345 }
346
347 if ($zero){
348 die "Option '--zero' is only available if '--cytosine_report' is specified as well. Please respecify\n" unless ($cytosine_report);
349 }
350
351 if ($CX_context){
352 die "Option '--CX_context' is only available if '--cytosine_report' or '--bedGraph' is specified as well. Please respecify\n" unless ($cytosine_report or $bedGraph);
353 }
354
355 if ($cytosine_report){
356
357 ### GENOME folder
358 if ($genome_folder){
359 unless ($genome_folder =~/\/$/){
360 $genome_folder =~ s/$/\//;
361 }
362 }
363 else{
364 die "Please specify a genome folder to proceed (full path only)\n";
365 }
366
367 unless ($bedGraph){
368 warn "Setting the option '--bedGraph' since this is required for the genome-wide cytosine report\n";
369 $bedGraph = 1;
370 }
371 unless ($counts){
372 warn "Setting the option '--counts' since this is required for the genome-wide cytosine report\n";
373 $counts = 1;
374 }
375 warn "\n";
376 }
377
378 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);
379 }
380
381
382 sub process_Bismark_results_file{
383 my $filename = shift;
384
385 warn "\nNow reading in Bismark result file $filename\n\n";
386
387 if ($filename =~ /\.gz$/) {
388 open (IN,"zcat $filename |") or die "Can't open gzipped file $filename: $!\n";
389 } else {
390 open (IN,$filename) or die "Can't open file $filename: $!\n";
391 }
392
393 ### Vanilla and SAM output need to read different numbers of header lines
394 if ($vanilla) {
395 my $bismark_version = <IN>; ## discarding the Bismark version info
396 chomp $bismark_version;
397 $bismark_version =~ s/\r//; # replaces \r line feed
398 $bismark_version =~ s/Bismark version: //;
399 if ($bismark_version =~ /^\@/) {
400 warn "Detected \@ as the first character of the version information. Is it possible that the file is in SAM format?\n\n";
401 sleep (2);
402 }
403
404 unless ($version eq $bismark_version){
405 die "The methylation extractor and Bismark itself need to be of the same version!\n\nVersions used:\nmethylation extractor: '$version'\nBismark: '$bismark_version'\n";
406 }
407 } else {
408 # If the read is in SAM format (default) it can either start with @ header lines or start with alignments directly.
409 # We are reading from it further down
410 }
411
412 my $output_filename = (split (/\//,$filename))[-1];
413
414 ### OPENING OUTPUT-FILEHANDLES
415 if ($report) {
416 my $report_filename = $output_filename;
417 $report_filename =~ s/[\.sam|\.txt]$//;
418 $report_filename =~ s/$/_splitting_report.txt/;
419 $report_filename = $output_dir . $report_filename;
420 open (REPORT,'>',$report_filename) or die "Failed to write to file $report_filename $!\n";
421 }
422
423 if ($report) {
424 print REPORT "$output_filename\n\n";
425 print REPORT "Parameters used to extract methylation information:\n";
426 if ($paired) {
427 if ($vanilla) {
428 print REPORT "Bismark result file: paired-end (vanilla Bismark format)\n";
429 } else {
430 print REPORT "Bismark result file: paired-end (SAM format)\n"; # default
431 }
432 }
433
434 if ($single) {
435 if ($vanilla) {
436 print REPORT "Bismark result file: single-end (vanilla Bismark format)\n";
437 } else {
438 print REPORT "Bismark result file: single-end (SAM format)\n"; # default
439 }
440 }
441
442 if ($ignore) {
443 print REPORT "Ignoring first $ignore bases\n";
444 }
445
446 if ($full) {
447 print REPORT "Output specified: comprehensive\n";
448 } else {
449 print REPORT "Output specified: strand-specific (default)\n";
450 }
451
452 if ($no_overlap) {
453 print REPORT "No overlapping methylation calls specified\n";
454 }
455 if ($genomic_fasta) {
456 print REPORT "Genomic equivalent sequences will be printed out in FastA format\n";
457 }
458 if ($merge_non_CpG) {
459 print REPORT "Methylation in CHG and CHH context will be merged into \"non-CpG context\" output\n";
460 }
461
462 print REPORT "\n";
463 }
464
465 ### CpG-context and non-CpG context. THIS SECTION IS OPTIONAL
466 ### if --comprehensive AND --merge_non_CpG was specified we are only writing out one CpG-context and one Any-Other-context result file
467 if ($full and $merge_non_CpG) {
468 my $cpg_output = my $other_c_output = $output_filename;
469 ### C in CpG context
470 $cpg_output =~ s/^/CpG_context_/;
471 $cpg_output =~ s/sam$/txt/;
472 $cpg_output =~ s/$/.txt/ unless ($cpg_output =~ /\.txt$/);
473 $cpg_output = $output_dir . $cpg_output;
474 push @sorting_files,$cpg_output;
475 open ($fhs{CpG_context},'>',$cpg_output) or die "Failed to write to $cpg_output $! \n";
476 print "Writing result file containing methylation information for C in CpG context to $cpg_output\n";
477
478 unless ($no_header) {
479 print {$fhs{CpG_context}} "Bismark methylation extractor version $version\n";
480 }
481
482 ### C in any other context than CpG
483 $other_c_output =~ s/^/Non_CpG_context_/;
484 $other_c_output =~ s/sam$/txt/;
485 $other_c_output =~ s/$/.txt/ unless ($other_c_output =~ /\.txt$/);
486 $other_c_output = $output_dir . $other_c_output;
487 push @sorting_files,$other_c_output;
488 open ($fhs{other_context},'>',$other_c_output) or die "Failed to write to $other_c_output $!\n";
489 print "Writing result file containing methylation information for C in any other context to $other_c_output\n";
490
491 unless ($no_header) {
492 print {$fhs{other_context}} "Bismark methylation extractor version $version\n";
493 }
494 }
495
496 ### 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
497 elsif ($merge_non_CpG) {
498
499 my $cpg_ot = my $cpg_ctot = my $cpg_ctob = my $cpg_ob = $output_filename;
500
501 ### For cytosines in CpG context
502 $cpg_ot =~ s/^/CpG_OT_/;
503 $cpg_ot =~ s/sam$/txt/;
504 $cpg_ot =~ s/$/.txt/ unless ($cpg_ot =~ /\.txt$/);
505 $cpg_ot = $output_dir . $cpg_ot;
506 push @sorting_files,$cpg_ot;
507 open ($fhs{0}->{CpG},'>',$cpg_ot) or die "Failed to write to $cpg_ot $!\n";
508 print "Writing result file containing methylation information for C in CpG context from the original top strand to $cpg_ot\n";
509
510 unless($no_header){
511 print {$fhs{0}->{CpG}} "Bismark methylation extractor version $version\n";
512 }
513
514 $cpg_ctot =~ s/^/CpG_CTOT_/;
515 $cpg_ctot =~ s/sam$/txt/;
516 $cpg_ctot =~ s/$/.txt/ unless ($cpg_ctot =~ /\.txt$/);
517 $cpg_ctot = $output_dir . $cpg_ctot;
518 push @sorting_files,$cpg_ctot;
519 open ($fhs{1}->{CpG},'>',$cpg_ctot) or die "Failed to write to $cpg_ctot $!\n";
520 print "Writing result file containing methylation information for C in CpG context from the complementary to original top strand to $cpg_ctot\n";
521
522 unless($no_header){
523 print {$fhs{1}->{CpG}} "Bismark methylation extractor version $version\n";
524 }
525
526 $cpg_ctob =~ s/^/CpG_CTOB_/;
527 $cpg_ctob =~ s/sam$/txt/;
528 $cpg_ctob =~ s/$/.txt/ unless ($cpg_ctob =~ /\.txt$/);
529 $cpg_ctob = $output_dir . $cpg_ctob;
530 push @sorting_files,$cpg_ctob;
531 open ($fhs{2}->{CpG},'>',$cpg_ctob) or die "Failed to write to $cpg_ctob $!\n";
532 print "Writing result file containing methylation information for C in CpG context from the complementary to original bottom strand to $cpg_ctob\n";
533
534 unless($no_header){
535 print {$fhs{2}->{CpG}} "Bismark methylation extractor version $version\n";
536 }
537
538 $cpg_ob =~ s/^/CpG_OB_/;
539 $cpg_ob =~ s/sam$/txt/;
540 $cpg_ob =~ s/$/.txt/ unless ($cpg_ob =~ /\.txt$/);
541 $cpg_ob = $output_dir . $cpg_ob;
542 push @sorting_files,$cpg_ob;
543 open ($fhs{3}->{CpG},'>',$cpg_ob) or die "Failed to write to $cpg_ob $!\n";
544 print "Writing result file containing methylation information for C in CpG context from the original bottom strand to $cpg_ob\n\n";
545
546 unless($no_header){
547 print {$fhs{3}->{CpG}} "Bismark methylation extractor version $version\n";
548 }
549
550 ### For cytosines in Non-CpG (CC, CT or CA) context
551 my $other_c_ot = my $other_c_ctot = my $other_c_ctob = my $other_c_ob = $output_filename;
552
553 $other_c_ot =~ s/^/Non_CpG_OT_/;
554 $other_c_ot =~ s/sam$/txt/;
555 $other_c_ot =~ s/$/.txt/ unless ($other_c_ot =~ /\.txt$/);
556 $other_c_ot = $output_dir . $other_c_ot;
557 push @sorting_files,$other_c_ot;
558 open ($fhs{0}->{other_c},'>',$other_c_ot) or die "Failed to write to $other_c_ot $!\n";
559 print "Writing result file containing methylation information for C in any other context from the original top strand to $other_c_ot\n";
560
561 unless($no_header){
562 print {$fhs{0}->{other_c}} "Bismark methylation extractor version $version\n";
563 }
564
565 $other_c_ctot =~ s/^/Non_CpG_CTOT_/;
566 $other_c_ctot =~ s/sam$/txt/;
567 $other_c_ctot =~ s/$/.txt/ unless ($other_c_ctot =~ /\.txt$/);
568 $other_c_ctot = $output_dir . $other_c_ctot;
569 push @sorting_files,$other_c_ctot;
570 open ($fhs{1}->{other_c},'>',$other_c_ctot) or die "Failed to write to $other_c_ctot $!\n";
571 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";
572
573 unless($no_header){
574 print {$fhs{1}->{other_c}} "Bismark methylation extractor version $version\n";
575 }
576
577 $other_c_ctob =~ s/^/Non_CpG_CTOB_/;
578 $other_c_ctob =~ s/sam$/txt/;
579 $other_c_ctob =~ s/$/.txt/ unless ($other_c_ctob =~ /\.txt$/);
580 $other_c_ctob = $output_dir . $other_c_ctob;
581 push @sorting_files,$other_c_ctob;
582 open ($fhs{2}->{other_c},'>',$other_c_ctob) or die "Failed to write to $other_c_ctob $!\n";
583 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";
584
585 unless($no_header){
586 print {$fhs{2}->{other_c}} "Bismark methylation extractor version $version\n";
587 }
588
589 $other_c_ob =~ s/^/Non_CpG_OB_/;
590 $other_c_ob =~ s/sam$/txt/;
591 $other_c_ob =~ s/$/.txt/ unless ($other_c_ob =~ /\.txt$/);
592 $other_c_ob = $output_dir . $other_c_ob;
593 push @sorting_files,$other_c_ob;
594 open ($fhs{3}->{other_c},'>',$other_c_ob) or die "Failed to write to $other_c_ob $!\n";
595 print "Writing result file containing methylation information for C in any other context from the original bottom strand to $other_c_ob\n\n";
596
597 unless($no_header){
598 print {$fhs{3}->{other_c}} "Bismark methylation extractor version $version\n";
599 }
600 }
601 ### THIS SECTION IS THE DEFAULT (CpG, CHG and CHH context)
602
603 ### if --comprehensive was specified we are only writing one file per context
604 elsif ($full) {
605 my $cpg_output = my $chg_output = my $chh_output = $output_filename;
606 ### C in CpG context
607 $cpg_output =~ s/^/CpG_context_/;
608 $cpg_output =~ s/sam$/txt/;
609 $cpg_output =~ s/$/.txt/ unless ($cpg_output =~ /\.txt$/);
610 $cpg_output = $output_dir . $cpg_output;
611 push @sorting_files,$cpg_output;
612 open ($fhs{CpG_context},'>',$cpg_output) or die "Failed to write to $cpg_output $! \n";
613 print "Writing result file containing methylation information for C in CpG context to $cpg_output\n";
614
615 unless($no_header){
616 print {$fhs{CpG_context}} "Bismark methylation extractor version $version\n";
617 }
618
619 ### C in CHG context
620 $chg_output =~ s/^/CHG_context_/;
621 $chg_output =~ s/sam$/txt/;
622 $chg_output =~ s/$/.txt/ unless ($chg_output =~ /\.txt$/);
623 $chg_output = $output_dir . $chg_output;
624 push @sorting_files,$chg_output;
625 open ($fhs{CHG_context},'>',$chg_output) or die "Failed to write to $chg_output $!\n";
626 print "Writing result file containing methylation information for C in CHG context to $chg_output\n";
627
628 unless($no_header){
629 print {$fhs{CHG_context}} "Bismark methylation extractor version $version\n";
630 }
631
632 ### C in CHH context
633 $chh_output =~ s/^/CHH_context_/;
634 $chh_output =~ s/sam$/txt/;
635 $chh_output =~ s/$/.txt/ unless ($chh_output =~ /\.txt$/);
636 $chh_output = $output_dir . $chh_output;
637 push @sorting_files, $chh_output;
638 open ($fhs{CHH_context},'>',$chh_output) or die "Failed to write to $chh_output $!\n";
639 print "Writing result file containing methylation information for C in CHH context to $chh_output\n";
640
641 unless($no_header){
642 print {$fhs{CHH_context}} "Bismark methylation extractor version $version\n";
643 }
644 }
645 ### else we will write out 12 different output files, depending on where the (first) unique best alignment was found
646 else {
647 my $cpg_ot = my $cpg_ctot = my $cpg_ctob = my $cpg_ob = $output_filename;
648
649 ### For cytosines in CpG context
650 $cpg_ot =~ s/^/CpG_OT_/;
651 $cpg_ot =~ s/sam$/txt/;
652 $cpg_ot =~ s/$/.txt/ unless ($cpg_ot =~ /\.txt$/);
653 $cpg_ot = $output_dir . $cpg_ot;
654 push @sorting_files,$cpg_ot;
655 open ($fhs{0}->{CpG},'>',$cpg_ot) or die "Failed to write to $cpg_ot $!\n";
656 print "Writing result file containing methylation information for C in CpG context from the original top strand to $cpg_ot\n";
657
658 unless($no_header){
659 print {$fhs{0}->{CpG}} "Bismark methylation extractor version $version\n";
660 }
661
662 $cpg_ctot =~ s/^/CpG_CTOT_/;
663 $cpg_ctot =~ s/sam$/txt/;
664 $cpg_ctot =~ s/$/.txt/ unless ($cpg_ctot =~ /\.txt$/);
665 $cpg_ctot = $output_dir . $cpg_ctot;
666 push @sorting_files,$cpg_ctot;
667 open ($fhs{1}->{CpG},'>',$cpg_ctot) or die "Failed to write to $cpg_ctot $!\n";
668 print "Writing result file containing methylation information for C in CpG context from the complementary to original top strand to $cpg_ctot\n";
669
670 unless($no_header){
671 print {$fhs{1}->{CpG}} "Bismark methylation extractor version $version\n";
672 }
673
674 $cpg_ctob =~ s/^/CpG_CTOB_/;
675 $cpg_ctob =~ s/sam$/txt/;
676 $cpg_ctob =~ s/$/.txt/ unless ($cpg_ctob =~ /\.txt$/);
677 $cpg_ctob = $output_dir . $cpg_ctob;
678 push @sorting_files,$cpg_ctob;
679 open ($fhs{2}->{CpG},'>',$cpg_ctob) or die "Failed to write to $cpg_ctob $!\n";
680 print "Writing result file containing methylation information for C in CpG context from the complementary to original bottom strand to $cpg_ctob\n";
681
682 unless($no_header){
683 print {$fhs{2}->{CpG}} "Bismark methylation extractor version $version\n";
684 }
685
686 $cpg_ob =~ s/^/CpG_OB_/;
687 $cpg_ob =~ s/sam$/txt/;
688 $cpg_ob =~ s/$/.txt/ unless ($cpg_ob =~ /\.txt$/);
689 $cpg_ob = $output_dir . $cpg_ob;
690 push @sorting_files,$cpg_ob;
691 open ($fhs{3}->{CpG},'>',$cpg_ob) or die "Failed to write to $cpg_ob $!\n";
692 print "Writing result file containing methylation information for C in CpG context from the original bottom strand to $cpg_ob\n\n";
693
694 unless($no_header){
695 print {$fhs{3}->{CpG}} "Bismark methylation extractor version $version\n";
696 }
697
698 ### For cytosines in CHG context
699 my $chg_ot = my $chg_ctot = my $chg_ctob = my $chg_ob = $output_filename;
700
701 $chg_ot =~ s/^/CHG_OT_/;
702 $chg_ot =~ s/sam$/txt/;
703 $chg_ot =~ s/$/.txt/ unless ($chg_ot =~ /\.txt$/);
704 $chg_ot = $output_dir . $chg_ot;
705 push @sorting_files,$chg_ot;
706 open ($fhs{0}->{CHG},'>',$chg_ot) or die "Failed to write to $chg_ot $!\n";
707 print "Writing result file containing methylation information for C in CHG context from the original top strand to $chg_ot\n";
708
709 unless($no_header){
710 print {$fhs{0}->{CHG}} "Bismark methylation extractor version $version\n";
711 }
712
713 $chg_ctot =~ s/^/CHG_CTOT_/;
714 $chg_ctot =~ s/sam$/txt/;
715 $chg_ctot =~ s/$/.txt/ unless ($chg_ctot =~ /\.txt$/);
716 $chg_ctot = $output_dir . $chg_ctot;
717 push @sorting_files,$chg_ctot;
718 open ($fhs{1}->{CHG},'>',$chg_ctot) or die "Failed to write to $chg_ctot $!\n";
719 print "Writing result file containing methylation information for C in CHG context from the complementary to original top strand to $chg_ctot\n";
720
721 unless($no_header){
722 print {$fhs{1}->{CHG}} "Bismark methylation extractor version $version\n";
723 }
724
725 $chg_ctob =~ s/^/CHG_CTOB_/;
726 $chg_ctob =~ s/sam$/txt/;
727 $chg_ctob =~ s/$/.txt/ unless ($chg_ctob =~ /\.txt$/);
728 $chg_ctob = $output_dir . $chg_ctob;
729 push @sorting_files,$chg_ctob;
730 open ($fhs{2}->{CHG},'>',$chg_ctob) or die "Failed to write to $chg_ctob $!\n";
731 print "Writing result file containing methylation information for C in CHG context from the complementary to original bottom strand to $chg_ctob\n";
732
733 unless($no_header){
734 print {$fhs{2}->{CHG}} "Bismark methylation extractor version $version\n";
735 }
736
737 $chg_ob =~ s/^/CHG_OB_/;
738 $chg_ob =~ s/sam$/txt/;
739 $chg_ob =~ s/$/.txt/ unless ($chg_ob =~ /\.txt$/);
740 $chg_ob = $output_dir . $chg_ob;
741 push @sorting_files,$chg_ob;
742 open ($fhs{3}->{CHG},'>',$chg_ob) or die "Failed to write to $chg_ob $!\n";
743 print "Writing result file containing methylation information for C in CHG context from the original bottom strand to $chg_ob\n\n";
744
745 unless($no_header){
746 print {$fhs{3}->{CHG}} "Bismark methylation extractor version $version\n";
747 }
748
749 ### For cytosines in CHH context
750 my $chh_ot = my $chh_ctot = my $chh_ctob = my $chh_ob = $output_filename;
751
752 $chh_ot =~ s/^/CHH_OT_/;
753 $chh_ot =~ s/sam$/txt/;
754 $chh_ot =~ s/$/.txt/ unless ($chh_ot =~ /\.txt$/);
755 $chh_ot = $output_dir . $chh_ot;
756 push @sorting_files,$chh_ot;
757 open ($fhs{0}->{CHH},'>',$chh_ot) or die "Failed to write to $chh_ot $!\n";
758 print "Writing result file containing methylation information for C in CHH context from the original top strand to $chh_ot\n";
759
760 unless($no_header){
761 print {$fhs{0}->{CHH}} "Bismark methylation extractor version $version\n";
762 }
763
764 $chh_ctot =~ s/^/CHH_CTOT_/;
765 $chh_ctot =~ s/sam$/txt/;
766 $chh_ctot =~ s/$/.txt/ unless ($chh_ctot =~ /\.txt$/);
767 $chh_ctot = $output_dir . $chh_ctot;
768 push @sorting_files,$chh_ctot;
769 open ($fhs{1}->{CHH},'>',$chh_ctot) or die "Failed to write to $chh_ctot $!\n";
770 print "Writing result file containing methylation information for C in CHH context from the complementary to original top strand to $chh_ctot\n";
771
772 unless($no_header){
773 print {$fhs{1}->{CHH}} "Bismark methylation extractor version $version\n";
774 }
775
776 $chh_ctob =~ s/^/CHH_CTOB_/;
777 $chh_ctob =~ s/sam$/txt/;
778 $chh_ctob =~ s/$/.txt/ unless ($chh_ctob =~ /\.txt$/);
779 $chh_ctob = $output_dir . $chh_ctob;
780 push @sorting_files,$chh_ctob;
781 open ($fhs{2}->{CHH},'>',$chh_ctob) or die "Failed to write to $chh_ctob $!\n";
782 print "Writing result file containing methylation information for C in CHH context from the complementary to original bottom strand to $chh_ctob\n";
783
784 unless($no_header){
785 print {$fhs{2}->{CHH}} "Bismark methylation extractor version $version\n";
786 }
787
788 $chh_ob =~ s/^/CHH_OB_/;
789 $chh_ob =~ s/sam$/txt/;
790 $chh_ob =~ s/$/.txt/ unless ($chh_ob =~ /\.txt$/);
791 $chh_ob = $output_dir . $chh_ob;
792 push @sorting_files,$chh_ob;
793 open ($fhs{3}->{CHH},'>',$chh_ob) or die "Failed to write to $chh_ob $!\n";
794 print "Writing result file containing methylation information for C in CHH context from the original bottom strand to $chh_ob\n\n";
795
796 unless($no_header){
797 print {$fhs{3}->{CHH}} "Bismark methylation extractor version $version\n";
798 }
799 }
800
801 my $methylation_call_strings_processed = 0;
802 my $line_count = 0;
803
804 ### proceeding differently now for single-end or paired-end Bismark files
805
806 ### PROCESSING SINGLE-END RESULT FILES
807 if ($single) {
808
809 ### also proceeding differently now for SAM format or vanilla Bismark format files
810 if ($vanilla) { # old vanilla Bismark output format
811 while (<IN>) {
812 ++$line_count;
813 warn "Processed lines: $line_count\n" if ($line_count%500000==0);
814
815 ### $seq here is the chromosomal sequence (to use for the repeat analysis for example)
816 my ($id,$strand,$chrom,$start,$seq,$meth_call,$read_conversion,$genome_conversion) = (split("\t"))[0,1,2,3,6,7,8,9];
817
818 ### 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
819 ### last position
820 chomp $genome_conversion;
821
822 my $index;
823 if ($meth_call) {
824
825 if ($read_conversion eq 'CT' and $genome_conversion eq 'CT') { ## original top strand
826 $index = 0;
827 } elsif ($read_conversion eq 'GA' and $genome_conversion eq 'CT') { ## complementary to original top strand
828 $index = 1;
829 } elsif ($read_conversion eq 'CT' and $genome_conversion eq 'GA') { ## original bottom strand
830 $index = 3;
831 } elsif ($read_conversion eq 'GA' and $genome_conversion eq 'GA') { ## complementary to original bottom strand
832 $index = 2;
833 } else {
834 die "Unexpected combination of read and genome conversion: '$read_conversion' / '$genome_conversion'\n";
835 }
836
837 ### Clipping off the first <int> number of bases from the methylation call string as specified with --ignore <int>
838 if ($ignore) {
839 $meth_call = substr($meth_call,$ignore,length($meth_call)-$ignore);
840
841 ### If we are clipping off some bases at the start we need to adjust the start position of the alignments accordingly!
842 if ($strand eq '+') {
843 $start += $ignore;
844 } elsif ($strand eq '-') {
845 $start += length($meth_call)-1; ## $meth_call is already shortened!
846 } else {
847 die "Alignment did not have proper strand information: $strand\n";
848 }
849 }
850 ### printing out the methylation state of every C in the read
851 print_individual_C_methylation_states_single_end($meth_call,$chrom,$start,$id,$strand,$index);
852
853 ++$methylation_call_strings_processed; # 1 per single-end result
854 }
855 }
856 } else { # processing single-end SAM format (default)
857 while (<IN>) {
858 ### SAM format can either start with header lines (starting with @) or start with alignments directly
859 if (/^\@/) { # skipping header lines (starting with @)
860 warn "skipping SAM header line:\t$_";
861 next;
862 }
863
864 ++$line_count;
865 warn "Processed lines: $line_count\n" if ($line_count%500000==0);
866
867 # example read in SAM format
868 # 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
869 ###
870
871 # < 0.7.6 my ($id,$chrom,$start,$meth_call,$read_conversion,$genome_conversion) = (split("\t"))[0,2,3,13,14,15];
872 # < 0.7.6 $meth_call =~ s/^XM:Z://;
873 # < 0.7.6 $read_conversion =~ s/^XR:Z://;
874 # < 0.7.6 $genome_conversion =~ s/^XG:Z://;
875
876 my ($id,$chrom,$start,$cigar) = (split("\t"))[0,2,3,5];
877
878 ### detecting the following SAM flags in case the SAM entry was shuffled by CRAM or Goby compression/decompression
879 my $meth_call; ### Thanks to Zachary Zeno for this solution
880 my $read_conversion;
881 my $genome_conversion;
882
883 while ( /(XM|XR|XG):Z:([^\t]+)/g ) {
884 my $tag = $1;
885 my $value = $2;
886
887 if ($tag eq "XM") {
888 $meth_call = $value;
889 $meth_call =~ s/\r//;
890 } elsif ($tag eq "XR") {
891 $read_conversion = $value;
892 $read_conversion =~ s/\r//;
893 } elsif ($tag eq "XG") {
894 $genome_conversion = $value;
895 $genome_conversion =~ s/\r//;
896 }
897 }
898
899 my $strand;
900 chomp $genome_conversion;
901 # print "$meth_call\n$read_conversion\n$genome_conversion\n";
902
903 my $index;
904 if ($meth_call) {
905 if ($read_conversion eq 'CT' and $genome_conversion eq 'CT') { ## original top strand
906 $index = 0;
907 $strand = '+';
908 } elsif ($read_conversion eq 'GA' and $genome_conversion eq 'CT') { ## complementary to original top strand
909 $index = 1;
910 $strand = '-';
911 } elsif ($read_conversion eq 'GA' and $genome_conversion eq 'GA') { ## complementary to original bottom strand
912 $index = 2;
913 $strand = '+';
914 } elsif ($read_conversion eq 'CT' and $genome_conversion eq 'GA') { ## original bottom strand
915 $index = 3;
916 $strand = '-';
917 } else {
918 die "Unexpected combination of read and genome conversion: '$read_conversion' / '$genome_conversion'\n";
919 }
920
921 ### If the read is in SAM format we need to reverse the methylation call if the read has been reverse-complemented for the output
922 if ($strand eq '-') {
923 $meth_call = reverse $meth_call;
924 }
925
926 ### Clipping off the first <int> number of bases from the methylation call string as specified with --ignore <int>
927 if ($ignore) {
928 # print "\n\n$meth_call\n";
929 $meth_call = substr($meth_call,$ignore,length($meth_call)-$ignore);
930 # print "$meth_call\n";
931 ### If we are ignoring a part of the sequence we also need to adjust the cigar string accordingly
932
933 my @len = split (/\D+/,$cigar); # storing the length per operation
934 my @ops = split (/\d+/,$cigar); # storing the operation
935 shift @ops; # remove the empty first element
936 die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len == scalar @ops);
937
938 my @comp_cigar; # building an array with all CIGAR operations
939 foreach my $index (0..$#len) {
940 foreach (1..$len[$index]) {
941 # print "$ops[$index]";
942 push @comp_cigar, $ops[$index];
943 }
944 }
945 # print "original CIGAR: $cigar\n";
946 # print "original CIGAR: @comp_cigar\n";
947
948 ### If we are clipping off some bases at the start we need to adjust the start position of the alignments accordingly!
949 if ($strand eq '+') {
950
951 my $D_count = 0; # counting all deletions that affect the ignored genomic position, i.e. Deletions and insertions
952 my $I_count = 0;
953
954 for (1..$ignore) {
955 my $op = shift @comp_cigar; # adjusting composite CIGAR string by removing $ignore operations from the start
956 # print "$_ deleted $op\n";
957
958 while ($op eq 'D') { # repeating this for deletions (D)
959 $D_count++;
960 $op = shift @comp_cigar;
961 # print "$_ deleted $op\n";
962 }
963 if ($op eq 'I') { # adjusting the genomic position for insertions (I)
964 $I_count++;
965 }
966 }
967 $start += $ignore + $D_count - $I_count;
968 # print "start $start\t ignore: $ignore\t D count: $D_count I_count: $I_count\n";
969 } elsif ($strand eq '-') {
970
971 for (1..$ignore) {
972 my $op = pop @comp_cigar; # adjusting composite CIGAR string by removing $ignore operations, here the last value of the array
973 while ($op eq 'D') { # repeating this for deletions (D)
974 $op = pop @comp_cigar;
975 }
976 }
977
978 ### For reverse strand alignments we need to determine the number of matching bases (M) or deletions (D) in the read from the CIGAR
979 ### string to be able to work out the starting position of the read which is on the 3' end of the sequence
980 my $MD_count = 0; # counting all operations that affect the genomic position, i.e. M and D. Insertions do not affect the start position
981 foreach (@comp_cigar) {
982 ++$MD_count if ($_ eq 'M' or $_ eq 'D');
983 }
984 $start += $MD_count - 1;
985 }
986
987 ### reconstituting shortened CIGAR string
988 my $new_cigar;
989 my $count = 0;
990 my $last_op;
991 # print "ignore adjusted: @comp_cigar\n";
992 foreach my $op (@comp_cigar) {
993 unless (defined $last_op){
994 $last_op = $op;
995 ++$count;
996 next;
997 }
998 if ($last_op eq $op) {
999 ++$count;
1000 } else {
1001 $new_cigar .= "$count$last_op";
1002 $last_op = $op;
1003 $count = 1;
1004 }
1005 }
1006 $new_cigar .= "$count$last_op"; # appending the last operation and count
1007 $cigar = $new_cigar;
1008 # print "ignore adjusted scalar: $cigar\n";
1009 }
1010 }
1011 ### printing out the methylation state of every C in the read
1012 print_individual_C_methylation_states_single_end($meth_call,$chrom,$start,$id,$strand,$index,$cigar);
1013
1014 ++$methylation_call_strings_processed; # 1 per single-end result
1015 }
1016 }
1017 }
1018
1019 ### PROCESSING PAIRED-END RESULT FILES
1020 elsif ($paired) {
1021
1022 ### proceeding differently now for SAM format or vanilla Bismark format files
1023 if ($vanilla) { # old vanilla Bismark paired-end output format
1024 while (<IN>) {
1025 ++$line_count;
1026 warn "processed line: $line_count\n" if ($line_count%500000==0);
1027
1028 ### $seq here is the chromosomal sequence (to use for the repeat analysis for example)
1029 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];
1030
1031 my $index;
1032 chomp $genome_conversion;
1033
1034 if ($first_read_conversion eq 'CT' and $genome_conversion eq 'CT') {
1035 $index = 0; ## this is OT
1036 } elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'GA') {
1037 $index = 2; ## this is CTOB!!!
1038 } elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'CT') {
1039 $index = 1; ## this is CTOT!!!
1040 } elsif ($first_read_conversion eq 'CT' and $genome_conversion eq 'GA') {
1041 $index = 3; ## this is OB
1042 } else {
1043 die "Unexpected combination of read and genome conversion: $first_read_conversion / $genome_conversion\n";
1044 }
1045
1046 if ($meth_call_1 and $meth_call_2) {
1047 ### Clipping off the first <int> number of bases from the methylation call strings as specified with '--ignore <int>'
1048 if ($ignore) {
1049 $meth_call_1 = substr($meth_call_1,$ignore,length($meth_call_1)-$ignore);
1050 $meth_call_2 = substr($meth_call_2,$ignore,length($meth_call_2)-$ignore);
1051
1052 ### we also need to adjust the start and end positions of the alignments accordingly if '--ignore' was specified
1053 $start_read_1 += $ignore;
1054 $end_read_2 -= $ignore;
1055 }
1056 my $end_read_1;
1057 my $start_read_2;
1058
1059 if ($strand eq '+') {
1060
1061 $end_read_1 = $start_read_1+length($meth_call_1)-1;
1062 $start_read_2 = $end_read_2-length($meth_call_2)+1;
1063
1064 ## we first pass the first read which is in + orientation on the forward strand
1065 print_individual_C_methylation_states_paired_end_files($meth_call_1,$chrom,$start_read_1,$id,'+',$index,0,0);
1066
1067 # we next pass the second read which is in - orientation on the reverse strand
1068 ### 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
1069 print_individual_C_methylation_states_paired_end_files($meth_call_2,$chrom,$end_read_2,$id,'-',$index,$no_overlap,$end_read_1);
1070 } else {
1071
1072 $end_read_1 = $start_read_1+length($meth_call_2)-1; # read 1 is the second reported read!
1073 $start_read_2 = $end_read_2-length($meth_call_1)+1; # read 2 is the first reported read!
1074
1075 ## we first pass the first read which is in - orientation on the reverse strand
1076 print_individual_C_methylation_states_paired_end_files($meth_call_1,$chrom,$end_read_2,$id,'-',$index,0,0);
1077
1078 # we next pass the second read which is in + orientation on the forward strand
1079 ### 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
1080 print_individual_C_methylation_states_paired_end_files($meth_call_2,$chrom,$start_read_1,$id,'+',$index,$no_overlap,$start_read_2);
1081 }
1082
1083 $methylation_call_strings_processed += 2; # paired-end = 2 methylation calls
1084 }
1085 }
1086 } else { # Bismark paired-end SAM output format (default)
1087 while (<IN>) {
1088 ### SAM format can either start with header lines (starting with @) or start with alignments directly
1089 if (/^\@/) { # skipping header lines (starting with @)
1090 warn "skipping SAM header line:\t$_";
1091 next;
1092 }
1093
1094 ++$line_count;
1095 warn "Processed lines: $line_count\n" if ($line_count%500000==0);
1096
1097 # example paired-end reads in SAM format (2 consecutive lines)
1098 # 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
1099 # 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
1100
1101 # < 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];
1102
1103 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
1104 my $meth_call_1;
1105 my $first_read_conversion;
1106 my $genome_conversion;
1107
1108 while ( /(XM|XR|XG):Z:([^\t]+)/g ) {
1109 my $tag = $1;
1110 my $value = $2;
1111
1112 if ($tag eq "XM") {
1113 $meth_call_1 = $value;
1114 $meth_call_1 =~ s/\r//;
1115 } elsif ($tag eq "XR") {
1116 $first_read_conversion = $value;
1117 $first_read_conversion =~ s/\r//;
1118 } elsif ($tag eq "XG") {
1119 $genome_conversion = $value;
1120 $genome_conversion =~ s/\r//;
1121 }
1122 }
1123
1124 $_ = <IN>; # reading in the paired read
1125
1126 # < version 0.7.6 my ($id_2,$start_read_2,$meth_call_2,$second_read_conversion) = (split("\t"))[0,3,13,14];
1127 # < version 0.7.6 $meth_call_1 =~ s/^XM:Z://;
1128 # < version 0.7.6 $meth_call_2 =~ s/^XM:Z://;
1129 # < version 0.7.6 $first_read_conversion =~ s/^XR:Z://;
1130 # < version 0.7.6 $second_read_conversion =~ s/^XR:Z://;
1131
1132 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
1133
1134 my $meth_call_2;
1135 my $second_read_conversion;
1136
1137 while ( /(XM|XR):Z:([^\t]+)/g ) {
1138 my $tag = $1;
1139 my $value = $2;
1140
1141 if ($tag eq "XM") {
1142 $meth_call_2 = $value;
1143 $meth_call_2 =~ s/\r//;
1144 } elsif ($tag eq "XR") {
1145 $second_read_conversion = $value;
1146 $second_read_conversion = s/\r//;
1147 }
1148 }
1149
1150 # < version 0.7.6 $genome_conversion =~ s/^XG:Z://;
1151 chomp $genome_conversion; # in case it captured a new line character
1152
1153 # print join ("\t",$meth_call_1,$meth_call_2,$first_read_conversion,$second_read_conversion,$genome_conversion),"\n";
1154
1155 my $index;
1156 my $strand;
1157
1158 if ($first_read_conversion eq 'CT' and $genome_conversion eq 'CT') {
1159 $index = 0; ## this is OT
1160 $strand = '+';
1161 } elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'CT') {
1162 $index = 1; ## this is CTOT
1163 $strand = '-';
1164 } elsif ($first_read_conversion eq 'GA' and $genome_conversion eq 'GA') {
1165 $index = 2; ## this is CTOB
1166 $strand = '+';
1167 } elsif ($first_read_conversion eq 'CT' and $genome_conversion eq 'GA') {
1168 $index = 3; ## this is OB
1169 $strand = '-';
1170 } else {
1171 die "Unexpected combination of read and genome conversion: $first_read_conversion / $genome_conversion\n";
1172 }
1173
1174 ### reversing the methylation call of the read that was reverse-complemented
1175 if ($strand eq '+') {
1176 $meth_call_2 = reverse $meth_call_2;
1177 } else {
1178 $meth_call_1 = reverse $meth_call_1;
1179 }
1180
1181 if ($meth_call_1 and $meth_call_2) {
1182
1183 my $end_read_1;
1184
1185 ### READ 1
1186 my @len_1 = split (/\D+/,$cigar_1); # storing the length per operation
1187 my @ops_1 = split (/\d+/,$cigar_1); # storing the operation
1188 shift @ops_1; # remove the empty first element
1189 die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len_1 == scalar @ops_1);
1190
1191 my @comp_cigar_1; # building an array with all CIGAR operations
1192 foreach my $index (0..$#len_1) {
1193 foreach (1..$len_1[$index]) {
1194 # print "$ops_1[$index]";
1195 push @comp_cigar_1, $ops_1[$index];
1196 }
1197 }
1198 # print "original CIGAR read 1: $cigar_1\n";
1199 # print "original CIGAR read 1: @comp_cigar_1\n";
1200
1201 ### READ 2
1202 my @len_2 = split (/\D+/,$cigar_2); # storing the length per operation
1203 my @ops_2 = split (/\d+/,$cigar_2); # storing the operation
1204 shift @ops_2; # remove the empty first element
1205 die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len_2 == scalar @ops_2);
1206 my @comp_cigar_2; # building an array with all CIGAR operations for read 2
1207 foreach my $index (0..$#len_2) {
1208 foreach (1..$len_2[$index]) {
1209 # print "$ops_2[$index]";
1210 push @comp_cigar_2, $ops_2[$index];
1211 }
1212 }
1213 # print "original CIGAR read 2: $cigar_2\n";
1214 # print "original CIGAR read 2: @comp_cigar_2\n";
1215
1216 if ($ignore) {
1217 ### Clipping off the first <int> number of bases from the methylation call strings as specified with '--ignore <int>'
1218 ### the methylation calls have already been reversed where necessary
1219 $meth_call_1 = substr($meth_call_1,$ignore,length($meth_call_1)-$ignore);
1220 $meth_call_2 = substr($meth_call_2,$ignore,length($meth_call_2)-$ignore);
1221
1222 ### If we are ignoring a part of the sequence we also need to adjust the cigar string accordingly
1223
1224 if ($strand eq '+') {
1225
1226 ### if the (read 1) strand information is '+', read 1 needs to be trimmed from the start
1227 my $D_count_1 = 0; # counting all deletions that affect the ignored genomic position for read 1, i.e. Deletions and insertions
1228 my $I_count_1 = 0;
1229
1230 for (1..$ignore) {
1231 my $op = shift @comp_cigar_1; # adjusting composite CIGAR string of read 1 by removing $ignore operations from the start
1232 # print "$_ deleted $op\n";
1233
1234 while ($op eq 'D') { # repeating this for deletions (D)
1235 $D_count_1++;
1236 $op = shift @comp_cigar_1;
1237 # print "$_ deleted $op\n";
1238 }
1239 if ($op eq 'I') { # adjusting the genomic position for insertions (I)
1240 $I_count_1++;
1241 }
1242 }
1243
1244 $start_read_1 += $ignore + $D_count_1 - $I_count_1;
1245 # print "start read 1 $start_read_1\t ignore: $ignore\t D count 1: $D_count_1\tI_count 1: $I_count_1\n";
1246
1247 ### if the (read 1) strand information is '+', read 2 needs to be trimmed from the back
1248
1249 for (1..$ignore) {
1250 my $op = pop @comp_cigar_2; # adjusting composite CIGAR string by removing $ignore operations, here the last value of the array
1251 while ($op eq 'D') { # repeating this for deletions (D)
1252 $op = pop @comp_cigar_2;
1253 }
1254 }
1255 # the start position of reads mapping to the reverse strand is being adjusted further below
1256 } elsif ($strand eq '-') {
1257
1258 ### if the (read 1) strand information is '-', read 1 needs to be trimmed from the back
1259 for (1..$ignore) {
1260 my $op = pop @comp_cigar_1; # adjusting composite CIGAR string by removing $ignore operations, here the last value of the array
1261 while ($op eq 'D') { # repeating this for deletions (D)
1262 $op = pop @comp_cigar_1;
1263 }
1264 }
1265 # the start position of reads mapping to the reverse strand is being adjusted further below
1266
1267 ### if the (read 1) strand information is '-', read 2 needs to be trimmed from the start
1268 my $D_count_2 = 0; # counting all deletions that affect the ignored genomic position for read 2, i.e. Deletions and insertions
1269 my $I_count_2 = 0;
1270
1271 for (1..$ignore) {
1272 my $op = shift @comp_cigar_2; # adjusting composite CIGAR string of read 2 by removing $ignore operations from the start
1273 # print "$_ deleted $op\n";
1274
1275 while ($op eq 'D') { # repeating this for deletions (D)
1276 $D_count_2++;
1277 $op = shift @comp_cigar_2;
1278 # print "$_ deleted $op\n";
1279 }
1280 if ($op eq 'I') { # adjusting the genomic position for insertions (I)
1281 $I_count_2++;
1282 }
1283 }
1284
1285 $start_read_2 += $ignore + $D_count_2 - $I_count_2;
1286 # print "start read 2 $start_read_2\t ignore: $ignore\t D count 2: $D_count_2\tI_count 2: $I_count_2\n";
1287
1288 }
1289
1290 ### reconstituting shortened CIGAR string 1
1291 my $new_cigar_1;
1292 my $count_1 = 0;
1293 my $last_op_1;
1294 # print "ignore adjusted CIGAR 1: @comp_cigar_1\n";
1295 foreach my $op (@comp_cigar_1) {
1296 unless (defined $last_op_1){
1297 $last_op_1 = $op;
1298 ++$count_1;
1299 next;
1300 }
1301 if ($last_op_1 eq $op) {
1302 ++$count_1;
1303 } else {
1304 $new_cigar_1 .= "$count_1$last_op_1";
1305 $last_op_1 = $op;
1306 $count_1 = 1;
1307 }
1308 }
1309 $new_cigar_1 .= "$count_1$last_op_1"; # appending the last operation and count
1310 $cigar_1 = $new_cigar_1;
1311 # print "ignore adjusted CIGAR 1 scalar: $cigar_1\n";
1312
1313 ### reconstituting shortened CIGAR string 2
1314 my $new_cigar_2;
1315 my $count_2 = 0;
1316 my $last_op_2;
1317 # print "ignore adjusted CIGAR 2: @comp_cigar_2\n";
1318 foreach my $op (@comp_cigar_2) {
1319 unless (defined $last_op_2){
1320 $last_op_2 = $op;
1321 ++$count_2;
1322 next;
1323 }
1324 if ($last_op_2 eq $op) {
1325 ++$count_2;
1326 } else {
1327 $new_cigar_2 .= "$count_2$last_op_2";
1328 $last_op_2 = $op;
1329 $count_2 = 1;
1330 }
1331 }
1332 $new_cigar_2 .= "$count_2$last_op_2"; # appending the last operation and count
1333 $cigar_2 = $new_cigar_2;
1334 # print "ignore adjusted CIGAR 2 scalar: $cigar_2\n";
1335
1336 }
1337
1338 if ($strand eq '+') {
1339 ### adjusting the start position for all reads mapping to the reverse strand, in this case read 2
1340 @comp_cigar_2 = reverse@comp_cigar_2; # the CIGAR string needs to be reversed for all reads aligning to the reverse strand, too
1341 # print "reverse: @comp_cigar_2\n";
1342
1343 my $MD_count_1 = 0;
1344 foreach (@comp_cigar_1) {
1345 ++$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
1346 }
1347
1348 my $MD_count_2 = 0;
1349 foreach (@comp_cigar_2) {
1350 ++$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
1351 }
1352
1353 $end_read_1 = $start_read_1 + $MD_count_1 - 1;
1354 $start_read_2 += $MD_count_2 - 1; ## Passing on the start position on the reverse strand
1355 } else {
1356 ### adjusting the start position for all reads mapping to the reverse strand, in this case read 1
1357
1358 @comp_cigar_1 = reverse@comp_cigar_1; # the CIGAR string needs to be reversed for all reads aligning to the reverse strand, too
1359 # print "reverse: @comp_cigar_1\n";
1360
1361 my $MD_count_1 = 0;
1362 foreach (@comp_cigar_1) {
1363 ++$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
1364 }
1365
1366 $end_read_1 = $start_read_1;
1367 $start_read_1 += $MD_count_1 - 1; ### Passing on the start position on the reverse strand
1368
1369 }
1370
1371 if ($strand eq '+') {
1372 ## we first pass the first read which is in + orientation on the forward strand
1373 print_individual_C_methylation_states_paired_end_files($meth_call_1,$chrom,$start_read_1,$id_1,'+',$index,0,0,$cigar_1);
1374
1375 # we next pass the second read which is in - orientation on the reverse strand
1376 ### 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
1377 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);
1378 } else {
1379 ## we first pass the first read which is in - orientation on the reverse strand
1380 print_individual_C_methylation_states_paired_end_files($meth_call_1,$chrom,$start_read_1,$id_1,'-',$index,0,0,$cigar_1);
1381
1382 # we next pass the second read which is in + orientation on the forward strand
1383 ### 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
1384 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);
1385 }
1386
1387 $methylation_call_strings_processed += 2; # paired-end = 2 methylation calls
1388 }
1389 }
1390 }
1391 } else {
1392 die "Single-end or paired-end reads not specified properly\n";
1393 }
1394
1395 print "\n\nProcessed $line_count lines from $filename in total\n";
1396 print "Total number of methylation call strings processed: $methylation_call_strings_processed\n\n";
1397 if ($report) {
1398 print REPORT "Total number of methylation call strings processed: $methylation_call_strings_processed\n\n";
1399 }
1400 print_splitting_report ();
1401 }
1402
1403
1404
1405 sub print_splitting_report{
1406
1407 ### Calculating methylation percentages if applicable
1408
1409 my $percent_meCpG;
1410 if (($counting{total_meCpG_count}+$counting{total_unmethylated_CpG_count}) > 0){
1411 $percent_meCpG = sprintf("%.1f",100*$counting{total_meCpG_count}/($counting{total_meCpG_count}+$counting{total_unmethylated_CpG_count}));
1412 }
1413
1414 my $percent_meCHG;
1415 if (($counting{total_meCHG_count}+$counting{total_unmethylated_CHG_count}) > 0){
1416 $percent_meCHG = sprintf("%.1f",100*$counting{total_meCHG_count}/($counting{total_meCHG_count}+$counting{total_unmethylated_CHG_count}));
1417 }
1418
1419 my $percent_meCHH;
1420 if (($counting{total_meCHH_count}+$counting{total_unmethylated_CHH_count}) > 0){
1421 $percent_meCHH = sprintf("%.1f",100*$counting{total_meCHH_count}/($counting{total_meCHH_count}+$counting{total_unmethylated_CHH_count}));
1422 }
1423
1424 my $percent_non_CpG_methylation;
1425 if ($merge_non_CpG){
1426 if ( ($counting{total_meCHH_count}+$counting{total_unmethylated_CHH_count}+$counting{total_meCHG_count}+$counting{total_unmethylated_CHG_count}) > 0){
1427 $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} ) );
1428 }
1429 }
1430
1431 if ($report){
1432 ### detailed information about Cs analysed
1433 print REPORT "Final Cytosine Methylation Report\n",'='x33,"\n";
1434
1435 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};
1436 print REPORT "Total number of C's analysed:\t$total_number_of_C\n\n";
1437
1438 print REPORT "Total methylated C's in CpG context:\t$counting{total_meCpG_count}\n";
1439 print REPORT "Total methylated C's in CHG context:\t$counting{total_meCHG_count}\n";
1440 print REPORT "Total methylated C's in CHH context:\t$counting{total_meCHH_count}\n\n";
1441
1442 print REPORT "Total C to T conversions in CpG context:\t$counting{total_unmethylated_CpG_count}\n";
1443 print REPORT "Total C to T conversions in CHG context:\t$counting{total_unmethylated_CHG_count}\n";
1444 print REPORT "Total C to T conversions in CHH context:\t$counting{total_unmethylated_CHH_count}\n\n";
1445
1446 ### calculating methylated CpG percentage if applicable
1447 if ($percent_meCpG){
1448 print REPORT "C methylated in CpG context:\t${percent_meCpG}%\n";
1449 }
1450 else{
1451 print REPORT "Can't determine percentage of methylated Cs in CpG context if value was 0\n";
1452 }
1453
1454 ### 2-Context Output
1455 if ($merge_non_CpG){
1456 if ($percent_non_CpG_methylation){
1457 print REPORT "C methylated in non-CpG context:\t${percent_non_CpG_methylation}%\n\n\n";
1458 }
1459 else{
1460 print REPORT "Can't determine percentage of methylated Cs in non-CpG context if value was 0\n\n\n";
1461 }
1462 }
1463
1464 ### 3 Context Output
1465 else{
1466 ### calculating methylated CHG percentage if applicable
1467 if ($percent_meCHG){
1468 print REPORT "C methylated in CHG context:\t${percent_meCHG}%\n";
1469 }
1470 else{
1471 print REPORT "Can't determine percentage of methylated Cs in CHG context if value was 0\n";
1472 }
1473
1474 ### calculating methylated CHH percentage if applicable
1475 if ($percent_meCHH){
1476 print REPORT "C methylated in CHH context:\t${percent_meCHH}%\n\n\n";
1477 }
1478 else{
1479 print REPORT "Can't determine percentage of methylated Cs in CHH context if value was 0\n\n\n";
1480 }
1481 }
1482 }
1483
1484 ### detailed information about Cs analysed for on-screen report
1485 print "Final Cytosine Methylation Report\n",'='x33,"\n";
1486
1487 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};
1488 print "Total number of C's analysed:\t$total_number_of_C\n\n";
1489
1490 print "Total methylated C's in CpG context:\t$counting{total_meCpG_count}\n";
1491 print "Total methylated C's in CHG context:\t$counting{total_meCHG_count}\n";
1492 print "Total methylated C's in CHH context:\t$counting{total_meCHH_count}\n\n";
1493
1494 print "Total C to T conversions in CpG context:\t$counting{total_unmethylated_CpG_count}\n";
1495 print "Total C to T conversions in CHG context:\t$counting{total_unmethylated_CHG_count}\n";
1496 print "Total C to T conversions in CHH context:\t$counting{total_unmethylated_CHH_count}\n\n";
1497
1498 ### printing methylated CpG percentage if applicable
1499 if ($percent_meCpG){
1500 print "C methylated in CpG context:\t${percent_meCpG}%\n";
1501 }
1502 else{
1503 print "Can't determine percentage of methylated Cs in CpG context if value was 0\n";
1504 }
1505
1506 ### 2-Context Output
1507 if ($merge_non_CpG){
1508 if ($percent_non_CpG_methylation){
1509 print "C methylated in non-CpG context:\t${percent_non_CpG_methylation}%\n\n\n";
1510 }
1511 else{
1512 print "Can't determine percentage of methylated Cs in non-CpG context if value was 0\n\n\n";
1513 }
1514 }
1515
1516 ### 3-Context Output
1517 else{
1518 ### printing methylated CHG percentage if applicable
1519 if ($percent_meCHG){
1520 print "C methylated in CHG context:\t${percent_meCHG}%\n";
1521 }
1522 else{
1523 print "Can't determine percentage of methylated Cs in CHG context if value was 0\n";
1524 }
1525
1526 ### printing methylated CHH percentage if applicable
1527 if ($percent_meCHH){
1528 print "C methylated in CHH context:\t${percent_meCHH}%\n\n\n";
1529 }
1530 else{
1531 print "Can't determine percentage of methylated Cs in CHH context if value was 0\n\n\n";
1532 }
1533 }
1534 }
1535
1536
1537
1538
1539
1540 sub print_individual_C_methylation_states_paired_end_files{
1541
1542 my ($meth_call,$chrom,$start,$id,$strand,$filehandle_index,$no_overlap,$end_read_1,$cigar) = @_;
1543 my @methylation_calls = split(//,$meth_call);
1544
1545 #################################################################
1546 ### . for bases not involving cytosines ###
1547 ### X for methylated C in CHG context (was protected) ###
1548 ### x for not methylated C in CHG context (was converted) ###
1549 ### H for methylated C in CHH context (was protected) ###
1550 ### h for not methylated C in CHH context (was converted) ###
1551 ### Z for methylated C in CpG context (was protected) ###
1552 ### z for not methylated C in CpG context (was converted) ###
1553 #################################################################
1554
1555 my $methyl_CHG_count = 0;
1556 my $methyl_CHH_count = 0;
1557 my $methyl_CpG_count = 0;
1558 my $unmethylated_CHG_count = 0;
1559 my $unmethylated_CHH_count = 0;
1560 my $unmethylated_CpG_count = 0;
1561
1562 my @len;
1563 my @ops;
1564 my $pos_offset = 0; # this is only relevant for SAM reads with insertions or deletions
1565 my $cigar_offset = 0; # again, this is only relevant for SAM reads containing indels
1566 my @comp_cigar;
1567
1568 if ($cigar){ # parsing CIGAR string
1569 @len = split (/\D+/,$cigar); # storing the length per operation
1570 @ops = split (/\d+/,$cigar); # storing the operation
1571 shift @ops; # remove the empty first element
1572 die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len == scalar @ops);
1573
1574 foreach my $index (0..$#len){
1575 foreach (1..$len[$index]){
1576 # print "$ops[$index]";
1577 push @comp_cigar, $ops[$index];
1578 }
1579 }
1580 # warn "\nDetected CIGAR string: $cigar\n";
1581 # warn "Length of methylation call: ",length $meth_call,"\n";
1582 # warn "number of operations: ",scalar @ops,"\n";
1583 # warn "number of length digits: ",scalar @len,"\n\n";
1584 # print @comp_cigar,"\n";
1585 # print "$meth_call\n\n";
1586 # sleep (1);
1587 }
1588
1589
1590 if ($strand eq '-') {
1591
1592 ### the CIGAR string needs to be reversed, the methylation call has already been reversed above
1593 @comp_cigar = reverse@comp_cigar; # the CIGAR string needs to be reversed for all reads aligning to the reverse strand, too
1594 # print "reverse CIGAR string: @comp_cigar\n";
1595
1596 ### the start position of paired-end files has already been corrected, see above
1597 }
1598
1599 ### THIS IS AN OPTIONAL 2-CONTEXT (CpG and non-CpG) SECTION IF --merge_non_CpG was specified
1600
1601 if ($merge_non_CpG) {
1602
1603 if ($no_overlap) {
1604
1605 ### single-file CpG and non-CpG context output
1606 if ($full) {
1607 if ($strand eq '+') {
1608 for my $index (0..$#methylation_calls) {
1609
1610 if ($cigar){ # only needed for SAM files
1611 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
1612 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
1613 $cigar_offset += $cigar_mod;
1614 $pos_offset += $pos_mod;
1615 }
1616
1617 ### Returning as soon as the methylation calls start overlapping
1618 if ($start+$index+$pos_offset >= $end_read_1) {
1619 return;
1620 }
1621
1622 if ($methylation_calls[$index] eq 'X') {
1623 $counting{total_meCHG_count}++;
1624 print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1625 } elsif ($methylation_calls[$index] eq 'x') {
1626 $counting{total_unmethylated_CHG_count}++;
1627 print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1628 } elsif ($methylation_calls[$index] eq 'Z') {
1629 $counting{total_meCpG_count}++;
1630 print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1631 } elsif ($methylation_calls[$index] eq 'z') {
1632 $counting{total_unmethylated_CpG_count}++;
1633 print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1634 } elsif ($methylation_calls[$index] eq 'H') {
1635 $counting{total_meCHH_count}++;
1636 print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1637 } elsif ($methylation_calls[$index] eq 'h') {
1638 $counting{total_unmethylated_CHH_count}++;
1639 print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1640 }
1641 elsif ($methylation_calls[$index] eq '.'){}
1642 else{
1643 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
1644 }
1645 }
1646 } elsif ($strand eq '-') {
1647 for my $index (0..$#methylation_calls) {
1648
1649 if ($cigar){ # only needed for SAM files
1650 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
1651 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
1652 $cigar_offset += $cigar_mod;
1653 $pos_offset += $pos_mod;
1654 }
1655
1656 ### Returning as soon as the methylation calls start overlapping
1657 if ($start-$index+$pos_offset <= $end_read_1) {
1658 return;
1659 }
1660
1661 if ($methylation_calls[$index] eq 'X') {
1662 $counting{total_meCHG_count}++;
1663 print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1664 } elsif ($methylation_calls[$index] eq 'x') {
1665 $counting{total_unmethylated_CHG_count}++;
1666 print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1667 } elsif ($methylation_calls[$index] eq 'Z') {
1668 $counting{total_meCpG_count}++;
1669 print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1670 } elsif ($methylation_calls[$index] eq 'z') {
1671 $counting{total_unmethylated_CpG_count}++;
1672 print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1673 } elsif ($methylation_calls[$index] eq 'H') {
1674 $counting{total_meCHH_count}++;
1675 print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1676 } elsif ($methylation_calls[$index] eq 'h') {
1677 $counting{total_unmethylated_CHH_count}++;
1678 print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1679 }
1680 elsif ($methylation_calls[$index] eq '.') {}
1681 else{
1682 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
1683 }
1684 }
1685 } else {
1686 die "The read orientation was neither + nor -: '$strand'\n";
1687 }
1688 }
1689
1690 ### strand-specific methylation output
1691 else {
1692 if ($strand eq '+') {
1693 for my $index (0..$#methylation_calls) {
1694
1695 if ($cigar){ # only needed for SAM files
1696 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
1697 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
1698 $cigar_offset += $cigar_mod;
1699 $pos_offset += $pos_mod;
1700 }
1701
1702 ### Returning as soon as the methylation calls start overlapping
1703 if ($start+$index+$pos_offset >= $end_read_1) {
1704 return;
1705 }
1706
1707 if ($methylation_calls[$index] eq 'X') {
1708 $counting{total_meCHG_count}++;
1709 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1710 } elsif ($methylation_calls[$index] eq 'x') {
1711 $counting{total_unmethylated_CHG_count}++;
1712 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1713 } elsif ($methylation_calls[$index] eq 'Z') {
1714 $counting{total_meCpG_count}++;
1715 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1716 } elsif ($methylation_calls[$index] eq 'z') {
1717 $counting{total_unmethylated_CpG_count}++;
1718 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1719 } elsif ($methylation_calls[$index] eq 'H') {
1720 $counting{total_meCHH_count}++;
1721 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1722 } elsif ($methylation_calls[$index] eq 'h') {
1723 $counting{total_unmethylated_CHH_count}++;
1724 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1725 }
1726 elsif ($methylation_calls[$index] eq '.') {}
1727 else{
1728 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
1729 }
1730 }
1731 } elsif ($strand eq '-') {
1732 for my $index (0..$#methylation_calls) {
1733
1734 if ($cigar){ # only needed for SAM files
1735 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
1736 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
1737 $cigar_offset += $cigar_mod;
1738 $pos_offset += $pos_mod;
1739 }
1740
1741 ### Returning as soon as the methylation calls start overlapping
1742 if ($start-$index+$pos_offset <= $end_read_1) {
1743 return;
1744 }
1745
1746 if ($methylation_calls[$index] eq 'X') {
1747 $counting{total_meCHG_count}++;
1748 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1749 } elsif ($methylation_calls[$index] eq 'x') {
1750 $counting{total_unmethylated_CHG_count}++;
1751 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1752 } elsif ($methylation_calls[$index] eq 'Z') {
1753 $counting{total_meCpG_count}++;
1754 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1755 } elsif ($methylation_calls[$index] eq 'z') {
1756 $counting{total_unmethylated_CpG_count}++;
1757 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1758 } elsif ($methylation_calls[$index] eq 'H') {
1759 $counting{total_meCHH_count}++;
1760 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1761 } elsif ($methylation_calls[$index] eq 'h') {
1762 $counting{total_unmethylated_CHH_count}++;
1763 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1764 }
1765 elsif ($methylation_calls[$index] eq '.') {}
1766 else{
1767 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
1768 }
1769 }
1770 } else {
1771 die "The strand orientation was neither + nor -: '$strand'/n";
1772 }
1773 }
1774 }
1775
1776 ### this is the default paired-end procedure allowing overlaps and using every single C position
1777 ### Still within the 2-CONTEXT ONLY optional section
1778 else {
1779 ### single-file CpG and non-CpG context output
1780 if ($full) {
1781 if ($strand eq '+') {
1782 for my $index (0..$#methylation_calls) {
1783
1784 if ($cigar){ # only needed for SAM files
1785 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
1786 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
1787 $cigar_offset += $cigar_mod;
1788 $pos_offset += $pos_mod;
1789 }
1790
1791 if ($methylation_calls[$index] eq 'X') {
1792 $counting{total_meCHG_count}++;
1793 print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1794 } elsif ($methylation_calls[$index] eq 'x') {
1795 $counting{total_unmethylated_CHG_count}++;
1796 print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1797 } elsif ($methylation_calls[$index] eq 'Z') {
1798 $counting{total_meCpG_count}++;
1799 print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1800 } elsif ($methylation_calls[$index] eq 'z') {
1801 $counting{total_unmethylated_CpG_count}++;
1802 print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1803 } elsif ($methylation_calls[$index] eq 'H') {
1804 $counting{total_meCHH_count}++;
1805 print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1806 } elsif ($methylation_calls[$index] eq 'h') {
1807 $counting{total_unmethylated_CHH_count}++;
1808 print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1809 }
1810 elsif ($methylation_calls[$index] eq '.') {}
1811 else{
1812 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
1813 }
1814 }
1815 } elsif ($strand eq '-') {
1816 for my $index (0..$#methylation_calls) {
1817
1818 if ($cigar){ # only needed for SAM files
1819 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
1820 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
1821 $cigar_offset += $cigar_mod;
1822 $pos_offset += $pos_mod;
1823 }
1824
1825 if ($methylation_calls[$index] eq 'X') {
1826 $counting{total_meCHG_count}++;
1827 print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1828 } elsif ($methylation_calls[$index] eq 'x') {
1829 $counting{total_unmethylated_CHG_count}++;
1830 print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1831 } elsif ($methylation_calls[$index] eq 'Z') {
1832 $counting{total_meCpG_count}++;
1833 print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1834 } elsif ($methylation_calls[$index] eq 'z') {
1835 $counting{total_unmethylated_CpG_count}++;
1836 print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1837 } elsif ($methylation_calls[$index] eq 'H') {
1838 $counting{total_meCHH_count}++;
1839 print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1840 } elsif ($methylation_calls[$index] eq 'h') {
1841 $counting{total_unmethylated_CHH_count}++;
1842 print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1843 }
1844 elsif ($methylation_calls[$index] eq '.') {}
1845 else{
1846 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
1847 }
1848 }
1849 } else {
1850 die "The strand orientation as neither + nor -: '$strand'\n";
1851 }
1852 }
1853
1854 ### strand-specific methylation output
1855 ### still within the 2-CONTEXT optional section
1856 else {
1857 if ($strand eq '+') {
1858 for my $index (0..$#methylation_calls) {
1859
1860 if ($cigar){ # only needed for SAM files
1861 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
1862 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
1863 $cigar_offset += $cigar_mod;
1864 $pos_offset += $pos_mod;
1865 }
1866
1867 if ($methylation_calls[$index] eq 'X') {
1868 $counting{total_meCHG_count}++;
1869 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1870 } elsif ($methylation_calls[$index] eq 'x') {
1871 $counting{total_unmethylated_CHG_count}++;
1872 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1873 } elsif ($methylation_calls[$index] eq 'Z') {
1874 $counting{total_meCpG_count}++;
1875 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1876 } elsif ($methylation_calls[$index] eq 'z') {
1877 $counting{total_unmethylated_CpG_count}++;
1878 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1879 } elsif ($methylation_calls[$index] eq 'H') {
1880 $counting{total_meCHH_count}++;
1881 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1882 } elsif ($methylation_calls[$index] eq 'h') {
1883 $counting{total_unmethylated_CHH_count}++;
1884 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1885 }
1886 elsif ($methylation_calls[$index] eq '.') {}
1887 else{
1888 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
1889 }
1890 }
1891 } elsif ($strand eq '-') {
1892 for my $index (0..$#methylation_calls) {
1893
1894 if ($cigar){ # only needed for SAM files
1895 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
1896 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
1897 $cigar_offset += $cigar_mod;
1898 $pos_offset += $pos_mod;
1899 }
1900
1901 if ($methylation_calls[$index] eq 'X') {
1902 $counting{total_meCHG_count}++;
1903 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1904 } elsif ($methylation_calls[$index] eq 'x') {
1905 $counting{total_unmethylated_CHG_count}++;
1906 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1907 } elsif ($methylation_calls[$index] eq 'Z') {
1908 $counting{total_meCpG_count}++;
1909 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1910 } elsif ($methylation_calls[$index] eq 'z') {
1911 $counting{total_unmethylated_CpG_count}++;
1912 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1913 } elsif ($methylation_calls[$index] eq 'H') {
1914 $counting{total_meCHH_count}++;
1915 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1916 } elsif ($methylation_calls[$index] eq 'h') {
1917 $counting{total_unmethylated_CHH_count}++;
1918 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1919 }
1920 elsif ($methylation_calls[$index] eq '.') {}
1921 else{
1922 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
1923 }
1924 }
1925 } else {
1926 die "The strand orientation as neither + nor -: '$strand'\n";
1927 }
1928 }
1929 }
1930 }
1931
1932 ############################################
1933 ### THIS IS THE DEFAULT 3-CONTEXT OUTPUT ###
1934 ############################################
1935
1936 elsif ($no_overlap) {
1937 ### single-file CpG, CHG and CHH context output
1938 if ($full) {
1939 if ($strand eq '+') {
1940 for my $index (0..$#methylation_calls) {
1941
1942 if ($cigar){ # only needed for SAM files
1943 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
1944 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
1945 $cigar_offset += $cigar_mod;
1946 $pos_offset += $pos_mod;
1947 }
1948
1949 ### Returning as soon as the methylation calls start overlapping
1950 if ($start+$index+$pos_offset >= $end_read_1) {
1951 return;
1952 }
1953
1954 if ($methylation_calls[$index] eq 'X') {
1955 $counting{total_meCHG_count}++;
1956 print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1957 } elsif ($methylation_calls[$index] eq 'x') {
1958 $counting{total_unmethylated_CHG_count}++;
1959 print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1960 } elsif ($methylation_calls[$index] eq 'Z') {
1961 $counting{total_meCpG_count}++;
1962 print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1963 } elsif ($methylation_calls[$index] eq 'z') {
1964 $counting{total_unmethylated_CpG_count}++;
1965 print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1966 } elsif ($methylation_calls[$index] eq 'H') {
1967 $counting{total_meCHH_count}++;
1968 print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1969 } elsif ($methylation_calls[$index] eq 'h') {
1970 $counting{total_unmethylated_CHH_count}++;
1971 print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
1972 }
1973 elsif ($methylation_calls[$index] eq '.') {}
1974 else{
1975 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
1976 }
1977 }
1978 } elsif ($strand eq '-') {
1979 for my $index (0..$#methylation_calls) {
1980
1981 if ($cigar){ # only needed for SAM files
1982 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
1983 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
1984 $cigar_offset += $cigar_mod;
1985 $pos_offset += $pos_mod;
1986 }
1987
1988 ### Returning as soon as the methylation calls start overlapping
1989 if ($start-$index+$pos_offset <= $end_read_1) {
1990 return;
1991 }
1992
1993 if ($methylation_calls[$index] eq 'X') {
1994 $counting{total_meCHG_count}++;
1995 print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1996 } elsif ($methylation_calls[$index] eq 'x') {
1997 $counting{total_unmethylated_CHG_count}++;
1998 print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
1999 } elsif ($methylation_calls[$index] eq 'Z') {
2000 $counting{total_meCpG_count}++;
2001 print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2002 } elsif ($methylation_calls[$index] eq 'z') {
2003 $counting{total_unmethylated_CpG_count}++;
2004 print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2005 } elsif ($methylation_calls[$index] eq 'H') {
2006 $counting{total_meCHH_count}++;
2007 print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2008 } elsif ($methylation_calls[$index] eq 'h') {
2009 $counting{total_unmethylated_CHH_count}++;
2010 print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2011 }
2012 elsif ($methylation_calls[$index] eq '.') {}
2013 else{
2014 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2015 }
2016 }
2017 } else {
2018 die "The strand orientation as neither + nor -: '$strand'\n";
2019 }
2020 }
2021
2022 ### strand-specific methylation output
2023 else {
2024 if ($strand eq '+') {
2025 for my $index (0..$#methylation_calls) {
2026
2027 if ($cigar){ # only needed for SAM files
2028 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2029 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
2030 $cigar_offset += $cigar_mod;
2031 $pos_offset += $pos_mod;
2032 }
2033
2034 ### Returning as soon as the methylation calls start overlapping
2035 if ($start+$index+$pos_offset >= $end_read_1) {
2036 return;
2037 }
2038
2039 if ($methylation_calls[$index] eq 'X') {
2040 $counting{total_meCHG_count}++;
2041 print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2042 } elsif ($methylation_calls[$index] eq 'x') {
2043 $counting{total_unmethylated_CHG_count}++;
2044 print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2045 } elsif ($methylation_calls[$index] eq 'Z') {
2046 $counting{total_meCpG_count}++;
2047 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2048 } elsif ($methylation_calls[$index] eq 'z') {
2049 $counting{total_unmethylated_CpG_count}++;
2050 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2051 } elsif ($methylation_calls[$index] eq 'H') {
2052 $counting{total_meCHH_count}++;
2053 print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2054 } elsif ($methylation_calls[$index] eq 'h') {
2055 $counting{total_unmethylated_CHH_count}++;
2056 print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2057 }
2058 elsif ($methylation_calls[$index] eq '.') {}
2059 else{
2060 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2061 }
2062 }
2063 } elsif ($strand eq '-') {
2064 for my $index (0..$#methylation_calls) {
2065
2066 if ($cigar){ # only needed for SAM files
2067 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
2068 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2069 $cigar_offset += $cigar_mod;
2070 $pos_offset += $pos_mod;
2071 }
2072
2073 ### Returning as soon as the methylation calls start overlapping
2074 if ($start-$index+$pos_offset <= $end_read_1) {
2075 return;
2076 }
2077
2078 if ($methylation_calls[$index] eq 'X') {
2079 $counting{total_meCHG_count}++;
2080 print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2081 } elsif ($methylation_calls[$index] eq 'x') {
2082 $counting{total_unmethylated_CHG_count}++;
2083 print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2084 } elsif ($methylation_calls[$index] eq 'Z') {
2085 $counting{total_meCpG_count}++;
2086 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2087 } elsif ($methylation_calls[$index] eq 'z') {
2088 $counting{total_unmethylated_CpG_count}++;
2089 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2090 } elsif ($methylation_calls[$index] eq 'H') {
2091 $counting{total_meCHH_count}++;
2092 print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2093 } elsif ($methylation_calls[$index] eq 'h') {
2094 $counting{total_unmethylated_CHH_count}++;
2095 print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2096 }
2097 elsif ($methylation_calls[$index] eq '.') {}
2098 else{
2099 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2100 }
2101 }
2102 } else {
2103 die "The strand orientation as neither + nor -: '$strand'\n";
2104 }
2105 }
2106 }
2107
2108 ### this is the default paired-end procedure allowing overlaps and using every single C position
2109 else {
2110 ### single-file CpG, CHG and CHH context output
2111 if ($full) {
2112 if ($strand eq '+') {
2113 for my $index (0..$#methylation_calls) {
2114
2115 if ($cigar){ # only needed for SAM files
2116 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2117 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
2118 $cigar_offset += $cigar_mod;
2119 $pos_offset += $pos_mod;
2120 }
2121
2122 if ($methylation_calls[$index] eq 'X') {
2123 $counting{total_meCHG_count}++;
2124 print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2125 } elsif ($methylation_calls[$index] eq 'x') {
2126 $counting{total_unmethylated_CHG_count}++;
2127 print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2128 } elsif ($methylation_calls[$index] eq 'Z') {
2129 $counting{total_meCpG_count}++;
2130 print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2131 } elsif ($methylation_calls[$index] eq 'z') {
2132 $counting{total_unmethylated_CpG_count}++;
2133 print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2134 } elsif ($methylation_calls[$index] eq 'H') {
2135 $counting{total_meCHH_count}++;
2136 print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2137 } elsif ($methylation_calls[$index] eq 'h') {
2138 $counting{total_unmethylated_CHH_count}++;
2139 print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2140 }
2141 elsif ($methylation_calls[$index] eq '.') {}
2142 else{
2143 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2144 }
2145 }
2146 } elsif ($strand eq '-') {
2147 for my $index (0..$#methylation_calls) {
2148
2149 if ($cigar){ # only needed for SAM files
2150 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
2151 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2152 $cigar_offset += $cigar_mod;
2153 $pos_offset += $pos_mod;
2154 }
2155
2156 if ($methylation_calls[$index] eq 'X') {
2157 $counting{total_meCHG_count}++;
2158 print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2159 } elsif ($methylation_calls[$index] eq 'x') {
2160 $counting{total_unmethylated_CHG_count}++;
2161 print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2162 } elsif ($methylation_calls[$index] eq 'Z') {
2163 $counting{total_meCpG_count}++;
2164 print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2165 } elsif ($methylation_calls[$index] eq 'z') {
2166 $counting{total_unmethylated_CpG_count}++;
2167 print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2168 } elsif ($methylation_calls[$index] eq 'H') {
2169 $counting{total_meCHH_count}++;
2170 print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2171 } elsif ($methylation_calls[$index] eq 'h') {
2172 $counting{total_unmethylated_CHH_count}++;
2173 print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2174 }
2175 elsif ($methylation_calls[$index] eq '.') {}
2176 else{
2177 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2178 }
2179 }
2180 } else {
2181 die "The strand orientation as neither + nor -: '$strand'\n";
2182 }
2183 }
2184
2185 ### strand-specific methylation output
2186 else {
2187 if ($strand eq '+') {
2188 for my $index (0..$#methylation_calls) {
2189
2190 if ($cigar){ # only needed for SAM files
2191 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2192 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
2193 $cigar_offset += $cigar_mod;
2194 $pos_offset += $pos_mod;
2195 }
2196
2197 if ($methylation_calls[$index] eq 'X') {
2198 $counting{total_meCHG_count}++;
2199 print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2200 } elsif ($methylation_calls[$index] eq 'x') {
2201 $counting{total_unmethylated_CHG_count}++;
2202 print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2203 } elsif ($methylation_calls[$index] eq 'Z') {
2204 $counting{total_meCpG_count}++;
2205 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2206 } elsif ($methylation_calls[$index] eq 'z') {
2207 $counting{total_unmethylated_CpG_count}++;
2208 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2209 } elsif ($methylation_calls[$index] eq 'H') {
2210 $counting{total_meCHH_count}++;
2211 print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2212 } elsif ($methylation_calls[$index] eq 'h') {
2213 $counting{total_unmethylated_CHH_count}++;
2214 print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2215 }
2216 elsif ($methylation_calls[$index] eq '.') {}
2217 else{
2218 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2219 }
2220 }
2221 } elsif ($strand eq '-') {
2222 for my $index (0..$#methylation_calls) {
2223
2224 if ($cigar){ # only needed for SAM files
2225 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
2226 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2227 $cigar_offset += $cigar_mod;
2228 $pos_offset += $pos_mod;
2229 }
2230
2231 if ($methylation_calls[$index] eq 'X') {
2232 $counting{total_meCHG_count}++;
2233 print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2234 } elsif ($methylation_calls[$index] eq 'x') {
2235 $counting{total_unmethylated_CHG_count}++;
2236 print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2237 } elsif ($methylation_calls[$index] eq 'Z') {
2238 $counting{total_meCpG_count}++;
2239 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2240 } elsif ($methylation_calls[$index] eq 'z') {
2241 $counting{total_unmethylated_CpG_count}++;
2242 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2243 } elsif ($methylation_calls[$index] eq 'H') {
2244 $counting{total_meCHH_count}++;
2245 print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2246 } elsif ($methylation_calls[$index] eq 'h') {
2247 $counting{total_unmethylated_CHH_count}++;
2248 print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2249 }
2250 elsif ($methylation_calls[$index] eq '.') {}
2251 else{
2252 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2253 }
2254 }
2255 } else {
2256 die "The strand orientation as neither + nor -: '$strand'\n";
2257 }
2258 }
2259 }
2260 }
2261
2262 sub check_cigar_string {
2263 my ($index,$cigar_offset,$pos_offset,$strand,$comp_cigar) = @_;
2264 # print "$index\t$cigar_offset\t$pos_offset\t$strand\t";
2265 my ($new_cigar_offset,$new_pos_offset) = (0,0);
2266
2267 if ($strand eq '+') {
2268 # print "### $strand strand @$comp_cigar[$index + $cigar_offset]\t";
2269
2270 if (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'M'){ # sequence position matches the genomic position
2271 # warn "position needs no adjustment\n";
2272 }
2273
2274 elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'I'){ # insertion in the read sequence
2275 $new_pos_offset -= 1; # we need to subtract the length of inserted bases from the genomic position
2276 # warn "adjusted genomic position by -1 bp (insertion)\n";
2277 }
2278
2279 elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'D'){ # deletion in the read sequence
2280 $new_cigar_offset += 1; # the composite cigar string does no longer match the methylation call index
2281 $new_pos_offset += 1; # we need to add the length of deleted bases to get the genomic position
2282 # warn "adjusted genomic position by +1 bp (deletion). Now looping through the CIGAR string until we hit another M or I\n";
2283
2284 while ( ($index + $cigar_offset + $new_cigar_offset) < (scalar @$comp_cigar) ){
2285 if (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'M'){ # sequence position matches the genomic position
2286 # warn "position needs no adjustment\n";
2287 last;
2288 }
2289 elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'I'){
2290 $new_pos_offset -= 1; # we need to subtract the length of inserted bases from the genomic position
2291 # warn "adjusted genomic position by another -1 bp (insertion)\n";
2292 last;
2293 }
2294 elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'D'){ # deletion in the read sequence
2295 $new_cigar_offset += 1; # the composite cigar string does no longer match the methylation call index
2296 $new_pos_offset += 1; # we need to add the length of deleted bases to get the genomic position
2297 # warn "adjusted genomic position by another +1 bp (deletion)\n";
2298 }
2299 else{
2300 die "The CIGAR string contained undefined operations in addition to 'M', 'I' and 'D': '@$comp_cigar[$index + $cigar_offset + $new_cigar_offset]'\n";
2301 }
2302 }
2303 }
2304 else{
2305 die "The CIGAR string contained undefined operations in addition to 'M', 'I' and 'D': '@$comp_cigar[$index + $cigar_offset + $new_cigar_offset]'\n";
2306 }
2307 }
2308
2309 elsif ($strand eq '-') {
2310 # print "### $strand strand @$comp_cigar[$index + $cigar_offset]\t";
2311
2312 if (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'M'){ # sequence position matches the genomic position
2313 # warn "position needs no adjustment\n";
2314 }
2315
2316 elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'I'){ # insertion in the read sequence
2317 $new_pos_offset += 1; # we need to add the length of inserted bases to the genomic position
2318 # warn "adjusted genomic position by +1 bp (insertion)\n";
2319 }
2320
2321 elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'D'){ # deletion in the read sequence
2322 $new_cigar_offset += 1; # the composite cigar string does no longer match the methylation call index
2323 $new_pos_offset -= 1; # we need to subtract the length of deleted bases to get the genomic position
2324 # warn "adjusted genomic position by -1 bp (deletion). Now looping through the CIGAR string until we hit another M or I\n";
2325
2326 while ( ($index + $cigar_offset + $new_cigar_offset) < (scalar @$comp_cigar) ){
2327 if (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'M'){ # sequence position matches the genomic position
2328 # warn "Found new 'M' operation; position needs no adjustment\n";
2329 last;
2330 }
2331 elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'I'){
2332 $new_pos_offset += 1; # we need to subtract the length of inserted bases from the genomic position
2333 # warn "Found new 'I' operation; adjusted genomic position by another +1 bp (insertion)\n";
2334 last;
2335 }
2336 elsif (@$comp_cigar[$index + $cigar_offset + $new_cigar_offset] eq 'D'){ # deletion in the read sequence
2337 $new_cigar_offset += 1; # the composite cigar string does no longer match the methylation call index
2338 $new_pos_offset -= 1; # we need to subtract the length of deleted bases to get the genomic position
2339 # warn "adjusted genomic position by another -1 bp (deletion)\n";
2340 }
2341 else{
2342 die "The CIGAR string contained undefined operations in addition to 'M', 'I' and 'D': '@$comp_cigar[$index + $cigar_offset + $new_cigar_offset]'\n";
2343 }
2344 }
2345 }
2346 else{
2347 die "The CIGAR string contained undefined operations in addition to 'M', 'I' and 'D': '@$comp_cigar[$index + $cigar_offset + $new_cigar_offset]'\n";
2348 }
2349 }
2350 # print "new cigar offset: $new_cigar_offset\tnew pos offset: $new_pos_offset\n";
2351 return ($new_cigar_offset,$new_pos_offset);
2352 }
2353
2354 sub print_individual_C_methylation_states_single_end{
2355
2356 my ($meth_call,$chrom,$start,$id,$strand,$filehandle_index,$cigar) = @_;
2357 my @methylation_calls = split(//,$meth_call);
2358
2359 #################################################################
2360 ### . for bases not involving cytosines ###
2361 ### X for methylated C in CHG context (was protected) ###
2362 ### x for not methylated C in CHG context (was converted) ###
2363 ### H for methylated C in CHH context (was protected) ###
2364 ### h for not methylated C in CHH context (was converted) ###
2365 ### Z for methylated C in CpG context (was protected) ###
2366 ### z for not methylated C in CpG context (was converted) ###
2367 #################################################################
2368
2369 my $methyl_CHG_count = 0;
2370 my $methyl_CHH_count = 0;
2371 my $methyl_CpG_count = 0;
2372 my $unmethylated_CHG_count = 0;
2373 my $unmethylated_CHH_count = 0;
2374 my $unmethylated_CpG_count = 0;
2375
2376 my @len;
2377 my @ops;
2378 my $pos_offset = 0; # this is only relevant for SAM reads with insertions or deletions
2379 my $cigar_offset = 0; # again, this is only relevant for SAM reads containing indels
2380
2381 my @comp_cigar;
2382
2383 if ($cigar){ # parsing CIGAR string
2384 @len = split (/\D+/,$cigar); # storing the length per operation
2385 @ops = split (/\d+/,$cigar); # storing the operation
2386 shift @ops; # remove the empty first element
2387 die "CIGAR string contained a non-matching number of lengths and operations\n" unless (scalar @len == scalar @ops);
2388
2389 foreach my $index (0..$#len){
2390 foreach (1..$len[$index]){
2391 # print "$ops[$index]";
2392 push @comp_cigar, $ops[$index];
2393 }
2394 }
2395 # warn "\nDetected CIGAR string: $cigar\n";
2396 # warn "Length of methylation call: ",length $meth_call,"\n";
2397 # warn "number of operations: ",scalar @ops,"\n";
2398 # warn "number of length digits: ",scalar @len,"\n\n";
2399 # print @comp_cigar,"\n";
2400 # print "$meth_call\n\n";
2401 # sleep (1);
2402 }
2403
2404 ### adjusting the start position for all reads mapping to the reverse strand
2405 if ($strand eq '-') {
2406
2407 @comp_cigar = reverse@comp_cigar; # the CIGAR string needs to be reversed for all reads aligning to the reverse strand, too
2408 # print @comp_cigar,"\n";
2409
2410 unless ($ignore){ ### if --ignore was specified the start position has already been corrected
2411
2412 if ($cigar){ ### SAM format
2413 my $MD_count = 0;
2414 foreach (@comp_cigar){
2415 ++$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
2416 }
2417 $start += $MD_count - 1;
2418 }
2419 else{ ### vanilla format
2420 $start += length($meth_call)-1;
2421 }
2422 }
2423 }
2424
2425 ### THIS IS THE CpG and Non-CpG SECTION (OPTIONAL)
2426
2427 ### single-file CpG and other-context output
2428 if ($full and $merge_non_CpG) {
2429 if ($strand eq '+') {
2430 for my $index (0..$#methylation_calls) {
2431
2432 if ($cigar){ # only needed for SAM files
2433 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2434 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition+index: ",$start+$index,"\t";
2435 $cigar_offset += $cigar_mod;
2436 $pos_offset += $pos_mod;
2437 }
2438
2439 ### methylated Cs (any context) will receive a forward (+) orientation
2440 ### not methylated Cs (any context) will receive a reverse (-) orientation
2441 if ($methylation_calls[$index] eq 'X') {
2442 $counting{total_meCHG_count}++;
2443 print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2444 }
2445 elsif ($methylation_calls[$index] eq 'x') {
2446 $counting{total_unmethylated_CHG_count}++;
2447 print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2448 }
2449 elsif ($methylation_calls[$index] eq 'Z') {
2450 $counting{total_meCpG_count}++;
2451 print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2452 }
2453 elsif ($methylation_calls[$index] eq 'z') {
2454 $counting{total_unmethylated_CpG_count}++;
2455 print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2456 }
2457 elsif ($methylation_calls[$index] eq 'H') {
2458 $counting{total_meCHH_count}++;
2459 print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2460 }
2461 elsif ($methylation_calls[$index] eq 'h') {
2462 $counting{total_unmethylated_CHH_count}++;
2463 print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2464 }
2465 elsif ($methylation_calls[$index] eq '.') {
2466 }
2467 else{
2468 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2469 }
2470 }
2471 }
2472 elsif ($strand eq '-') {
2473
2474 for my $index (0..$#methylation_calls) {
2475 ### methylated Cs (any context) will receive a forward (+) orientation
2476 ### not methylated Cs (any context) will receive a reverse (-) orientation
2477
2478 if ($cigar){ # only needed for SAM files
2479 # print "index: $index\tmethylation_call: $methylation_calls[$index]\tposition-index: ",$start-$index,"\t";
2480 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2481 $cigar_offset += $cigar_mod;
2482 $pos_offset += $pos_mod;
2483 }
2484
2485 if ($methylation_calls[$index] eq 'X') {
2486 $counting{total_meCHG_count}++;
2487 print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2488 }
2489 elsif ($methylation_calls[$index] eq 'x') {
2490 $counting{total_unmethylated_CHG_count}++;
2491 print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2492 }
2493 elsif ($methylation_calls[$index] eq 'Z') {
2494 $counting{total_meCpG_count}++;
2495 print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2496 }
2497 elsif ($methylation_calls[$index] eq 'z') {
2498 $counting{total_unmethylated_CpG_count}++;
2499 print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2500 }
2501 elsif ($methylation_calls[$index] eq 'H') {
2502 $counting{total_meCHH_count}++;
2503 print {$fhs{other_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2504 }
2505 elsif ($methylation_calls[$index] eq 'h') {
2506 $counting{total_unmethylated_CHH_count}++;
2507 print {$fhs{other_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2508 }
2509 elsif ($methylation_calls[$index] eq '.'){
2510 }
2511 else{
2512 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2513 }
2514 }
2515 }
2516 else {
2517 die "The strand information was neither + nor -: $strand\n";
2518 }
2519 }
2520
2521 ### strand-specific methylation output
2522 elsif ($merge_non_CpG) {
2523 if ($strand eq '+') {
2524 for my $index (0..$#methylation_calls) {
2525 ### methylated Cs (any context) will receive a forward (+) orientation
2526 ### not methylated Cs (any context) will receive a reverse (-) orientation
2527
2528 if ($cigar){ # only needed for SAM files
2529 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2530 $cigar_offset += $cigar_mod;
2531 $pos_offset += $pos_mod;
2532 }
2533
2534 if ($methylation_calls[$index] eq 'X') {
2535 $counting{total_meCHG_count}++;
2536 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2537 }
2538 elsif ($methylation_calls[$index] eq 'x') {
2539 $counting{total_unmethylated_CHG_count}++;
2540 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2541 }
2542 elsif ($methylation_calls[$index] eq 'Z') {
2543 $counting{total_meCpG_count}++;
2544 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2545 }
2546 elsif ($methylation_calls[$index] eq 'z') {
2547 $counting{total_unmethylated_CpG_count}++;
2548 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2549 }
2550 elsif ($methylation_calls[$index] eq 'H') {
2551 $counting{total_meCHH_count}++;
2552 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2553 }
2554 elsif ($methylation_calls[$index] eq 'h') {
2555 $counting{total_unmethylated_CHH_count}++;
2556 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2557 }
2558 elsif ($methylation_calls[$index] eq '.') {
2559 }
2560 else{
2561 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2562 }
2563 }
2564 }
2565 elsif ($strand eq '-') {
2566
2567 for my $index (0..$#methylation_calls) {
2568 ### methylated Cs (any context) will receive a forward (+) orientation
2569 ### not methylated Cs (any context) will receive a reverse (-) orientation
2570
2571 if ($cigar){ # only needed for SAM files
2572 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2573 $cigar_offset += $cigar_mod;
2574 $pos_offset += $pos_mod;
2575 }
2576
2577 if ($methylation_calls[$index] eq 'X') {
2578 $counting{total_meCHG_count}++;
2579 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2580 }
2581 elsif ($methylation_calls[$index] eq 'x') {
2582 $counting{total_unmethylated_CHG_count}++;
2583 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2584 }
2585 elsif ($methylation_calls[$index] eq 'Z') {
2586 $counting{total_meCpG_count}++;
2587 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2588 }
2589 elsif ($methylation_calls[$index] eq 'z') {
2590 $counting{total_unmethylated_CpG_count}++;
2591 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2592 }
2593 elsif ($methylation_calls[$index] eq 'H') {
2594 $counting{total_meCHH_count}++;
2595 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2596 }
2597 elsif ($methylation_calls[$index] eq 'h') {
2598 $counting{total_unmethylated_CHH_count}++;
2599 print {$fhs{$filehandle_index}->{other_c}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2600 }
2601 elsif ($methylation_calls[$index] eq '.') {
2602 }
2603 else{
2604 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2605 }
2606 }
2607 }
2608 else {
2609 die "The strand information was neither + nor -: $strand\n";
2610 }
2611 }
2612
2613 ### THIS IS THE 3-CONTEXT (CpG, CHG and CHH) DEFAULT SECTION
2614
2615 elsif ($full) {
2616 if ($strand eq '+') {
2617 for my $index (0..$#methylation_calls) {
2618 ### methylated Cs (any context) will receive a forward (+) orientation
2619 ### not methylated Cs (any context) will receive a reverse (-) orientation
2620
2621 if ($cigar){ # only needed for SAM files
2622 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2623 $cigar_offset += $cigar_mod;
2624 $pos_offset += $pos_mod;
2625 }
2626
2627 if ($methylation_calls[$index] eq 'X') {
2628 $counting{total_meCHG_count}++;
2629 print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2630 } elsif ($methylation_calls[$index] eq 'x') {
2631 $counting{total_unmethylated_CHG_count}++;
2632 print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2633 } elsif ($methylation_calls[$index] eq 'Z') {
2634 $counting{total_meCpG_count}++;
2635 print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2636 } elsif ($methylation_calls[$index] eq 'z') {
2637 $counting{total_unmethylated_CpG_count}++;
2638 print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2639 } elsif ($methylation_calls[$index] eq 'H') {
2640 $counting{total_meCHH_count}++;
2641 print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2642 } elsif ($methylation_calls[$index] eq 'h') {
2643 $counting{total_unmethylated_CHH_count}++;
2644 print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2645 }
2646 elsif ($methylation_calls[$index] eq '.') {}
2647 else{
2648 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2649 }
2650 }
2651 }
2652 elsif ($strand eq '-') {
2653
2654 for my $index (0..$#methylation_calls) {
2655 ### methylated Cs (any context) will receive a forward (+) orientation
2656 ### not methylated Cs (any context) will receive a reverse (-) orientation
2657
2658 if ($cigar){ # only needed for SAM files
2659 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2660 $cigar_offset += $cigar_mod;
2661 $pos_offset += $pos_mod;
2662 }
2663
2664 if ($methylation_calls[$index] eq 'X') {
2665 $counting{total_meCHG_count}++;
2666 print {$fhs{CHG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2667 } elsif ($methylation_calls[$index] eq 'x') {
2668 $counting{total_unmethylated_CHG_count}++;
2669 print {$fhs{CHG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2670 } elsif ($methylation_calls[$index] eq 'Z') {
2671 $counting{total_meCpG_count}++;
2672 print {$fhs{CpG_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2673 } elsif ($methylation_calls[$index] eq 'z') {
2674 $counting{total_unmethylated_CpG_count}++;
2675 print {$fhs{CpG_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2676 } elsif ($methylation_calls[$index] eq 'H') {
2677 $counting{total_meCHH_count}++;
2678 print {$fhs{CHH_context}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2679 } elsif ($methylation_calls[$index] eq 'h') {
2680 $counting{total_unmethylated_CHH_count}++;
2681 print {$fhs{CHH_context}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2682 }
2683 elsif ($methylation_calls[$index] eq '.') {}
2684 else{
2685 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2686 }
2687 }
2688 }
2689 else {
2690 die "The read had a strand orientation which was neither + nor -: $strand\n";
2691 }
2692 }
2693
2694 ### strand-specific methylation output
2695 else {
2696 if ($strand eq '+') {
2697 for my $index (0..$#methylation_calls) {
2698 ### methylated Cs (any context) will receive a forward (+) orientation
2699 ### not methylated Cs (any context) will receive a reverse (-) orientation
2700
2701 if ($cigar){ # only needed for SAM files
2702 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2703 $cigar_offset += $cigar_mod;
2704 $pos_offset += $pos_mod;
2705 }
2706
2707 if ($methylation_calls[$index] eq 'X') {
2708 $counting{total_meCHG_count}++;
2709 print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2710 } elsif ($methylation_calls[$index] eq 'x') {
2711 $counting{total_unmethylated_CHG_count}++;
2712 print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2713 } elsif ($methylation_calls[$index] eq 'Z') {
2714 $counting{total_meCpG_count}++;
2715 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2716 } elsif ($methylation_calls[$index] eq 'z') {
2717 $counting{total_unmethylated_CpG_count}++;
2718 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2719 } elsif ($methylation_calls[$index] eq 'H') {
2720 $counting{total_meCHH_count}++;
2721 print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2722 } elsif ($methylation_calls[$index] eq 'h') {
2723 $counting{total_unmethylated_CHH_count}++;
2724 print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start+$index+$pos_offset,$methylation_calls[$index]),"\n";
2725 }
2726 elsif ($methylation_calls[$index] eq '.') {}
2727 else{
2728 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2729 }
2730 }
2731 }
2732 elsif ($strand eq '-') {
2733
2734 for my $index (0..$#methylation_calls) {
2735 ### methylated Cs (any context) will receive a forward (+) orientation
2736 ### not methylated Cs (any context) will receive a reverse (-) orientation
2737
2738 if ($cigar){ # only needed for SAM files
2739 my ($cigar_mod,$pos_mod) = check_cigar_string($index,$cigar_offset,$pos_offset,$strand,\@comp_cigar);
2740 $cigar_offset += $cigar_mod;
2741 $pos_offset += $pos_mod;
2742 }
2743
2744 if ($methylation_calls[$index] eq 'X') {
2745 $counting{total_meCHG_count}++;
2746 print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2747 } elsif ($methylation_calls[$index] eq 'x') {
2748 $counting{total_unmethylated_CHG_count}++;
2749 print {$fhs{$filehandle_index}->{CHG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2750 } elsif ($methylation_calls[$index] eq 'Z') {
2751 $counting{total_meCpG_count}++;
2752 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2753 } elsif ($methylation_calls[$index] eq 'z') {
2754 $counting{total_unmethylated_CpG_count}++;
2755 print {$fhs{$filehandle_index}->{CpG}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2756 } elsif ($methylation_calls[$index] eq 'H') {
2757 $counting{total_meCHH_count}++;
2758 print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'+',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2759 } elsif ($methylation_calls[$index] eq 'h') {
2760 $counting{total_unmethylated_CHH_count}++;
2761 print {$fhs{$filehandle_index}->{CHH}} join ("\t",$id,'-',$chrom,$start-$index+$pos_offset,$methylation_calls[$index]),"\n";
2762 }
2763 elsif ($methylation_calls[$index] eq '.') {}
2764 else{
2765 die "The methylation call string contained the following unrecognised character: $methylation_calls[$index]\n";
2766 }
2767 }
2768 }
2769 else {
2770 die "The strand information was neither + nor -: $strand\n";
2771 }
2772 }
2773 }
2774
2775
2776
2777 #######################################################################################################################################
2778 ### bismark2bedGaph section - START
2779 #######################################################################################################################################
2780
2781 sub process_bedGraph_output{
2782 warn "="x64,"\n";
2783 warn "Methylation information will now be written into a bedGraph file\n";
2784 warn "="x64,"\n\n";
2785 sleep (2);
2786
2787 ### Closing all filehandles so that the Bismark methtylation extractor output doesn't get truncated due to buffering issues
2788 foreach my $fh (keys %fhs) {
2789 if ($fh =~ /^[1230]$/) {
2790 foreach my $context (keys %{$fhs{$fh}}) {
2791 close $fhs{$fh}->{$context} or die $!;
2792 }
2793 } else {
2794 close $fhs{$fh} or die $!;
2795 }
2796 }
2797
2798 ### deciding which files to use for bedGraph conversion
2799 foreach my $filename (@sorting_files){
2800 # warn "$filename\n";
2801 if ($filename =~ /\//){ # if files are in a different output folder we extract the filename again
2802 $filename =~ s/.*\///; # replacing everything up to the last slash in the filename
2803 # warn "$filename\n";
2804 }
2805
2806 if ($CX_context){
2807 push @bedfiles,$filename;
2808 }
2809 else{ ## CpG context only (default)
2810 if ($filename =~ /^CpG_/){
2811 push @bedfiles,$filename;
2812 }
2813 else{
2814 # skipping CHH or CHG files
2815 }
2816 }
2817 }
2818
2819 warn "Using the following files as Input:\n";
2820 print join ("\t",@bedfiles),"\n\n";
2821 sleep (2);
2822
2823 my %temp_fhs;
2824 my @temp_files; # writing all context files (default CpG only) to these files prior to sorting
2825
2826 ### changing to the output directory
2827 unless ($output_dir eq ''){ # default
2828 chdir $output_dir or die "Failed to change directory to $output_dir\n";
2829 warn "Changed directory to $output_dir\n";
2830 }
2831
2832 foreach my $infile (@bedfiles) {
2833
2834 if ($remove) {
2835 warn "Now replacing whitespaces in the sequence ID field of the Bismark methylation extractor output $infile prior to bedGraph conversion\n\n";
2836
2837 open (READ,$infile) or die $!;
2838
2839 my $removed_spaces_outfile = $infile;
2840 $removed_spaces_outfile =~ s/$/.spaces_removed.txt/;
2841
2842 open (REM,'>',$output_dir.$removed_spaces_outfile) or die "Couldn't write to file $removed_spaces_outfile: $!\n";
2843
2844 unless ($no_header){
2845 $_ = <READ>; ### Bismark version header
2846 print REM $_; ### Bismark version header
2847 }
2848
2849 while (<READ>) {
2850 chomp;
2851 my ($id,$strand,$chr,$pos,$context) = (split (/\t/));
2852 $id =~ s/\s+/_/g;
2853 print REM join ("\t",$id,$strand,$chr,$pos,$context),"\n";
2854 }
2855
2856 close READ or die $!;
2857 close REM or die $!;
2858
2859 ### changing the infile name to the new file without spaces
2860 $infile = $removed_spaces_outfile;
2861 }
2862
2863 warn "Now writing methylation information for file $infile to individual files for each chromosome\n";
2864 open (IN,$infile) or die $!;
2865
2866 ## always ignoring the version header
2867 unless ($no_header){
2868 $_ = <IN>; ### Bismark version header
2869 }
2870
2871 while (<IN>) {
2872 chomp;
2873 my ($chr) = (split (/\t/))[2];
2874
2875 unless (exists $temp_fhs{$chr}) {
2876 open ($temp_fhs{$chr},'>','chr'.$chr.'.meth_extractor.temp') or die "Failed to open filehandle: $!";
2877 }
2878 print {$temp_fhs{$chr}} "$_\n";
2879 }
2880
2881 warn "Finished writing out individual chromosome files for $infile\n";
2882 }
2883 warn "\n";
2884
2885 @temp_files = <*.meth_extractor.temp>;
2886
2887 warn "Collecting temporary chromosome file information...\n";
2888 sleep (1);
2889 warn "processing the following input file(s):\n";
2890 warn join ("\n",@temp_files),"\n\n";
2891 sleep (1);
2892
2893 foreach my $in (@temp_files) {
2894 warn "Sorting input file $in by positions\n";
2895 open my $ifh, "sort -k3,3 -k4,4n $in |" or die "Input file could not be sorted. $!";
2896 # print "Chromosome\tStart Position\tEnd Position\tMethylation Percentage\n";
2897
2898 ############################################# m.a.bentley - moved the variables out of the while loop to hold the current line data {
2899
2900 my $name;
2901 my $meth_state;
2902 my $chr = "";
2903 my $pos = 0;
2904 my $meth_state2;
2905
2906 my $last_pos;
2907 my $last_chr;
2908
2909 ############################################# }
2910
2911 while (my $line = <$ifh>) {
2912 next if $line =~ /^Bismark/;
2913 chomp $line;
2914
2915 ########################################### 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 {
2916
2917 $last_chr = $chr;
2918 $last_pos = $pos;
2919 ($name, $meth_state, $chr, $pos, $meth_state2) = split "\t", $line;
2920
2921 if (($last_pos ne $pos) || ($last_chr ne $chr)) {
2922 generate_output($last_chr,$last_pos) if $methylcalls[2] > 0;
2923 @methylcalls = qw (0 0 0);
2924 }
2925
2926 ############################################# }
2927
2928 my $validated = validate_methylation_call($meth_state, $meth_state2);
2929 unless($validated){
2930 warn "Methylation state of sequence ($name) in file ($in) on line $. is inconsistent (meth_state is $meth_state, meth_state2 = $meth_state2)\n";
2931 next;
2932 }
2933 if ($meth_state eq "+") {
2934 $methylcalls[0]++;
2935 $methylcalls[2]++;
2936 } else {
2937 $methylcalls[1]++;
2938 $methylcalls[2]++;
2939 }
2940 }
2941
2942 ############################################# 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) {
2943
2944 $last_chr = $chr;
2945 $last_pos = $pos;
2946 if ($methylcalls[2] > 0) {
2947 generate_output($last_chr,$last_pos) if $methylcalls[2] > 0;
2948 }
2949 ############################################# }
2950
2951 close $ifh or die $!;
2952
2953 @methylcalls = qw (0 0 0); # resetting @methylcalls
2954
2955 ### deleting temporary files
2956 my $delete = unlink $in;
2957 if ($delete) {
2958 warn "Successfully deleted the temporary input file $in\n\n";
2959 }
2960 else {
2961 warn "The temporary inputfile $in could not be deleted $!\n\n";
2962 }
2963 }
2964 }
2965
2966 sub generate_output{
2967 my $methcount = $methylcalls[0];
2968 my $nonmethcount = $methylcalls[1];
2969 my $totalcount = $methylcalls[2];
2970 my $last_chr = shift;
2971 my $last_pos = shift;
2972 croak "Should not be generating output if there's no reads to this region" unless $totalcount > 0;
2973 croak "Total counts ($totalcount) is not the sum of the methylated ($methcount) and unmethylated ($nonmethcount) counts" if $totalcount != ($methcount + $nonmethcount);
2974
2975 ############################################# 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 {
2976
2977 my $bed_pos = $last_pos -1; ### Bismark coordinates are 1 based whereas bedGraph coordinates are 0 based.
2978 my $meth_percentage;
2979 ($totalcount >= $coverage_threshold) ? ($meth_percentage = ($methcount/$totalcount) * 100) : ($meth_percentage = undef);
2980 # $meth_percentage =~ s/(\.\d\d).+$/$1/ unless $meth_percentage =~ /^Below/;
2981 if (defined $meth_percentage){
2982 if ($counts){
2983 print OUT "$last_chr\t$bed_pos\t$bed_pos\t$meth_percentage\t$methcount\t$nonmethcount\n";
2984 }
2985 else{
2986 print OUT "$last_chr\t$bed_pos\t$bed_pos\t$meth_percentage\n";
2987 }
2988 }
2989 ############################################# }
2990 }
2991
2992 sub validate_methylation_call{
2993 my $meth_state = shift;
2994 croak "Missing (+/-) methylation call" unless defined $meth_state;
2995 my $meth_state2 = shift;
2996 croak "Missing alphabetical methylation call" unless defined $meth_state2;
2997 my $is_consistent;
2998 ($meth_state2 =~ /^z/i) ? ($is_consistent = check_CpG_methylation_call($meth_state, $meth_state2))
2999 : ($is_consistent = check_nonCpG_methylation_call($meth_state,$meth_state2));
3000 return 1 if $is_consistent;
3001 return 0;
3002 }
3003
3004 sub check_CpG_methylation_call{
3005 my $meth1 = shift;
3006 my $meth2 = shift;
3007 return 1 if($meth1 eq "+" && $meth2 eq "Z");
3008 return 1 if($meth1 eq "-" && $meth2 eq "z");
3009 return 0;
3010 }
3011
3012 sub check_nonCpG_methylation_call{
3013 my $meth1 = shift;
3014 my $meth2 = shift;
3015 return 1 if($meth1 eq "+" && $meth2 eq "C");
3016 return 1 if($meth1 eq "+" && $meth2 eq "X");
3017 return 1 if($meth1 eq "+" && $meth2 eq "H");
3018 return 1 if($meth1 eq "-" && $meth2 eq "c");
3019 return 1 if($meth1 eq "-" && $meth2 eq "x");
3020 return 1 if($meth1 eq "-" && $meth2 eq "h");
3021 return 0;
3022 }
3023
3024 #######################################################################################################################################
3025 ### bismark2bedGaph section - END
3026 #######################################################################################################################################
3027
3028
3029
3030
3031
3032
3033 #######################################################################################################################################
3034 ### genome-wide cytosine methylation report - START
3035 #######################################################################################################################################
3036
3037 sub generate_genome_wide_cytosine_report {
3038
3039 warn "="x78,"\n";
3040 warn "Methylation information will now be written into a genome-wide cytosine report\n";
3041 warn "="x78,"\n\n";
3042 sleep (2);
3043
3044 ### changing to the output directory again
3045 unless ($output_dir eq ''){ # default
3046 chdir $output_dir or die "Failed to change directory to $output_dir\n";
3047 # warn "Changed directory to $output_dir\n";
3048 }
3049
3050 my $in = shift;
3051 open (IN,$in) or die $!;
3052
3053 my $cytosine_out = shift;
3054
3055 if ($CX_context){
3056 $cytosine_out =~ s/$/genome-wide_CX_report.txt/;
3057 }
3058 else{
3059 $cytosine_out =~ s/$/genome_wide_CpG_report.txt/;
3060 }
3061
3062 ### note: we are still in the folder: $output_dir, so we do not have to include this into the open commands
3063 unless ($split_by_chromosome){ ### writing all output to a single file (default)
3064 open (CYT,'>',$cytosine_out) or die $!;
3065 warn "Writing genome-wide cytosine report to: $cytosine_out\n";
3066 sleep (3);
3067 }
3068
3069 my $last_chr;
3070 my %chr; # storing reads for one chromosome at a time
3071
3072 my $count = 0;
3073 while (<IN>){
3074 chomp;
3075 ++$count;
3076 my ($chr,$start,$end,undef,$meth,$nonmeth) = (split /\t/);
3077
3078 # defining the first chromosome
3079 unless (defined $last_chr){
3080 $last_chr = $chr;
3081 # warn "Storing all covered cytosine positions for chromosome: $chr\n";
3082 }
3083
3084 if ($chr eq $last_chr){
3085 $chr{$chr}->{$start}->{meth} = $meth;
3086 $chr{$chr}->{$start}->{nonmeth} = $nonmeth;
3087 }
3088 else{
3089 warn "Writing cytosine reports for chromosome $last_chr (stored ",scalar keys %{$chr{$last_chr}}," different covered positions)\n";
3090
3091 if ($split_by_chromosome){ ## writing output to 1 file per chromosome
3092 my $chromosome_out = $cytosine_out;
3093 $chromosome_out =~ s/txt$/chr${last_chr}.txt/;
3094 open (CYT,'>',$chromosome_out) or die $!;
3095 }
3096
3097 while ( $chromosomes{$last_chr} =~ /([CG])/g){
3098
3099 my $tri_nt = '';
3100 my $context = '';
3101 my $pos = pos$chromosomes{$last_chr};
3102
3103 my $strand;
3104 my $meth = 0;
3105 my $nonmeth = 0;
3106
3107 if ($1 eq 'C'){ # C on forward strand
3108 $tri_nt = substr ($chromosomes{$last_chr},($pos-1),3); # positions are 0-based!
3109 $strand = '+';
3110 }
3111 elsif ($1 eq 'G'){ # C on reverse strand
3112 $tri_nt = substr ($chromosomes{$last_chr},($pos-3),3); # positions are 0-based!
3113 $tri_nt = reverse $tri_nt;
3114 $tri_nt =~ tr/ACTG/TGAC/;
3115 $strand = '-';
3116 }
3117 next if (length$tri_nt < 3); # trinucleotide sequence could not be extracted
3118
3119 if (exists $chr{$last_chr}->{($pos-1)}){ # stored positions are 0-based!
3120 $meth = $chr{$last_chr}->{$pos-1}->{meth};
3121 $nonmeth = $chr{$last_chr}->{$pos-1}->{nonmeth};
3122 }
3123
3124 ### determining cytosine context
3125 if ($tri_nt =~ /^CG/){
3126 $context = 'CG';
3127 }
3128 elsif ($tri_nt =~ /^C.{1}G$/){
3129 $context = 'CHG';
3130 }
3131 elsif ($tri_nt =~ /^C.{2}$/){
3132 $context = 'CHH';
3133 }
3134 else{ # if the context can't be determined the positions will not be printed (it will equally not have been reported by Bismark)
3135 warn "The sequence context could not be determined (found: '$tri_nt'). Skipping.\n";
3136 next;
3137 }
3138
3139 if ($CpG_only){
3140 if ($tri_nt =~ /^CG/){ # CpG context is the default
3141 if ($zero){ # zero based coordinates
3142 $pos -= 1;
3143 print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
3144 }
3145 else{ # default
3146 print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
3147 }
3148 }
3149 }
3150 else{ ## all cytosines, specified with --CX
3151 if ($zero){ # zero based coordinates
3152 $pos -= 1;
3153 print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
3154 }
3155 else{ # default
3156 print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
3157 }
3158 }
3159 }
3160
3161 %chr = (); # resetting the hash
3162
3163 # new first entry
3164 $last_chr = $chr;
3165 $chr{$chr}->{$start}->{meth} = $meth;
3166 $chr{$chr}->{$start}->{nonmeth} = $nonmeth;
3167 }
3168 }
3169
3170 # Last found chromosome
3171 warn "Writing cytosine reports for chromosome $last_chr (stored ",scalar keys %{$chr{$last_chr}}," different covered positions)\n";
3172
3173 if ($split_by_chromosome){ ## writing output to 1 file per chromosome
3174 my $chromosome_out = $cytosine_out;
3175 $chromosome_out =~ s/txt$/chr${last_chr}.txt/;
3176 open (CYT,'>',$chromosome_out) or die $!;
3177 }
3178
3179 while ( $chromosomes{$last_chr} =~ /([CG])/g){
3180
3181 my $tri_nt;
3182 my $context;
3183 my $pos = pos$chromosomes{$last_chr};
3184
3185 my $strand;
3186 my $meth = 0;
3187 my $nonmeth = 0;
3188
3189 if ($1 eq 'C'){ # C on forward strand
3190 $tri_nt = substr ($chromosomes{$last_chr},($pos-1),3); # positions are 0-based!
3191 $strand = '+';
3192 }
3193 elsif ($1 eq 'G'){ # C on reverse strand
3194 $tri_nt = substr ($chromosomes{$last_chr},($pos-3),3); # positions are 0-based!
3195 $tri_nt = reverse $tri_nt;
3196 $tri_nt =~ tr/ACTG/TGAC/;
3197 $strand = '-';
3198 }
3199
3200 if (exists $chr{$last_chr}->{($pos-1)}){ # stored positions are 0-based!
3201 $meth = $chr{$last_chr}->{$pos-1}->{meth};
3202 $nonmeth = $chr{$last_chr}->{$pos-1}->{nonmeth};
3203 }
3204
3205 next if (length$tri_nt < 3); # trinucleotide sequence could not be extracted
3206
3207 ### determining cytosine context
3208 if ($tri_nt =~ /^CG/){
3209 $context = 'CG';
3210 }
3211 elsif ($tri_nt =~ /^C.{1}G$/){
3212 $context = 'CHG';
3213 }
3214 elsif ($tri_nt =~ /^C.{2}$/){
3215 $context = 'CHH';
3216 }
3217 else{ # if the context can't be determined the positions will not be printed (it will equally not have been reported by Bismark)
3218 warn "The cytosine context could not be determined (found: '$tri_nt'). Skipping.\n";
3219 next;
3220 }
3221
3222 if ($CpG_only){
3223 if ($tri_nt =~ /^CG/){ # CpG context is the default
3224 if ($zero){ # zero-based coordinates
3225 $pos -= 1;
3226 print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
3227 }
3228 else{ # default
3229 print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
3230 }
3231 }
3232 }
3233 else{ ## all cytosines, specified with --CX
3234 if ($zero){ # zero based coordinates
3235 $pos -= 1;
3236 print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
3237 }
3238 else{ # default
3239 print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n";
3240 }
3241 }
3242 }
3243 close CYT or die $!;
3244 }
3245
3246
3247 sub read_genome_into_memory{
3248
3249 ## reading in and storing the specified genome in the %chromosomes hash
3250 chdir ($genome_folder) or die "Can't move to $genome_folder: $!";
3251 warn "Now reading in and storing sequence information of the genome specified in: $genome_folder\n\n";
3252
3253 my @chromosome_filenames = <*.fa>;
3254
3255 ### if there aren't any genomic files with the extension .fa we will look for files with the extension .fasta
3256 unless (@chromosome_filenames){
3257 @chromosome_filenames = <*.fasta>;
3258 }
3259 unless (@chromosome_filenames){
3260 die "The specified genome folder $genome_folder does not contain any sequence files in FastA format (with .fa or .fasta file extensions)\n";
3261 }
3262
3263 foreach my $chromosome_filename (@chromosome_filenames){
3264
3265 # skipping the tophat entire mouse genome fasta file
3266 next if ($chromosome_filename eq 'Mus_musculus.NCBIM37.fa');
3267
3268 open (CHR_IN,$chromosome_filename) or die "Failed to read from sequence file $chromosome_filename $!\n";
3269 ### first line needs to be a fastA header
3270 my $first_line = <CHR_IN>;
3271 chomp $first_line;
3272 $first_line =~ s/\r//; # removing /r carriage returns
3273
3274 ### Extracting chromosome name from the FastA header
3275 my $chromosome_name = extract_chromosome_name($first_line);
3276
3277 my $sequence;
3278 while (<CHR_IN>){
3279 chomp;
3280 $_ =~ s/\r//; # removing /r carriage returns
3281
3282 if ($_ =~ /^>/){
3283 ### storing the previous chromosome in the %chromosomes hash, only relevant for Multi-Fasta-Files (MFA)
3284 if (exists $chromosomes{$chromosome_name}){
3285 warn "chr $chromosome_name (",length $sequence ," bp)\n";
3286 die "Exiting because chromosome name already exists. Please make sure all chromosomes have a unique name!\n";
3287 }
3288 else {
3289 if (length($sequence) == 0){
3290 warn "Chromosome $chromosome_name in the multi-fasta file $chromosome_filename did not contain any sequence information!\n";
3291 }
3292 warn "chr $chromosome_name (",length $sequence ," bp)\n";
3293 $chromosomes{$chromosome_name} = $sequence;
3294 }
3295 ### resetting the sequence variable
3296 $sequence = '';
3297 ### setting new chromosome name
3298 $chromosome_name = extract_chromosome_name($_);
3299 }
3300 else{
3301 $sequence .= uc$_;
3302 }
3303 }
3304
3305 if (exists $chromosomes{$chromosome_name}){
3306 warn "chr $chromosome_name (",length $sequence ," bp)\t";
3307 die "Exiting because chromosome name already exists. Please make sure all chromosomes have a unique name.\n";
3308 }
3309 else{
3310 if (length($sequence) == 0){
3311 warn "Chromosome $chromosome_name in the file $chromosome_filename did not contain any sequence information!\n";
3312 }
3313 warn "chr $chromosome_name (",length $sequence ," bp)\n";
3314 $chromosomes{$chromosome_name} = $sequence;
3315 }
3316 }
3317 warn "\n";
3318 chdir $parent_dir or die "Failed to move to directory $parent_dir\n";
3319 }
3320
3321 sub extract_chromosome_name {
3322 ## Bowtie extracts the first string after the inition > in the FASTA file, so we are doing this as well
3323 my $fasta_header = shift;
3324 if ($fasta_header =~ s/^>//){
3325 my ($chromosome_name) = split (/\s+/,$fasta_header);
3326 return $chromosome_name;
3327 }
3328 else{
3329 die "The specified chromosome ($fasta_header) file doesn't seem to be in FASTA format as required!\n";
3330 }
3331 }
3332
3333 #######################################################################################################################################
3334 ### genome-wide cytosine methylation report - END
3335 #######################################################################################################################################
3336
3337
3338
3339
3340 sub print_helpfile{
3341
3342 print << 'HOW_TO';
3343
3344
3345 DESCRIPTION
3346
3347 The following is a brief description of all options to control the Bismark
3348 methylation extractor. The script reads in a bisulfite read alignment results file
3349 produced by the Bismark bisulfite mapper and extracts the methylation information
3350 for individual cytosines. This information is found in the methylation call field
3351 which can contain the following characters:
3352
3353 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3354 ~~~ X for methylated C in CHG context (was protected) ~~~
3355 ~~~ x for not methylated C CHG (was converted) ~~~
3356 ~~~ H for methylated C in CHH context (was protected) ~~~
3357 ~~~ h for not methylated C in CHH context (was converted) ~~~
3358 ~~~ Z for methylated C in CpG context (was protected) ~~~
3359 ~~~ z for not methylated C in CpG context (was converted) ~~~
3360 ~~~ . for any bases not involving cytosines ~~~
3361 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3362
3363 The methylation extractor outputs result files for cytosines in CpG, CHG and CHH
3364 context (this distinction is actually already made in Bismark itself). As the methylation
3365 information for every C analysed can produce files which easily have tens or even hundreds of
3366 millions of lines, file sizes can become very large and more difficult to handle. The C
3367 methylation info additionally splits cytosine methylation calls up into one of the four possible
3368 strands a given bisulfite read aligned against:
3369
3370 OT original top strand
3371 CTOT complementary to original top strand
3372
3373 OB original bottom strand
3374 CTOB complementary to original bottom strand
3375
3376 Thus, by default twelve individual output files are being generated per input file (unless
3377 --comprehensive is specified, see below). The output files can be imported into a genome
3378 viewer, such as SeqMonk, and re-combined into a single data group if desired (in fact
3379 unless the bisulfite reads were generated preserving directionality it doesn't make any
3380 sense to look at the data in a strand-specific manner). Strand-specific output files can
3381 optionally be skipped, in which case only three output files for CpG, CHG or CHH context
3382 will be generated. For both the strand-specific and comprehensive outputs there is also
3383 the option to merge both non-CpG contexts (CHG and CHH) into one single non-CpG context.
3384
3385
3386 The output files are in the following format (tab delimited):
3387
3388 <sequence_id> <strand> <chromosome> <position> <methylation call>
3389
3390
3391 USAGE: methylation_extractor [options] <filenames>
3392
3393
3394 ARGUMENTS:
3395
3396 <filenames> A space-separated list of Bismark result files in SAM format from
3397 which methylation information is extracted for every cytosine in
3398 the reads. For alignment files in the older custom Bismark output
3399 see option '--vanilla'.
3400
3401 OPTIONS:
3402
3403 -s/--single-end Input file(s) are Bismark result file(s) generated from single-end
3404 read data. Specifying either --single-end or --paired-end is
3405 mandatory.
3406
3407 -p/--paired-end Input file(s) are Bismark result file(s) generated from paired-end
3408 read data. Specifying either --paired-end or --single-end is
3409 mandatory.
3410
3411 --vanilla The Bismark result input file(s) are in the old custom Bismark format
3412 (up to version 0.5.x) and not in SAM format which is the default as
3413 of Bismark version 0.6.x or higher. Default: OFF.
3414
3415 --no_overlap For paired-end reads it is theoretically possible that read_1 and
3416 read_2 overlap. This option avoids scoring overlapping methylation
3417 calls twice (only methylation calls of read 1 are used for in the process
3418 since read 1 has historically higher quality basecalls than read 2).
3419 Whilst this option removes a bias towards more methylation calls
3420 in the center of sequenced fragments it may de facto remove a sizable
3421 proportion of the data. This option is highly recommended for paired-end
3422 data.
3423
3424 --ignore <int> Ignore the first <int> bp at the 5' end of each read when processing the
3425 methylation call string. This can remove e.g. a restriction enzyme site
3426 at the start of each read.
3427
3428 --comprehensive Specifying this option will merge all four possible strand-specific
3429 methylation info into context-dependent output files. The default
3430 contexts are:
3431 - CpG context
3432 - CHG context
3433 - CHH context
3434
3435 --merge_non_CpG This will produce two output files (in --comprehensive mode) or eight
3436 strand-specific output files (default) for Cs in
3437 - CpG context
3438 - non-CpG context
3439
3440 --report Prints out a short methylation summary as well as the paramaters used to run
3441 this script.
3442
3443 --no_header Suppresses the Bismark version header line in all output files for more convenient
3444 batch processing.
3445
3446 -o/--output DIR Allows specification of a different output directory (absolute or relative
3447 path). If not specified explicitely, the output will be written to the current directory.
3448
3449 --version Displays version information.
3450
3451 -h/--help Displays this help file and exits.
3452
3453
3454
3455 bedGraph specific options:
3456
3457 --bedGraph After finishing the methylation extraction, the methylation output is written into a
3458 sorted bedGraph file that reports the position of a given cytosine and its methylation
3459 state (in %, seem details below). The methylation extractor output is temporarily split up into
3460 temporary files, one per chromosome (written into the current directory or folder
3461 specified with -o/--output); these temp files are then used for sorting and deleted
3462 afterwards. By default, only cytosines in CpG context will be sorted. The option
3463 '--CX_context' may be used to report all cyosines irrespective of sequence context
3464 (this will take MUCH longer!).
3465
3466
3467 --cutoff [threshold] The minimum number of times a methylation state has to be seen for that nucleotide
3468 before its methylation percentage is reported. Default: 1.
3469
3470 --remove_spaces Replaces whitespaces in the sequence ID field with underscores to allow sorting.
3471
3472
3473 --counts Adds two additional columns to the output file to enable further calculations:
3474 col 5: number of methylated calls
3475 col 6: number of unmethylated calls
3476 This option is required if '--cytosine_report' is specified (and will be set automatically if
3477 necessary).
3478
3479 --CX/--CX_context The sorted bedGraph output file contains information on every single cytosine that was covered
3480 in the experiment irrespective of its sequence context. This applies to both forward and
3481 reverse strands. Please be aware that this option may generate large temporary and output files
3482 and may take a long time to sort (up to many hours). Default: OFF.
3483 (i.e. Default = CpG context only).
3484
3485
3486
3487 Genome-wide cytosine methylation report specific options:
3488
3489 --cytosine_report After the conversion to bedGraph has completed, the option '--cytosine_report' produces a
3490 genome-wide methylation report for all cytosines in the genome. By default, the output uses 1-based
3491 chromosome coordinates (zero-based cords are optional) and reports CpG context only (all
3492 cytosine context is optional). The output considers all Cs on both forward and reverse strands and
3493 reports their position, strand, trinucleotide content and methylation state (counts are 0 if not
3494 covered).
3495
3496 --CX/--CX_context The output file contains information on every single cytosine in the genome irrespective of
3497 its context. This applies to both forward and reverse strands. Please be aware that this will
3498 generate output files with > 1.1 billion lines for a mammalian genome such as human or mouse.
3499 Default: OFF (i.e. Default = CpG context only).
3500
3501 --zero_based Uses zero-based coordinates like used in e.g. bed files instead of 1-based coordinates. Default: OFF.
3502
3503 --genome_folder <path> Enter the genome folder you wish to use to extract sequences from (full path only). Accepted
3504 formats are FastA files ending with '.fa' or '.fasta'. Specifying a genome folder path is mandatory.
3505
3506 --split_by_chromosome Writes the output into individual files for each chromosome instead of a single output file. Files
3507 will be named to include the input filename and the chromosome number.
3508
3509
3510
3511 OUTPUT:
3512
3513 The bismark_methylation_extractor output is in the form:
3514 ========================================================
3515 <seq-ID> <methylation state*> <chromosome> <start position (= end position)> <methylation call>
3516
3517 * Methylated cytosines receive a '+' orientation,
3518 * Unmethylated cytosines receive a '-' orientation.
3519
3520
3521
3522 The bedGraph output (optional) looks like this (tab-delimited):
3523 ===============================================================
3524 <chromosome> <start position> <end position> <methylation percentage>
3525
3526
3527
3528 The genome-wide cytosine methylation output file is tab-delimited in the following format:
3529 ==========================================================================================
3530 <chromosome> <position> <strand> <count methylated> <count non-methylated> <C-context> <trinucleotide context>
3531
3532
3533
3534 This script was last modified on 02 Oct 2012.
3535
3536 HOW_TO
3537 }