#!/usr/bin/perl -w
## make easy usage of BisSNP
## author: Yaping Liu  lyping1986@gmail.com



use strict;
use Getopt::Long;




sub usage {

    print "\nUsage:\n";
    print "perl bissnp_easy_usage.pl [Options] BISSNP_JAR INPUT_BAM REF_GENOME DBSNP_VCF\n\n";

    print " Provide an easy way to operate BisSNP in one command line\n";
    print " It would automately generate: \n";
    print " 1) cpg.vcf (detailed CpGs methylation and genotyping information)\n";
    print " 2) snp.vcf (detailed SNPs information)\n";
    print " 3) cpg.bed (CpG methylation bed file) and cpg.coverage.bed (CpG C/T reads coverage bedGraph file)\n";
    print " 4) realigned and mark duplicated reads BAM file (if local indel realignment option is enabled)\n";
    print " 5) realigned, mark duplicated reads, base quality score recalibrated BAM file (if base quality score recalibration option is enabled)\n\n";

    
    print " All the BisSNP scripts should be put under the same directory. \n";
    print " It include the following steps:\n";
    print " 1). Indel realignment (Default: not enabled since its time consuming)\n";
    print " 2). Mark duplicated reads (Default: not enabled, if it is enabled, need to provide the path of MarkDuplicates.jar in picard tools)\n";
    print " 3). base quality recalibration (Default: not enabled since its time consuming)\n";
    print " 4). genotyping, get cpg.vcf and snp.vcf file\n";
    print " 5). filter out fake SNPs inside cpg.vcf and snp.vcf file\n";
    print " 6). generate cpg.methylation.bed(only homozygous CpGs, CpG adjacent in two strands are combined), \n";
    print " 	cpg.coverage.bedgraph\n\n";

	print "  BISSNP_JAR     Bis-SNP .jar name.\n";
    print "  INPUT_BAM     input BAM files to do genotyping. It should contain ReadGroup tag inside BAM.\n";
    print "  REF_GENOME    .fa/.fasta file, the reference genome file which should be the \n";
    print "              same when you do reads mapping to obtain BAM files.\n";
    print "  DBSNP_VCF     dbSNP.vcf file used for genotyping and base quality recalibration. Its chromosome order should be the same as REF_GENOME.\n\n";
	
	print "  [Options]:\n\n";
    print "  --nomeseq : enable NOMe-seq mode\n\n";
    print "  --lowCov :  enable low coverage mode (--qual default value will be set to 10)\n\n";
    print "  --indel :   enable indel local realignment\n\n";
    print "  --knownIndel INDEL_FILE : indel.vcf file name, when --indel enabled, this option need to be provided\n\n";
    print "  --recal :   enable base quality recalibration\n\n";
    print "  --duplicate DIR : directory of MarkDuplicates.jar in picard tools\n\n";
	print "  --qual NUM : genotyping quality score (Default: 20)\n\n";
	print "  --interval FILE : interval .bed file name\n\n";
	print "  --minCT NUM : minimum number of C or T reads to call methylation (Default: 1)\n\n";
	print "  --mmq NUM : minimum mapping quality score of reads used for genotyping (Default: 30)\n\n";
	print "  --mbq NUM : minimum base quality score of bases used for genotyping and methylation (Default: 5)\n\n";
	print "  --minConv NUM : minimum number of converted C in the begining of reads required\n";
	print "                  to begin count for genotyping and methylation (Default: 1)\n\n";
	print "  --nt NUM : number of CPU threads to use (Default: 1)\n\n";
	print "  --mem NUM : how much gigabytes memory to use (Default: 10)\n\n";
	print "  --allC : output all of cyosines in filtered VCF file into bed file, otherwise, it will only output CpG into bed file except for NOMe-seq mode\n";
	print "  --cytosines STRING : cytosine pattern to be analyzed (Default: --cytosine CG,1 --cytosine CH,1)\n";
	print "                      you can specify '--cytosine CG,1 --cytosine CHG,1 --cytosine CHH,1' to get CG, CHG and CHH information \n";
	print "                      NOMe-seq mode will change the default value to \n";
	print "                      '--cytosine GCH,2 --cytosine HCH,2 --cytosine HCG,2 --cytosine GCG,2'\n\n";
	print "  --outMode NUM : (Default: --outMode 1)\n";
	print "                 1) output CpG and SNPs\n";
	print "                 2) output Cytosine and SNPs\n";
	print "                 3) output Cytosines\n";
	print "                 4) output CpGs\n";
	print "                 5) output SNPs\n";
	print "                 6) output heterozygous SNPs\n";
	print "                 7) output all confident sites' genotype\n";
	print "                 8) output all callable sites' genotype\n\n";
    exit(1);
}

##default option setting
my $nomeseq = "";
my $lowCov = "";
my $indel = "";
my @knownIndels = ();
my $recal = "";
my $duplicate = "";
my $qual = 20;
my $interval = "";
my $minCT = 1;
my $mmq = 30;
my $mbq = 5;
my $minConv = 1;
my $nt = 1;
my $mem = 10;
my @cytosines = ();
my $outMode = 1;
my $allC = "";

GetOptions( "nomeseq" => \$nomeseq,
			"lowCov" => \$lowCov,
			"indel" => \$indel,
			"knownIndels=s" => \@knownIndels,
			"recal" => \$recal,
			"duplicate=s" => \$duplicate,
			"qual=i" => \$qual,
			"interval=s" => \$interval,
			"minCT=i" => \$minCT,
			"mmq=i" => \$mmq,
			"mbq=i" => \$mbq,
			"minConv=i" => \$minConv,
			"nt=i" => \$nt,
			"mem=i" => \$mem,
			"allC" => \$allC,
	    "cytosines=s" => \@cytosines,
	    "outMode=i" => \$outMode);

usage() if ( scalar(@ARGV) == 0 );

if ( scalar(@ARGV) != 4 ) {
    print "Wrong number of arguments\n";
    usage();
}

my $bissnp_jar = $ARGV[0];
my $input_file = $ARGV[1];
my $ref = $ARGV[2];
my $dbsnp = $ARGV[3];

my $PICARD = $duplicate;
my $BISSNP = $bissnp_jar;
my $JAVA = "java -Xmx$mem" . "G";

my $VCFTOOLS = "./sortByRefAndCor.pl";
my $VCF2BED = "./vcf2bed.pl";
my $VCF2COV = "./vcf2coverage.pl";
my $VCF2BEDGRAPH = "./vcf2bedGraph.pl";

unless (-e $BISSNP) {
 	print "\n $BISSNP jar file not exists!\n\n";
 	exit(1);
}

my %out_modes_hash = (
	1 => "DEFAULT_FOR_TCGA",
	2 => "EMIT_VARIANT_AND_CYTOSINES",
	3 => "EMIT_ALL_CYTOSINES",
	4 => "EMIT_ALL_CPG",
	5 => "EMIT_VARIANTS_ONLY",
	6 => "EMIT_HET_SNPS_ONLY",
	7 => "EMIT_ALL_CONFIDENT_SITES",
	8 => "EMIT_ALL_SITES",
);

if(scalar(@cytosines) == 0){
	@cytosines = ("CG,1", "CH,1");
}

if($allC ne ""){
		if($outMode == 1){
			$outMode = 2;
		}
		elsif($outMode == 4){
			$outMode = 3;
		}
}


if($nomeseq ne ""){
	@cytosines = ("GCH,2", "HCH,2","HCG,2", "GCG,2");
	$outMode = 2;
}

if($lowCov ne ""){
	if($qual == 20){
		$qual=10;
	}
}

## intermdediate file name 


my $tmp_dir="./";
my $original_input = $input_file;

if($indel ne ""){
	&bam_indel_realign();
}
if($duplicate ne ""){
	&bam_mdups();
}
if($recal ne ""){
	&bam_base_recalibration();
}

my $vcf_unsorted_cpg = $input_file;
my $vcf_unsorted_snp = $input_file;
$vcf_unsorted_cpg =~ s/\.bam//;
if($nomeseq ne "" or $outMode == 2 or $outMode == 3){
	$vcf_unsorted_cpg .= ".cytosine.raw.vcf";
}
elsif($outMode == 5 or $outMode == 6){
	$vcf_unsorted_cpg .= ".snp.raw.vcf";
}
elsif($outMode == 7 or $outMode == 8){
	$vcf_unsorted_cpg .= ".genotype.raw.vcf";
}
else{
	$vcf_unsorted_cpg .= ".cpg.raw.vcf";
}

$vcf_unsorted_snp =~ s/\.bam//;
$vcf_unsorted_snp .= ".snp.raw.vcf";

my $vcf_sorted_cpg = $vcf_unsorted_cpg;
my $vcf_sorted_snp = $vcf_unsorted_snp;
$vcf_sorted_cpg =~ s/.vcf//;
$vcf_sorted_cpg .= ".sort.vcf";
$vcf_sorted_snp =~ s/.vcf//;
$vcf_sorted_snp .= ".sort.vcf";

my $vcf_sorted_cpg_filtered = $vcf_sorted_cpg;
my $vcf_sorted_snp_filtered = $vcf_sorted_snp;
$vcf_sorted_cpg_filtered =~ s/.raw.sort.vcf//;
$vcf_sorted_cpg_filtered .= ".filtered.sort.vcf";
$vcf_sorted_snp_filtered =~ s/.raw.sort.vcf//;
$vcf_sorted_snp_filtered .= ".filtered.sort.vcf";

&bissnp();
&vcf_sort();
my $cmd="rm $vcf_unsorted_cpg";
system($cmd);
$cmd="rm $vcf_unsorted_snp";
system($cmd);
if($outMode != 3 and $outMode != 4 and $outMode != 7 and $outMode != 8 ){
	&vcf_filter();
}

&vcf2bed();



sub bam_indel_realign{
	if(scalar(@knownIndels) == 0){
		die "  Enabled indel realignemnt steps, but not provide known indels vcf file for the alignment!\n\n";
	}
	my $indel_target_interval = $input_file;
	$indel_target_interval =~ s/\.bam//;
	$indel_target_interval .= ".indels.intervals";
	my $header .= "$JAVA -jar $BISSNP -R $ref ";	
	$header .= "-I $input_file ";
	$header .= "-T BisulfiteRealignerTargetCreator ";
	if($interval ne ""){
		$header .= "-L $interval ";
	}
	foreach my $knownIndel(@knownIndels){
		$header .= "-known $knownIndel ";
	}
	
	$header .= "-o $indel_target_interval -nt $nt\n ";
										
	#realign 
	
	$header .= "$JAVA -jar $BISSNP -R $ref ";	
	$header .= "-I $input_file ";
	$header .= "-T BisulfiteIndelRealigner -targetIntervals $indel_target_interval ";
	foreach my $knownIndel(@knownIndels){
		$header .= "-known $knownIndel ";
	}
	if($lowCov ne ""){
		$header .= "-LOD 2.0 ";
	}
	$header .= "-compress 5 -cigar ";
	$input_file =~ s/\.bam//;
	$input_file .= ".realign.bam";
	$header .= "-o $input_file\n ";
	print "$header\n";
	system($header);
	
}

sub bam_mdups{
	my $remove_bam = $input_file;
	my $md = "$PICARD/MarkDuplicates.jar";
	unless (-e $md) {
 		print "\n MarkDuplicates.jar file not exists!\n\n";
 		exit(1);
	}
	my $header .= "$JAVA -jar $PICARD/MarkDuplicates.jar I=$input_file "; 
	$input_file =~ s/\.bam//;
	$input_file .= ".mdups.bam";
	$header .= "O=$input_file ";
	$header .= "CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT MAX_RECORDS_IN_RAM=3000000 ";
	my $metrics_file = $input_file.".dupsMetrics.txt";
	$header .= "METRICS_FILE=$metrics_file ";
	$header .= "TMP_DIR=$tmp_dir \n";
	print "$header\n";
	system($header);
	if($original_input ne $remove_bam){
		system("rm $remove_bam");
		$remove_bam =~ s/\.bam/.bai/;
		system("rm $remove_bam");
	}
	
	
}

sub bam_base_recalibration{


	my $remove_bam = $input_file;
	my $recalFile_before = $input_file;
	$recalFile_before =~ s/\.bam//;
	$recalFile_before .= ".beforeRecal.txt";
	my $recalFile_after=$input_file;
	$recalFile_after =~ s/\.bam//;
	$recalFile_after .= ".afterRecal.txt";
	
##1 countCovariant
	my $header .= "$JAVA -jar $BISSNP -R $ref ";	
	$header .= "-I $input_file ";
	$header .= "-T BisulfiteCountCovariates -nt $nt ";
	$header .= "-knownSites $dbsnp ";
	$header .= "-cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate ";
	$header .= "-recalFile $recalFile_before \n ";

##2 TableRecalibration
	$header .= "$JAVA -jar $BISSNP -R $ref ";	
	$header .= "-I $input_file ";
	$header .= "-T BisulfiteTableRecalibration ";
	$input_file =~ s/\.bam//;
	$input_file .= ".recal.bam";
	$header .= "-o $input_file ";
	$header .= "-recalFile $recalFile_before -maxQ 40 \n ";

##3 countCovariantAfterRecalibrate
#	$header .= "$JAVA -jar $BISSNP -R $ref ";	
#	$header .= "-I $input_bam_realign_mdups_recal ";
#	$header .= "-T BisulfiteCountCovariates -nt $nt ";
#	$header .= "-knownSites $dbsnp ";
#	$header .= "-cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate ";
#	$header .= "-recalFile $recalFile_after \n ";
	
	print "$header\n";
	system($header);
	if($original_input ne $remove_bam){
		system("rm $remove_bam");
		$remove_bam =~ s/\.bam/.bai/;
		system("rm $remove_bam");
	}
	
}

sub bissnp{
	my $header .= "$JAVA -jar $BISSNP -R $ref ";
	$header .= "-I $input_file ";									
	$header .= "-D $dbsnp -T BisulfiteGenotyper -vfn1 $vcf_unsorted_cpg -vfn2 $vcf_unsorted_snp ";
	if($nomeseq ne ""){
		$header .= "-out_modes $out_modes_hash{ $outMode } -sm GM ";
	}
	else{
		foreach my $cytosine(@cytosines){
			$header .= "-C $cytosine ";
		}
		
		$header .= "-out_modes $out_modes_hash{ $outMode } ";
	}
	if($interval ne ""){
		$header .= "-L $interval ";
	}
	$header .= "-stand_call_conf $qual -stand_emit_conf 0 -nt $nt ";
	$header .= "-minConv $minConv -vcfCache 1000000 ";
	$header .= "-mmq $mmq ";
	$header .= "-mbq $mbq\n";
	print "$header\n";
	system($header);

}
 
sub vcf_sort{
	unless (-e $VCFTOOLS) {
 		print "\n $VCFTOOLS file not exists!\n\n";
 		exit(1);
	}
	my $header .= "perl $VCFTOOLS --k 1 --c 2 ";
	$header .= "--tmp $tmp_dir ";
	$header .= "$vcf_unsorted_cpg ";
	$header .= "$ref";
	$header .= ".fai ";
	$header .= "> $vcf_sorted_cpg \n";
	$header .= "perl $VCFTOOLS --k 1 --c 2 ";
	$header .= "--tmp $tmp_dir ";
	$header .= "$vcf_unsorted_snp ";
	$header .= "$ref";
	$header .= ".fai ";
	$header .= "> $vcf_sorted_snp \n";
	print "$header\n";
	system($header);
	
}

sub vcf_filter{
	my $header .= "$JAVA -jar $BISSNP -R $ref -T VCFpostprocess ";
	$header .= "-qual $qual ";
	if($interval ne ""){
		$header .= "-L $interval ";
	}
	if($nomeseq ne ""){
		$header .= "-C GCH -C HCH -C GCG -C HCG ";
	}
	else{
		foreach my $cytosine(@cytosines){
			$cytosine =~ s/\,\d+//;
			$header .= "-C $cytosine ";
		}
	}
	if($lowCov ne ""){
		$header .= "-sb 0.0 ";
	}
	$header .= "-oldVcf $vcf_sorted_cpg ";
	$header .= "-snpVcf $vcf_sorted_snp ";
	$header .= "-newVcf $vcf_sorted_cpg_filtered ";
	$header .= "-o $vcf_sorted_cpg_filtered.cpgSummary.txt -minCT $minCT\n";
	
	$header .= "$JAVA -jar $BISSNP -R $ref -T VCFpostprocess ";
	$header .= "-qual $qual ";
	if($interval ne ""){
		$header .= "-L $interval ";
	}
	if($nomeseq ne ""){
		$header .= "-C GCH -C HCH -C GCG -C HCG ";
	}
	else{
		foreach my $cytosine(@cytosines){
			$cytosine =~ s/\,\d+//;
			$header .= "-C $cytosine ";
		}
	}
	if($lowCov ne ""){
		$header .= "-sb 0.0 ";
	}
	$header .= "-oldVcf $vcf_sorted_snp ";
	$header .= "-snpVcf $vcf_sorted_snp ";
	$header .= "-newVcf $vcf_sorted_snp_filtered ";
	$header .= "-o $vcf_sorted_snp_filtered.cpgSummary.txt \n";
	print "$header\n";
	system($header);
	
	
}

sub vcf2bed{
	unless (-e $VCF2BED) {
 		print "\n $VCF2BED file not exists! \n\n";
 		exit(1);
	}
	my $header .= "perl $VCF2BED $vcf_sorted_cpg_filtered ";
	if($nomeseq ne ""){
		$header .= "GCH \n";
		$header .= "perl $VCF2BED $vcf_sorted_cpg_filtered HCG \n";
	}
	else{
		if($allC eq ""){
			$header .= "CG \n";
		}
		else{
			$header = "";
			foreach my $cytosine(@cytosines){
				$cytosine =~ s/\,\d+//;
				$header .= "perl $VCF2BED $vcf_sorted_cpg_filtered $cytosine \n";
			}
			
		}
			
	}
	print "$header\n";
	system($header);
	
	unless (-e $VCF2COV) {
 		print "\n $VCF2COV file not exists! \n\n";
 		exit(1);
	}
	$header = "perl $VCF2COV $vcf_sorted_cpg_filtered ";
	if($nomeseq ne ""){
		$header .= "GCH \n";
		$header .= "perl $VCF2COV $vcf_sorted_cpg_filtered HCG \n";
	}
	else{
		if($allC eq ""){
			$header .= "CG \n";
		}
		else{
			$header = "";
			foreach my $cytosine(@cytosines){
				$cytosine =~ s/\,\d+//;
				$header .= "perl $VCF2COV $vcf_sorted_cpg_filtered $cytosine \n";
			}
			
		}
	}
	
	print "$header\n";
	system($header);
	
	unless (-e $VCF2BEDGRAPH) {
 		print "\n $VCF2BEDGRAPH file not exists! \n\n";
 		exit(1);
	}
	$header = "perl $VCF2BEDGRAPH $vcf_sorted_cpg_filtered ";
	if($nomeseq ne ""){
		$header .= "GCH \n";
		$header .= "perl $VCF2BEDGRAPH $vcf_sorted_cpg_filtered HCG \n";
	}
	else{
		if($allC eq ""){
			$header .= "CG \n";
		}
		else{
			$header = "";
			foreach my $cytosine(@cytosines){
				$cytosine =~ s/\,\d+//;
				$header .= "perl $VCF2BEDGRAPH $vcf_sorted_cpg_filtered $cytosine \n";
			}
			
		}
			
	}
	print "$header\n";
	system($header);
	
}

