Mercurial > repos > bgruening > bismark
diff bismark_methyl_extractor/coverage2cytosine @ 7:fcadce4d9a06 draft
planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/tools/bismark commit b'e6ee273f75fff61d1e419283fa8088528cf59470\n'
| author | bgruening |
|---|---|
| date | Sat, 06 May 2017 13:18:09 -0400 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bismark_methyl_extractor/coverage2cytosine Sat May 06 13:18:09 2017 -0400 @@ -0,0 +1,1254 @@ +#!/usr/bin/env perl +use warnings; +use strict; +$|++; +use Getopt::Long; +use Cwd; +use Carp; + +## This program is Copyright (C) 2010-16, Felix Krueger (felix.krueger@babraham.ac.uk) + +## This program is free software: you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation, either version 3 of the License, or +## (at your option) any later version. + +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. + +## You should have received a copy of the GNU General Public License +## along with this program. If not, see <http://www.gnu.org/licenses/>. + +my %chromosomes; # storing sequence information of all chromosomes/scaffolds +my %processed; # keeping a record of which chromosomes have been processed +my $coverage2cytosine_version = 'v0.16.3'; + +my ($output_dir,$genome_folder,$zero,$CpG_only,$CX_context,$split_by_chromosome,$parent_dir,$coverage_infile,$cytosine_out,$merge_CpGs,$gc_context,$gzip,$tetra) = process_commandline(); + +warn "Summary of parameters for genome-wide cytosine report:\n"; +warn '='x78,"\n"; +warn "Coverage infile:\t\t\t$coverage_infile\n"; +warn "Output directory:\t\t\t>$output_dir<\n"; +warn "Parent directory:\t\t\t>$parent_dir<\n"; +warn "Genome directory:\t\t\t>$genome_folder<\n"; + +if ($CX_context){ + warn "CX context:\t\t\t\tyes\n"; +} +else{ + warn "CX context:\t\t\t\tno (CpG context only, default)\n"; +} +if ($merge_CpGs){ + warn "Pooling CpG top/bottom strand evidence:\tyes\n"; +} +if($gc_context){ + warn "Optional GC context track:\t\tyes\n"; +} +if ($tetra){ + warn "Tetra/Penta nucleotide context:\t\tyes\n"; +} + +if ($zero){ + warn "Genome coordinates used:\t\t0-based (user specified)\n"; +} +else{ + warn "Genome coordinates used:\t\t1-based (default)\n"; +} + +if ($gzip){ + warn "GZIP compression:\t\t\tyes\n"; +} +else{ + warn "GZIP compression:\t\t\tno\n"; +} + +if ($split_by_chromosome){ + warn "Split by chromosome:\t\t\tyes\n\n\n"; +} +else{ + warn "Split by chromosome:\t\t\tno\n\n\n"; +} +sleep (3); + +read_genome_into_memory(); +warn "Stored sequence information of ",scalar keys %chromosomes," chromosomes/scaffolds in total\n\n"; + +generate_genome_wide_cytosine_report($coverage_infile); + +### 11 December 2014 + +# The following optional section re-reads the genome-wide report and merges methylation evidence of both top and bottom strand +# into a single CpG dinucleotide entity. This significantly simplifies downstream processing, e.g. by the bsseq R-/Bioconductor package +# which recommends this merging process to increase coverage per CpG and reduce the memory burden for its processing + +if ($merge_CpGs){ + # we only allow this operation if the report is limited to CpG context, and for a single report for the entire genome for the time being + combine_CpGs_to_single_CG_entity($cytosine_out); +} + +### 18 August 2015 + +# The following section reprocessed the genome to generate cytosine methylation output in GC context (e.g. when a GpC methylase had been deployed +if ($gc_context){ + generate_GC_context_report($coverage_infile); +} + + +sub combine_CpGs_to_single_CG_entity{ + my $CpG_report_file = shift; + warn "Now merging top and bottom strand CpGs into a single CG dinucleotide entity\n"; + + open (IN,$CpG_report_file) or die "Failed to open file $CpG_report_file: $!\n\n"; + my $pooled_CG = $CpG_report_file; + $pooled_CG =~ s/$/.merged_CpG_evidence.cov/; + open (OUT,'>',$pooled_CG) or die "Failed to write to file '$pooled_CG': $!\n\n"; + warn ">>> Writing a new coverage file with top and bottom strand CpG methylation evidence merged to $pooled_CG <<<\n\n"; + sleep(1); + + while (1){ + my $line1 = <IN>; + my $line2 = <IN>; + last unless ($line1 and $line2); + + chomp $line1; + chomp $line2; + + my ($chr1,$pos1,$strand1,$m1,$u1,$context1) = (split /\t/,$line1); + my ($chr2,$pos2,$strand2,$m2,$u2,$context2) = (split /\t/,$line2); + # print join ("\t",$chr1,$pos1,$strand1,$m1,$u1,$context1),"\n"; + # print join ("\t",$chr2,$pos2,$strand2,$m2,$u2,$context2),"\n"; sleep(1); + + if ($pos1 < 2){ + $line1 = $line2; + $line2 = <IN>; + chomp $line2; + ($chr1,$pos1,$strand1,$m1,$u1,$context1) = (split /\t/,$line1); + ($chr2,$pos2,$strand2,$m2,$u2,$context2) = (split /\t/,$line2); + } + + # a few sanity checks + die "The context of the first line was not CG:\t$line1" unless ($context1 eq 'CG'); + die "The context of the second line was not CG:\t$line2" unless ($context2 eq 'CG'); + + unless ($strand1 eq '+' and $strand2 eq'-'){ + die "The strand of line 1 and line 2 were not + and -:\nline1:\t$line1\nline2:\t$line2\n"; + } + unless ($pos2 eq ($pos1 + 1)){ + die "The reported position were not spaced 1bp apart:\nline1:\t$line1\nline2:\t$line2\n"; + } + unless($chr1 eq $chr2){ + die "The chromsome information for line 1 and 2 did not match:\nline1:\t$line1\nline2:\t$line2\n"; + } + + # This looks alright from what I can tell, lets pool the 2 strands + + my $pooled_m = $m1 + $m2; + my $pooled_u = $u1 + $u2; + + # since this is a new coverage file we only write it out if the coverage is at least 1 + next unless ($pooled_m + $pooled_u) > 0; + + my $pooled_percentage = sprintf ("%.6f",$pooled_m/ ($pooled_m + $pooled_u) *100); + # print join ("\t",$chr1,$pos1,$strand1,$m1,$u1,$context1),"\n"; + # print join ("\t",$chr2,$pos2,$strand2,$m2,$u2,$context2),"\n"; + # print join ("\t",$chr1,$pos1,$pos2,$pooled_percentage,$pooled_m,$pooled_u),"\n"; + # sleep(1); + print OUT join ("\t",$chr1,$pos1,$pos2,$pooled_percentage,$pooled_m,$pooled_u),"\n"; + } + + warn "CpG merging complete. Good luck finding DMRs with bsseq!\n\n"; + +} + + +sub process_commandline{ + my $help; + my $output_dir; + my $genome_folder; + my $zero; + my $CpG_only; + my $CX_context; + my $gzip; + + my $split_by_chromosome; + my $cytosine_out; + my $parent_dir; + my $version; + my $merge_CpGs; + my $gc_context; + my $tetra; + + my $command_line = GetOptions ('help|man' => \$help, + 'dir=s' => \$output_dir, + 'g|genome_folder=s' => \$genome_folder, + "zero_based" => \$zero, + "CX|CX_context" => \$CX_context, + "split_by_chromosome" => \$split_by_chromosome, + 'o|output=s' => \$cytosine_out, + 'parent_dir=s' => \$parent_dir, + 'version' => \$version, + 'merge_CpGs' => \$merge_CpGs, + 'GC|GC_context' => \$gc_context, + 'ff|qq' => \$tetra, + 'gzip' => \$gzip, + ); + + ### EXIT ON ERROR if there were errors with any of the supplied options + unless ($command_line){ + die "Please respecify command line options\n"; + } + + ### HELPFILE + if ($help){ + print_helpfile(); + exit; + } + + if ($version){ + print << "VERSION"; + + + Bismark Methylation Extractor Module - + coverage2cytosine + + Bismark Extractor Version: $coverage2cytosine_version + Copyright 2010-15 Felix Krueger, Babraham Bioinformatics + www.bioinformatics.babraham.ac.uk/projects/bismark/ + + +VERSION + exit; + } + + ### no files provided + unless (@ARGV){ + warn "You need to provide a Bismark coverage file (with counts methylated/unmethylated cytosines) to create an individual C methylation output. Please respecify!\n"; + sleep(2); + + print_helpfile(); + exit; + } + + my $coverage_infile = shift @ARGV; + + unless ($parent_dir){ + $parent_dir = getcwd(); + } + unless ($parent_dir =~ /\/$/){ + $parent_dir =~ s/$/\//; + } + + unless (defined $cytosine_out){ + die "Please provide the name of the output file using the option -o/--output filename\n"; + } + + ### OUTPUT DIR PATH + if (defined $output_dir){ + unless ($output_dir eq ''){ # if the output dir has been passed on by the methylation extractor and is an empty string we don't want to change it + unless ($output_dir =~ /\/$/){ + $output_dir =~ s/$/\//; + } + } + } + else{ + $output_dir = ''; + } + + unless ($CX_context){ + $CX_context = 0; + $CpG_only = 1; + } + + ### GENOME folder + if ($genome_folder){ + unless ($genome_folder =~/\/$/){ + $genome_folder =~ s/$/\//; + } + } + else{ + die "Please specify a genome folder to proceed (full path only)\n"; + } + + if ($merge_CpGs){ + if ($CX_context){ + die "Merging individual CpG calls into a single CpG dinucleotide entity is currently only supported if CpG-context is selected only (lose the option --CX)\n"; + } + if ($split_by_chromosome){ + die "Merging individual CpG calls into a single CpG dinucleotide entity is currently only supported if a single CpG report is written out (lose the option --split_by_chromosome)\n"; + } + } + + return ($output_dir,$genome_folder,$zero,$CpG_only,$CX_context,$split_by_chromosome,$parent_dir,$coverage_infile,$cytosine_out,$merge_CpGs,$gc_context,$gzip,$tetra); +} + + + +sub generate_genome_wide_cytosine_report { + + warn "="x78,"\n"; + warn "Methylation information will now be written into a genome-wide cytosine report\n"; + warn "="x78,"\n\n"; + sleep (2); + + my $number_processed = 0; + + ### changing to the output directory again + unless ($output_dir eq ''){ # default + chdir $output_dir or die "Failed to change directory to $output_dir\n"; + # warn "Changed directory to $output_dir\n"; + } + + my $in = shift; + # infiles handed over by the methylation extractor will be just the filename on their own. The directory should have been handed over with --dir + if ($in =~ /gz$/){ + open (IN,"gunzip -c $in |") or die "Failed to read from gzipped file $in: $!\n"; # changed from gunzip -c to gunzip -c 08 04 16 + } + else{ + open (IN,"$in") or die "Failed to read from file $in: $!\n"; + } + + + ### note: we are still in the folder: $output_dir, so we do not have to include this into the open commands + unless ($split_by_chromosome){ ### writing all output to a single file (default) + if ($gzip){ + unless ($cytosine_out =~ /\.gz$/){ + $cytosine_out .= '.gz'; + } + open (CYT,"| gzip -c - > $cytosine_out") or die "Failed to write to file $cytosine_out: $!\n"; + } + else{ + open (CYT,'>',$cytosine_out) or die "Failed to write to file $cytosine_out: $!\n"; + } + warn ">>> Writing genome-wide cytosine report to: $cytosine_out <<<\n\n"; + sleep (1); + } + + my $last_chr; + my %chr; # storing reads for one chromosome at a time + + my $count = 0; + while (<IN>){ + chomp; + ++$count; + my ($chr,$start,$end,undef,$meth,$nonmeth) = (split /\t/); + + # defining the first chromosome + unless (defined $last_chr){ + $last_chr = $chr; + ++$number_processed; + # warn "Storing all covered cytosine positions for chromosome: $chr\n"; + } + + ### As of version 0.9.1 the start positions are 1-based! + if ($chr eq $last_chr){ + $chr{$chr}->{$start}->{meth} = $meth; + $chr{$chr}->{$start}->{nonmeth} = $nonmeth; + } + else{ + warn "Writing cytosine report for chromosome $last_chr (stored ",scalar keys %{$chr{$last_chr}}," different covered positions)\n"; + ++$number_processed; + + if ($split_by_chromosome){ ## writing output to 1 file per chromosome + my $chromosome_out = $cytosine_out; + if ($chromosome_out =~ /\.txt$/){ + $chromosome_out =~ s/\.txt$/chr${last_chr}.txt/; + } + else{ + $chromosome_out =~ s/$/.chr${last_chr}.txt/; + } + open (CYT,'>',$chromosome_out) or die $!; + # warn "Writing output for $last_chr to $chromosome_out\n"; + } + + my $tri_nt; + my $tetra_nt; + my $penta_nt; + my $context; + + $processed{$last_chr} = 1; + while ( $chromosomes{$last_chr} =~ /([CG])/g){ + + $tri_nt = ''; + $context = ''; + + if ($tetra){ + $tetra_nt = ''; # clearing + $penta_nt = ''; + } + + my $pos = pos$chromosomes{$last_chr}; + + my $strand; + my $meth = 0; + my $nonmeth = 0; + + if ($1 eq 'C'){ # C on forward strand + $tri_nt = substr ($chromosomes{$last_chr},($pos-1),3); # positions are 0-based! + if ($tetra){ + if ( length($chromosomes{$last_chr}) >= ($pos - 1 + 4) ){ + $tetra_nt = substr ($chromosomes{$last_chr},($pos-1),4); + } + else{ + $tetra_nt = ''; + } + if ( length($chromosomes{$last_chr}) >= ($pos - 1 + 5) ){ + $penta_nt = substr ($chromosomes{$last_chr},($pos-1),5); + } + else{ + $penta_nt = ''; + } + } + $strand = '+'; + } + elsif ($1 eq 'G'){ # C on reverse strand + + $tri_nt = substr ($chromosomes{$last_chr},($pos-3),3); # positions are 0-based! + $tri_nt = reverse $tri_nt; + $tri_nt =~ tr/ACTG/TGAC/; + + if ($tetra){ + if ( $pos - 4 >= 0 ){ + $tetra_nt = substr ($chromosomes{$last_chr},($pos-4),4); + $tetra_nt = reverse $tetra_nt; + $tetra_nt =~ tr/ACTG/TGAC/; + } + else{ + $tetra_nt = ''; + } + if ( $pos - 5 >= 0 ){ + $penta_nt = substr ($chromosomes{$last_chr},($pos-5),5); + $penta_nt = reverse $penta_nt; + $penta_nt =~ tr/ACTG/TGAC/; + } + else{ + $penta_nt = ''; + } + } + $strand = '-'; + } + + next if (length$tri_nt < 3); # trinucleotide sequence could not be extracted + + # if (length$penta_nt < 5){ + # warn "$tri_nt\t$tetra_nt\t$penta_nt\n"; sleep(1); + # } + + if (exists $chr{$last_chr}->{($pos)}){ # stored positions are 1-based! (as of v0.9.1) + $meth = $chr{$last_chr}->{$pos}->{meth}; + $nonmeth = $chr{$last_chr}->{$pos}->{nonmeth}; + } + + + ### determining cytosine context + if ($tri_nt =~ /^CG/){ + $context = 'CG'; + } + elsif ($tri_nt =~ /^C.{1}G$/){ + $context = 'CHG'; + } + elsif ($tri_nt =~ /^C.{2}$/){ + $context = 'CHH'; + } + else{ # if the context can't be determined the positions will not be printed (it will equally not have been reported by Bismark) + warn "The sequence context could not be determined (found: '$tri_nt'). Skipping.\n"; + next; + } + + if ($CpG_only){ + if ($tri_nt =~ /^CG/){ # CpG context is the default + if ($zero){ # zero based coordinates + $pos -= 1; + if ($tetra){ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt,$tetra_nt,$penta_nt),"\n"; + } + else{ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n"; + } + } + else{ # default + if ($tetra){ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt,$tetra_nt,$penta_nt),"\n"; + } + else{ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n"; + } + } + } + } + else{ ## all cytosines, specified with --CX + if ($zero){ # zero based coordinates + $pos -= 1; + if ($tetra){ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt,$tetra_nt,$penta_nt),"\n"; + } + else{ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n" + } + } + else{ # default + if ($tetra){ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt,$tetra_nt,$penta_nt),"\n"; + } + else{ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n"; + } + } + } + } + + %chr = (); # resetting the hash + + # new first entry + $last_chr = $chr; + $chr{$chr}->{$start}->{meth} = $meth; + $chr{$chr}->{$start}->{nonmeth} = $nonmeth; + } + } + + # Last found chromosome + # If there never was a last chromosome then something must have gone wrong with reading the data in + unless (defined $last_chr){ + die "No last chromosome was defined, something must have gone wrong while reading the data in (e.g. specified wrong file path for a gzipped coverage file?). Please check your command!\n\n"; + } + + warn "Writing cytosine report for last chromosome $last_chr (stored ",scalar keys %{$chr{$last_chr}}," different covered positions)\n"; + $processed{$last_chr} = 1; + + if ($split_by_chromosome){ ## writing output to 1 file per chromosome + my $chromosome_out = $cytosine_out; + if ($chromosome_out =~ /\.txt$/){ # files passed on by the methylation extractor end in _report.txt + $chromosome_out =~ s/\.txt$/chr${last_chr}.txt/; + } + else{ # user specified output file name + $chromosome_out =~ s/$/.chr${last_chr}.txt/; + } + open (CYT,'>',$chromosome_out) or die $!; + # warn "Writing output for $last_chr to $chromosome_out\n"; + } + + my $tri_nt; + my $tetra_nt; + my $penta_nt; + my $context; + + while ( $chromosomes{$last_chr} =~ /([CG])/g){ + + $tri_nt = ''; + $context = ''; + + if ($tetra){ + $tetra_nt = ''; + $penta_nt = ''; + } + + my $pos = pos$chromosomes{$last_chr}; + + my $strand; + my $meth = 0; + my $nonmeth = 0; + + if ($1 eq 'C'){ # C on forward strand + $tri_nt = substr ($chromosomes{$last_chr},($pos-1),3); # positions are 0-based! + $strand = '+'; + + if ($tetra){ + if ( length($chromosomes{$last_chr}) >= ($pos - 1 + 4) ){ + $tetra_nt = substr ($chromosomes{$last_chr},($pos-1),4); + } + else{ + $tetra_nt = ''; + } + if ( length($chromosomes{$last_chr}) >= ($pos - 1 + 5) ){ + $penta_nt = substr ($chromosomes{$last_chr},($pos-1),5); + } + else{ + $penta_nt = ''; + } + } + + } + elsif ($1 eq 'G'){ # C on reverse strand + $tri_nt = substr ($chromosomes{$last_chr},($pos-3),3); # positions are 0-based! + $tri_nt = reverse $tri_nt; + $tri_nt =~ tr/ACTG/TGAC/; + $strand = '-'; + + if ($tetra){ + if ( $pos - 4 >= 0 ){ + $tetra_nt = substr ($chromosomes{$last_chr},($pos-4),4); + $tetra_nt = reverse $tetra_nt; + $tetra_nt =~ tr/ACTG/TGAC/; + } + else{ + $tetra_nt = ''; + } + + if ( $pos - 5 >= 0 ){ + $penta_nt = substr ($chromosomes{$last_chr},($pos-5),5); + $penta_nt = reverse $penta_nt; + $penta_nt =~ tr/ACTG/TGAC/; + } + else{ + $penta_nt = ''; + } + } + } + + if (exists $chr{$last_chr}->{($pos)}){ # stored positions are 1-based! as of v0.9.1 + $meth = $chr{$last_chr}->{$pos}->{meth}; + $nonmeth = $chr{$last_chr}->{$pos}->{nonmeth}; + } + + next if (length$tri_nt < 3); # trinucleotide sequence could not be extracted + # if (length$penta_nt < 5){ + # warn "$tri_nt\t$tetra_nt\t$penta_nt\n"; sleep(1); + # } + + ### determining cytosine context + if ($tri_nt =~ /^CG/){ + $context = 'CG'; + } + elsif ($tri_nt =~ /^C.{1}G$/){ + $context = 'CHG'; + } + elsif ($tri_nt =~ /^C.{2}$/){ + $context = 'CHH'; + } + else{ # if the context can't be determined the positions will not be printed (it will equally not have been reported by Bismark) + warn "The cytosine context could not be determined (found: '$tri_nt'). Skipping.\n"; + next; + } + + if ($CpG_only){ + if ($tri_nt =~ /^CG/){ # CpG context is the default + if ($zero){ # zero-based coordinates + $pos -= 1; + if ($tetra){ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt,$tetra_nt,$penta_nt),"\n"; + } + else{ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n"; + } + } + else{ # default + if ($tetra){ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt,$tetra_nt,$penta_nt),"\n"; + } + else{ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n"; + } + } + } + } + else{ ## all cytosines, specified with --CX + if ($zero){ # zero based coordinates + $pos -= 1; + if ($tetra){ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt,$tetra_nt,$penta_nt),"\n"; + } + else{ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n"; + } + } + else{ # default + if ($tetra){ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt,$tetra_nt,$penta_nt),"\n"; + } + else{ + print CYT join ("\t",$last_chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n"; + } + } + } + } + + warn "Finished writing out cytosine report for covered chromosomes (processed $number_processed chromosomes/scaffolds in total)\n\n"; + + ### Now processing chromosomes that were not covered in the coverage file + warn "Now processing chromosomes that were not covered by any methylation calls in the coverage file...\n"; + my $unprocessed = 0; + + foreach my $chr (sort keys %processed) { + unless ( $processed{$chr} ) { + ++$unprocessed; + ++$number_processed; + process_unprocessed_chromosomes($chr); + } + } + + if ($unprocessed == 0) { + warn "All chromosomes in the genome were covered by at least some reads. coverage2cytosine processing complete.\n\n"; + } + else{ + warn "Finished writing out cytosine report (processed $number_processed chromosomes/scaffolds in total). coverage2cytosine processing complete.\n\n"; + } + + close CYT or warn $!; + +} + + +#### GC CONTEXT - optional +#### + +sub generate_GC_context_report { + + warn "="x82,"\n"; + warn "Methylation information for GC context will now be written to a GpC-context report\n"; + warn "="x82,"\n\n"; + sleep (2); + + my $number_processed = 0; + + ### changing to the output directory again + unless ($output_dir eq ''){ # default + chdir $output_dir or die "Failed to change directory to $output_dir\n"; + # warn "Changed directory to $output_dir\n"; + } + + my $in = shift; + if ($in =~ /gz$/){ + open (IN,"gunzip -c $in |") or die "Failed to read from gzipped file $in: $!\n"; + } + else{ + open (IN,"$in") or die "Failed to read from file $in: $!\n"; + } + + my $gc_out = $cytosine_out.'.GpC_report.txt'; + my $gc_cov = $cytosine_out.'.GpC.cov'; + + ### note: we are still in the folder: $output_dir, so we do not have to include this into the open commands + open (GC,'>',$gc_out) or die "Failed to write to GpC file $gc_out: !\n\n"; + warn ">>> Writing genome-wide GpC cytosine report to: $gc_out <<<\n"; + open (GCCOV,'>',$gc_cov) or die "Failed to write to GpC coverage file $gc_cov: !\n\n"; + warn ">>> Writing genome-wide GpC coverage file to: $gc_cov <<<\n\n"; + + my $last_chr; + my %chr; # storing reads for one chromosome at a time + + my $count = 0; + while (<IN>){ + chomp; + ++$count; + my ($chr,$start,$end,undef,$meth,$nonmeth) = (split /\t/); + + # defining the first chromosome + unless (defined $last_chr){ + $last_chr = $chr; + ++$number_processed; + # warn "Storing all covered cytosine positions for chromosome: $chr\n"; + } + + ### start positions are 1-based + if ($chr eq $last_chr){ + $chr{$chr}->{$start}->{meth} = $meth; + $chr{$chr}->{$start}->{nonmeth} = $nonmeth; + } + else{ + warn "Writing cytosine report for chromosome $last_chr (stored ",scalar keys %{$chr{$last_chr}}," different covered positions)\n"; + ++$number_processed; + $processed{$last_chr} = 1; + while ( $chromosomes{$last_chr} =~ /(GC)/g){ + + # a GC on the top strand automatically means that there is a GC on the bottom strand as well, so we can process both at the same time + my $context_top = ''; + my $context_bottom = ''; + my $pos = pos$chromosomes{$last_chr}; + + my $meth_top = 0; + my $meth_bottom = 0; + my $nonmeth_top = 0; + my $nonmeth_bottom = 0; + + #warn "$1\n"; sleep(1); + # C on forward strand + my $tri_nt_top = substr ($chromosomes{$last_chr},($pos-1),3); # positions are 0-based! + my $strand_top = '+'; + + # C on reverse strand + my $tri_nt_bottom = substr ($chromosomes{$last_chr},($pos-4),3); # positions are 0-based! + $tri_nt_bottom = reverse $tri_nt_bottom; + $tri_nt_bottom =~ tr/ACTG/TGAC/; + my $strand_bottom = '-'; + + next if (length $tri_nt_top < 3); # trinucleotide sequence could not be extracted for the top strand + next if (length $tri_nt_bottom < 3); # trinucleotide sequence could not be extracted for the reverse strand + + # top strand + if (exists $chr{$last_chr}->{($pos)}){ # stored positions are 1-based! + $meth_top = $chr{$last_chr}->{$pos}->{meth}; + $nonmeth_top = $chr{$last_chr}->{$pos}->{nonmeth}; + } + # bottom strand + if (exists $chr{$last_chr}->{($pos-1)}){ # stored positions are 1-based! -1 for bottom strand Cs + $meth_bottom = $chr{$last_chr}->{$pos-1}->{meth}; + $nonmeth_bottom = $chr{$last_chr}->{$pos-1}->{nonmeth}; + } + + ### determining cytosine context top strand + if ($tri_nt_top =~ /^CG/){ + $context_top = 'CG'; + } + elsif ($tri_nt_top =~ /^C.{1}G$/){ + $context_top = 'CHG'; + } + elsif ($tri_nt_top =~ /^C.{2}$/){ + $context_top = 'CHH'; + } + else{ # if the context can't be determined the positions will not be printed (it will equally not have been reported by Bismark) + warn "The sequence context could not be determined for the top strand (found: '$tri_nt_top'). Skipping.\n"; + next; + } + + ### determining cytosine context bottom strand + if ($tri_nt_bottom =~ /^CG/){ + $context_bottom = 'CG'; + } + elsif ($tri_nt_bottom =~ /^C.{1}G$/){ + $context_bottom = 'CHG'; + } + elsif ($tri_nt_bottom =~ /^C.{2}$/){ + $context_bottom = 'CHH'; + } + else{ # if the context can't be determined the positions will not be printed + warn "The sequence context could not be determined for the bottom strand (found: '$tri_nt_bottom'). Skipping.\n"; + next; + } + + # if Cs were not covered at all they are not written out + unless ($meth_bottom + $nonmeth_bottom == 0){ + my $percentage = sprintf ("%.6f",$meth_bottom/ ($meth_bottom + $nonmeth_bottom) *100); + print GCCOV join ("\t",$last_chr,$pos-1,$pos-1,$percentage,$meth_bottom,$nonmeth_bottom),"\n"; + print GC join ("\t",$last_chr,$pos-1,$strand_bottom,$meth_bottom,$nonmeth_bottom,$context_bottom,$tri_nt_bottom),"\n"; + } + unless ($meth_top + $nonmeth_top == 0){ + my $percentage = sprintf ("%.6f",$meth_top/ ($meth_top + $nonmeth_top) *100); + print GCCOV join ("\t",$last_chr,$pos,$pos,$percentage,$meth_top,$nonmeth_top),"\n"; + print GC join ("\t",$last_chr,$pos,$strand_top,$meth_top,$nonmeth_top,$context_top,$tri_nt_top),"\n"; + } + + } + + %chr = (); # resetting the hash + + # new first entry + $last_chr = $chr; + $chr{$chr}->{$start}->{meth} = $meth; + $chr{$chr}->{$start}->{nonmeth} = $nonmeth; + } + } + + # Last found chromosome + warn "Writing cytosine GpC report for last chromosome $last_chr (stored ",scalar keys %{$chr{$last_chr}}," different covered positions)\n"; + $processed{$last_chr} = 1; + while ( $chromosomes{$last_chr} =~ /(GC)/g){ + + # a GC on the top strand automatically means that there is a GC on the bottom strand as well, so we can process both at the same time + my $context_top = ''; + my $context_bottom = ''; + my $pos = pos$chromosomes{$last_chr}; + + my $meth_top = 0; + my $meth_bottom = 0; + my $nonmeth_top = 0; + my $nonmeth_bottom = 0; + + # C on forward strand + my $tri_nt_top = substr ($chromosomes{$last_chr},($pos-1),3); # positions are 0-based! + my $strand_top = '+'; + + # C on reverse strand + my $tri_nt_bottom = substr ($chromosomes{$last_chr},($pos-4),3); # positions are 0-based! + $tri_nt_bottom = reverse $tri_nt_bottom; + $tri_nt_bottom =~ tr/ACTG/TGAC/; + my $strand_bottom = '-'; + + next if (length $tri_nt_top < 3); # trinucleotide sequence could not be extracted for the top strand + next if (length $tri_nt_bottom < 3); # trinucleotide sequence could not be extracted for the bottom strand + + ### determining cytosine context top strand + if ($tri_nt_top =~ /^CG/){ + $context_top = 'CG'; + } + elsif ($tri_nt_top =~ /^C.{1}G$/){ + $context_top = 'CHG'; + } + elsif ($tri_nt_top =~ /^C.{2}$/){ + $context_top = 'CHH'; + } + else{ # if the context can't be determined the positions will not be printed + warn "The sequence context could not be determined for the top strand (found: '$tri_nt_top'). Skipping.\n"; + next; + } + + ### determining cytosine context bottom strand + if ($tri_nt_bottom =~ /^CG/){ + $context_bottom = 'CG'; + } + elsif ($tri_nt_bottom =~ /^C.{1}G$/){ + $context_bottom = 'CHG'; + } + elsif ($tri_nt_bottom =~ /^C.{2}$/){ + $context_bottom = 'CHH'; + } + else{ # if the context can't be determined the positions will not be printed + warn "The sequence context could not be determined for the bottom strand (found: '$tri_nt_bottom'). Skipping.\n"; + next; + } + + # top strand + if (exists $chr{$last_chr}->{($pos)}){ # stored positions are 1-based! + $meth_top = $chr{$last_chr}->{$pos}->{meth}; + $nonmeth_top = $chr{$last_chr}->{$pos}->{nonmeth}; + } + # bottom strand + if (exists $chr{$last_chr}->{($pos-1)}){ # stored positions are 1-based! -1 for bottom strand Cs + $meth_bottom = $chr{$last_chr}->{$pos-1}->{meth}; + $nonmeth_bottom = $chr{$last_chr}->{$pos-1}->{nonmeth}; + } + + # if Cs were not covered at all they are not written out + unless ($meth_bottom + $nonmeth_bottom == 0){ + my $percentage = sprintf ("%.6f",$meth_bottom/ ($meth_bottom + $nonmeth_bottom) *100); + print GCCOV join ("\t",$last_chr,$pos-1,$pos-1,$percentage,$meth_bottom,$nonmeth_bottom),"\n"; + print GC join ("\t",$last_chr,$pos-1,$strand_bottom,$meth_bottom,$nonmeth_bottom,$context_bottom,$tri_nt_bottom),"\n"; + } + unless ($meth_top + $nonmeth_top == 0){ + my $percentage = sprintf ("%.6f",$meth_top/ ($meth_top + $nonmeth_top) *100); + print GCCOV join ("\t",$last_chr,$pos,$pos,$percentage,$meth_top,$nonmeth_top),"\n"; + print GC join ("\t",$last_chr,$pos,$strand_top,$meth_top,$nonmeth_top,$context_top,$tri_nt_top),"\n"; + } + } + + warn "Finished writing out GpC cytosine report for covered chromosomes (processed $number_processed chromosomes/scaffolds in total)\n\n"; + + close GC or warn $!; + +} + +#### +#### + + +sub process_unprocessed_chromosomes{ + + my $chr = shift; + + warn "Writing cytosine report for not covered chromosome $chr\n"; + $processed{$chr} = 1; + + if ($split_by_chromosome){ ## writing output to 1 file per chromosome + my $chromosome_out = $cytosine_out; + if ($chromosome_out =~ /txt$/){ # files passed on by the methylation extractor end in _report.txt + $chromosome_out =~ s/txt$/chr${chr}.txt/; + } + else{ # user specified output file name + $chromosome_out =~ s/$/.chr${chr}.txt/; + } + + open (CYT,'>',$chromosome_out) or die $!; + # warn "Writing output for $last_chr to $chromosome_out\n"; + } + + my $tri_nt; + my $tetra_nt; + my $penta_nt; + my $context; + + while ( $chromosomes{$chr} =~ /([CG])/g){ + + $tri_nt = ''; + $context = ''; + if ($tetra){ + $tetra_nt = ''; # clearing + $penta_nt = ''; + } + + my $pos = pos$chromosomes{$chr}; + + my $strand; + my $meth = 0; + my $nonmeth = 0; + + if ($1 eq 'C'){ # C on forward strand + $tri_nt = substr ($chromosomes{$chr},($pos-1),3); # positions are 0-based! + $strand = '+'; + if ($tetra){ + if ( length($chromosomes{$chr}) >= ($pos - 1 + 4) ){ + $tetra_nt = substr ($chromosomes{$chr},($pos-1),4); + } + else{ + $tetra_nt = ''; + } + if ( length($chromosomes{$chr}) >= ($pos - 1 + 5) ){ + $penta_nt = substr ($chromosomes{$chr},($pos-1),5); + } + else{ + $penta_nt = ''; + } + } + } + elsif ($1 eq 'G'){ # C on reverse strand + $tri_nt = substr ($chromosomes{$chr},($pos-3),3); # positions are 0-based! + $tri_nt = reverse $tri_nt; + $tri_nt =~ tr/ACTG/TGAC/; + $strand = '-'; + + if ($tetra){ + if ( $pos - 4 >= 0 ){ + $tetra_nt = substr ($chromosomes{$chr},($pos-4),4); + $tetra_nt = reverse $tetra_nt; + $tetra_nt =~ tr/ACTG/TGAC/; + } + else{ + $tetra_nt = ''; + } + if ( $pos - 5 >= 0 ){ + $penta_nt = substr ($chromosomes{$chr},($pos-5),5); + $penta_nt = reverse $penta_nt; + $penta_nt =~ tr/ACTG/TGAC/; + } + else{ + $penta_nt = ''; + } + } + $strand = '-'; + } + + next if (length$tri_nt < 3); # trinucleotide sequence could not be extracted + + ### determining cytosine context + if ($tri_nt =~ /^CG/){ + $context = 'CG'; + } + elsif ($tri_nt =~ /^C.{1}G$/){ + $context = 'CHG'; + } + elsif ($tri_nt =~ /^C.{2}$/){ + $context = 'CHH'; + } + else{ # if the context can't be determined the positions will not be printed (it will equally not have been reported by Bismark) + warn "The cytosine context could not be determined (found: '$tri_nt'). Skipping.\n"; + next; + } + + if ($CpG_only){ + if ($tri_nt =~ /^CG/){ # CpG context is the default + if ($zero){ # zero-based coordinates + $pos -= 1; + if ($tetra){ + print CYT join ("\t",$chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt,$tetra_nt,$penta_nt),"\n"; + } + else{ + print CYT join ("\t",$chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n"; + } + } + else{ # default + if ($tetra){ + print CYT join ("\t",$chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt,$tetra_nt,$penta_nt),"\n"; + } + else{ + print CYT join ("\t",$chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n"; + } + } + } + } + else{ ## all cytosines, specified with --CX + if ($zero){ # zero based coordinates + $pos -= 1; + if ($tetra){ + print CYT join ("\t",$chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt,$tetra_nt,$penta_nt),"\n"; + } + else{ + print CYT join ("\t",$chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n"; + } + } + else{ # default + if ($tetra){ + print CYT join ("\t",$chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt,$tetra_nt,$penta_nt),"\n"; + } + else{ + print CYT join ("\t",$chr,$pos,$strand,$meth,$nonmeth,$context,$tri_nt),"\n"; + } + } + } + } + + # close CYT or warn $!; + +} + + +sub read_genome_into_memory{ + + ## reading in and storing the specified genome in the %chromosomes hash + chdir ($genome_folder) or die "Can't move to $genome_folder: $!"; + warn "Now reading in and storing sequence information of the genome specified in: $genome_folder\n\n"; + + my @chromosome_filenames = <*.fa>; + + ### if there aren't any genomic files with the extension .fa we will look for files with the extension .fasta + unless (@chromosome_filenames){ + @chromosome_filenames = <*.fasta>; + } + unless (@chromosome_filenames){ + die "The specified genome folder $genome_folder does not contain any sequence files in FastA format (with .fa or .fasta file extensions)\n"; + } + + foreach my $chromosome_filename (@chromosome_filenames){ + + # skipping the tophat entire mouse genome fasta file + next if ($chromosome_filename eq 'Mus_musculus.NCBIM37.fa'); + + open (CHR_IN,$chromosome_filename) or die "Failed to read from sequence file $chromosome_filename $!\n"; + ### first line needs to be a fastA header + my $first_line = <CHR_IN>; + chomp $first_line; + $first_line =~ s/\r//; # removing /r carriage returns + + ### Extracting chromosome name from the FastA header + my $chromosome_name = extract_chromosome_name($first_line); + + my $sequence; + while (<CHR_IN>){ + chomp; + $_ =~ s/\r//; # removing /r carriage returns + + if ($_ =~ /^>/){ + ### storing the previous chromosome in the %chromosomes hash, only relevant for Multi-Fasta-Files (MFA) + if (exists $chromosomes{$chromosome_name}){ + warn "chr $chromosome_name (",length $sequence ," bp)\n"; + die "Exiting because chromosome name already exists. Please make sure all chromosomes have a unique name!\n"; + } + else { + if (length($sequence) == 0){ + warn "Chromosome $chromosome_name in the multi-fasta file $chromosome_filename did not contain any sequence information!\n"; + } + warn "chr $chromosome_name (",length $sequence ," bp)\n"; + $chromosomes{$chromosome_name} = $sequence; + $processed{$chromosome_name} = 0; # processed chromosomes will be set to 1 later to allow a record of which chromosome has been processed + } + ### resetting the sequence variable + $sequence = ''; + ### setting new chromosome name + $chromosome_name = extract_chromosome_name($_); + } + else{ + $sequence .= uc$_; + } + } + + if (exists $chromosomes{$chromosome_name}){ + warn "chr $chromosome_name (",length $sequence ," bp)\t"; + die "Exiting because chromosome name already exists. Please make sure all chromosomes have a unique name.\n"; + } + else{ + if (length($sequence) == 0){ + warn "Chromosome $chromosome_name in the file $chromosome_filename did not contain any sequence information!\n"; + } + warn "chr $chromosome_name (",length $sequence ," bp)\n"; + $chromosomes{$chromosome_name} = $sequence; + $processed{$chromosome_name} = 0; # processed chromosomes will be set to 1 later to allow a record of which chromosome has been processed + } + } + warn "\n"; + chdir $parent_dir or die "Failed to move to directory $parent_dir\n"; +} + +sub extract_chromosome_name { + ## Bowtie extracts the first string after the inition > in the FASTA file, so we are doing this as well + my $fasta_header = shift; + if ($fasta_header =~ s/^>//){ + my ($chromosome_name) = split (/\s+/,$fasta_header); + return $chromosome_name; + } + else{ + die "The specified chromosome ($fasta_header) file doesn't seem to be in FASTA format as required!\n"; + } +} + + +sub print_helpfile{ + + warn <<EOF + + SYNOPSIS: + + This script generates a cytosine methylation report for a genome of interest and a sorted methylation input file produced + by the script "bismark2bedGraph". By default, the output uses 1-based chromosome coordinates and reports CpG positions only + (for both strands individually and not merged in any way). Coordinates may be changed to 0-based using the option '--zero_based'. + + The input file needs to have been generated with the script bismark2bedGraph (the file is called *.cov) or otherwise be + sorted by position and exactly in the following format: + + <chromosome> <start position> <end position> <methylation percentage> <count methylated> <count unmethylated> + + The coordinates of the input file are expected to be 1-based throughout (do not use files ending in .zero.cov!). + + + USAGE: coverage2cytosine [options] --genome_folder <path> -o <output> [input] + + +-o/--output <filename> Name of the output file, mandatory. + +--dir Output directory. Output is written to the current directory if not specified explicitly. + +--genome_folder <path> Enter the genome folder you wish to use to extract sequences from (full path only). Accepted + formats are FastA files ending with '.fa' or '.fasta'. Specifying a genome folder path is mandatory. + +-CX/--CX_context The output file contains information on every single cytosine in the genome irrespective of + its context. This applies to both forward and reverse strands. Please be aware that this will + generate output files with > 1.1 billion lines for a mammalian genome such as human or mouse. + Default: OFF (i.e. Default = CpG context only). + +--merge_CpG Using this option will post-process the genome-wide report to write out an additional coverage + file (see above for the coverage file format) which has the top and bottom strand methylation + evidence pooled into a single CpG dinucleotide entity. This may be the desirable input format + for some downstream processing tools such as the R-package bsseq (by K.D. Hansen). An example would be: + + genome-wide CpG report (old) + gi|9626372|ref|NC_001422.1| 157 + 313 156 CG + gi|9626372|ref|NC_001422.1| 158 - 335 156 CG + merged CpG evidence coverage file (new) + gi|9626372|ref|NC_001422.1| 157 158 67.500000 648 312 + + This option is currently experimental, and only works if CpG context only and a single genome-wide report + were specified (i.e. it doesn't work with the options --CX or --split_by_chromosome). + +--gc/--gc_context In addition to normal processing this option will reprocess the genome to find methylation in + GpC context. This might be useful for specialist applications where GpC methylases had been + deployed. The output format is exactly the same as for the normal cytosine report, and only + positions covered by at least one read are reported (output file ends in .GpC_report.txt). In addition + this will write out a Bismark coverage file (ending in GpC.cov). + +--ff In addition to the standard output selecting --ff will also extract a four and five (tetra/penta) + nucleotide context for the cytosines in question. Too short sequences (e.g. at the edges of the + chromosome) will be left blank; sequences containing Ns are ignored. + +--zero_based Uses 0-based coordinates instead of 1-based coordinates throughout. Default: OFF. + +--split_by_chromosome Writes the output into individual files for each chromosome instead of a single output file. Files + will be named to include the input filename and the chromosome number. + +--gzip Output file will be GZIP compressed (ending in .gz). Only works for standard CpG- and CX-output. + +--help Displays this help message and exits + + + +OUTPUT FORMAT: + +The genome-wide cytosine methylation output file is tab-delimited in the following format (1-based coords): +=========================================================================================================== + +<chromosome> <position> <strand> <count methylated> <count non-methylated> <C-context> <trinucleotide context> + + + Script last modified: 04 April 2016 + +EOF + ; + exit 1; +} +