text2ui

#!/usr/bin/env perl

# Transforms the text annotations in Pokay's data dir into HTML pages and a NextClade JSON file for linking to them.

use strict;
use warnings;
use File::Copy;
use Time::Piece;

my $max_histogram_x = 9;

@ARGV == 4 or @ARGV == 6 or die "Usage: $0 <constellations of interest.tab> <input dir with txt files> <output html dir> <output.json> [b_epitopes.fasta t_epitopes.fasta]\n";

my $vocs_tab_file = shift @ARGV;
my $text_indir = shift @ARGV;
my $html_outdir = shift @ARGV;
my $json_outfile = shift @ARGV;
my $b_cell_epitope_fasta = @ARGV ? shift @ARGV : undef;
my $t_cell_epitope_fasta = @ARGV ? shift @ARGV : undef;

if(-e $html_outdir){
  die "Specified HTML output directory '$html_outdir' exists, but is not a directory\n" if not -d $html_outdir;
}
elsif(not mkdir $html_outdir, 0755){
  die "Specified HTML output directory '$html_outdir' could not be created: $!\n";
}

open(JSON_OUT, ">$json_outfile")
  or die "Cannot open $json_outfile for writing: $!\n";

# We may for convenience want to lump together all the variants for a lineage into one report from the functional perspective
# {name => {gene => [mutations]}}
my %lineage_mutations;
my %lineage_optional_mutations;
my %lineage_nickname;
my %lineage2epi_id;
open(TAB_FILE, $vocs_tab_file)
  or die "Cannot open $vocs_tab_file for reading: $!\n";
while(<TAB_FILE>){
	next if /^\s*#/; # comment lines
	chomp;
	my @fields = split /\t/, $_;
	if(@fields > 5){
		warn "Skipping line $. of $vocs_tab_file because it has more tab delimited fields (", scalar(@fields), ") than expected [2,5]\n";
		next;
	}
	my $lineage_name = $fields[0];
	my %gene_to_mutations;
	for my $mutation_spec (split /\s*;\s*/, $fields[1]){
		my ($gene, $mutation) = split /:p.|:/, $mutation_spec;
		$lineage_mutations{$lineage_name}->{$gene} //= [];
		push @{$lineage_mutations{$lineage_name}->{$gene}}, $mutation;
	}
	next if $#fields == 1; # no optional mutations in this lineage
        for my $mutation_spec (split /\s*;\s*/, $fields[2]){
                my ($gene, $mutation) = split /:p.|:/, $mutation_spec;
                $lineage_optional_mutations{$lineage_name}->{$gene} //= [];
                push @{$lineage_optional_mutations{$lineage_name}->{$gene}}, $mutation;
        }
	# This column is nick names for the lineage, like Alpha
	if($#fields >= 3){
		$lineage_nickname{$lineage_name} = $fields[3];
	}
	# This column is a GISAID genome representative of the lineage for epitope mismatch calculation purposes (a corresponding FastA file for it must exist in the data dir for that to work)
	if($#fields == 4){
		$lineage2epi_id{$lineage_name} = $fields[4];
	}

}
close(TAB_FILE);

my %constellations2categories; # { unique var combo name => { short effect desc => [lit references] } }
opendir(TXT_INDIR, $text_indir)
  or die "Cannot open directory $text_indir for reading: $!\n";
while(my $txt_file = readdir(TXT_INDIR)){
  next if $txt_file =~ /^\./; # hidden files
  next unless -f "$text_indir/$txt_file" and $txt_file =~ /^(\S+?)_(.+)\.txt$/; # bona fide text file with gene_short_effect_desc.txt naming convention
  my ($gene, $short_effect_desc) = ($1, $2);

  my $qualified_txt_file = "$text_indir/$txt_file";
  open(TXT_FILE, $qualified_txt_file)
    or die "Cannot open $qualified_txt_file for reading: $!\n";
  #print "Processing $qualified_txt_file\n";

  my @lit_references;
  my $lit_link_evidence_summary;
  while(<TXT_FILE>){
    if(/^[A-Z]+\d+(?:[A-Z]+|del|_[A-Z]\d+)/){
      tr/ //d; # remove any errant whitespace
      chomp $_;
      # Sort the constellation key atoms so the naming is consistent even if the file encoding of the combos was not
      $_ = join(";", sort {&constellation_cmp($a,$b)} split(/;/, $_));
      my $constellation_key = "$gene:$_";
      if(not @lit_references){
        warn "Ignoring variant '$&' without literature annotations ($qualified_txt_file line #$.)\n";
      }
      $constellations2categories{$constellation_key}->{$short_effect_desc} //= [];
      push @{$constellations2categories{$constellation_key}->{$short_effect_desc}}, @lit_references;
      @lit_references = ();
    }
    elsif(/^\s*#\s*(.*)/){
      $_ = $1;
      # Lit link
      if(m(^(.+?)\s+(https?://\S+|unpublished))){
        my ($authors, $link) = ($1, $2);
        push @lit_references, [$lit_link_evidence_summary, $authors, $link];
        $lit_link_evidence_summary = "";
      }
      # Otherwise assume part of summary evidence text for lit link to come.
      # Maybe a missing period between lines of evidence from same paper.
      elsif(/^[A-Z](?=[a-z])/ and defined $lit_link_evidence_summary and $lit_link_evidence_summary ne "" and $lit_link_evidence_summary !~ /\.\s*$/){
        $lit_link_evidence_summary .= ".<br/>".$_;
      }
      elsif(/\.\s*$/){ # Keep the line ending
        $lit_link_evidence_summary .= " ".$_."<br/>";
      }
      else{ # assume it's a sentence that runs over two lines
        $lit_link_evidence_summary .= " ".$_;
      }
      $lit_link_evidence_summary =~ s/\.\s*<br\/>\s*\.<br\/>/.<br\/>/g;
    }
    elsif(not /^\s*$/){
      warn "Ignoring line that is not a comment starting with '#', nor a properly formatted variant ($qualified_txt_file line #$.)\n";
      warn "Also dumping preceding comments.\n" if $lit_link_evidence_summary;
      $lit_link_evidence_summary = "";
    } # else it's blank
  }
  close(TXT_FILE);
}
close(TXT_INDIR);

my @json_entries;
my %constellation2html;
# Store summary data so we can generate summary Web pages per gene and evidence type afterwards.
my %category2gene2constellations;
my %gene2category2constellations;

for my $constellation_key (keys %constellations2categories){
  my ($gene, $constellation) = split /:/, $constellation_key;
  $gene2category2constellations{$gene} ||= {};
  my @substitutions;
  my @deletions;
  while($constellation =~ /([A-Z])(\d+)([A-Z])/g){
        push @substitutions, "{ \"refAA\": \"$1\", \"queryAA\": \"$3\", \"codon\": $2, \"gene\": \"$gene\" }";
  }
  while($constellation =~ /([A-Z])(\d+)del/g){
        push @deletions, "{ \"refAA\": \"$1\", \"codon\": $2, \"$gene\": \"$gene\"}";
  }
  my $substitutions = join(",\n", @substitutions);
  my $deletions = join(",\n", @deletions);

  $constellation2html{$constellation_key} = ();

  for my $short_effect_desc (sort keys %{$constellations2categories{$constellation_key}}){
    $category2gene2constellations{$short_effect_desc} //= {};
    $category2gene2constellations{$short_effect_desc}->{$gene} //= [];
    push @{$category2gene2constellations{$short_effect_desc}->{$gene}}, $constellation_key;
    $gene2category2constellations{$gene} //= {};
    $gene2category2constellations{$gene}->{$short_effect_desc} //= [];
    push @{$gene2category2constellations{$gene}->{$short_effect_desc}}, $constellation_key;

    my @summaries;
    for my $lit_ref (@{$constellations2categories{$constellation_key}->{$short_effect_desc}}) {
      push @summaries, textify_lit_ref($lit_ref); # Call fuunction to  get summary
    }
    my $text_summary = join("; ", @summaries); # Join all summaries
    
    push @json_entries, <<ENTRY;
  {
    "description": "$short_effect_desc",
    "url": "http://people.ucalgary.ca/~gordonp/$constellation_key.html#$short_effect_desc",
    "summary": "$text_summary",
    "substitutions": 
      [
        $substitutions
      ],
    "deletions": 
      [
        $deletions
      ]
  }
ENTRY

    my $short_effect_desc_pretty = prettify_effect_text($short_effect_desc);
    push @{$constellation2html{$constellation_key}}, "<div class=\"effect_section\"><h2 id=\"$short_effect_desc\">$short_effect_desc_pretty</h2>\n<ul class=\"lit_list\">\n";
    
    for my $lit_ref (@{$constellations2categories{$constellation_key}->{$short_effect_desc}}){
      push @{$constellation2html{$constellation_key}}, htmlify_lit_ref($lit_ref);
    }
    push @{$constellation2html{$constellation_key}}, "</ul></div>\n";
  }
  
}
print JSON_OUT "[\n", join(",\n", @json_entries), "]\n";
close(JSON_OUT);

# Using array for the front page ordering later
my @gene_pretty = ( nsp1 => "Host mRNA translation inhibiting [40S subunit binding] protein",
                    nsp2 => "Prohibitin-interacting non-structural protein",
                    PLpro => "Papain-like protease / ER lumenal double-membrane vesicle forming complex member (a.k.a. NSP3)",
                    nsp4 => "ER lumenal double-membrane vesicle forming complex member",
                    "3CL" => "3C-like protease (a.k.a. NSP5, Mpro, main protease)",
                    nsp6 => "ER lumenal double-membrane vesicle forming complex member / IFN-I signalling suppressor",
                    nsp7 => "Hexadecameric supercomplex member of the core replicase complex",
                    nsp8 => "Adenylyltransferase / hexadecameric supercomplex member of the core replicase complex",
                    nsp9 => "ssRNA-binding and orienting member of the replicase complex",
                    nsp10 => "2'-O-ribose type 1 cap methyltransferase cofactor",
                    nsp11 => "Intrinsically disordered protein",
                    RdRp => "RNA-dependent RNA polymerase member protein of the core replicase complex (a.k.a. NSP12)",
                    Hel => "Helicase / packaging signal",
                    ExoN => "Exonuclease / guanine N7 methyltransferase / packaging signal (a.k.a. NSP14)",
                    EndoU => "5' subgenomic antisense transcript uracil endonuclease (a.k.a. NSP15)",
                    nsp16 => "2'-O-ribose type 1 cap methyltransferase",
                    S => "Spike protein",
                    ORF3a => "Golgi-associating, caveolin-binding protein",
                    E => "Envelope protein",
                    M => "Membrane protein",
                    ORF6 => "Nuclear import / interferon response inhibiting protein",
                    ORF7a => "Cell cycle / program interferring protein",
                    ORF7b => "Non-essential transmembrane protein",
                    ORF8 => "Ig-like folding / IFN-I response interferring protein",
                    N => "Nucleocapsid"
                    );

my %gene_pretty = @gene_pretty;

my $front_rel_url = "front.html";
# Write the per-constellation pages
for my $constellation_key (keys %constellation2html){
  my ($plain_gene, $constellation) = split /:/, $constellation_key;
 my $gene = $plain_gene;
  if(exists $gene_pretty{$gene}){
    $gene = $gene_pretty{$gene}. " ($gene)";
  }
  my $constellation_key_pretty = hgvsify($constellation_key);

  open(HTML, ">$html_outdir/$constellation_key.html")
    or die "Cannot open $html_outdir/$constellation_key.html for writing: $!\n";
  print HTML <<DATA;
<html>
<head>
  <title>SARS-CoV-2 $gene gene $constellation literature reference collection</title>
  <link rel="stylesheet" href="constellations.css">
</head>
<body>
<h1><a href="$front_rel_url">SARS-CoV-2</a> <a href="$plain_gene.html">$gene</a> gene <span class="constellation">$constellation_key_pretty</span> literature reference collection</h1>
DATA
  print HTML @{$constellation2html{$constellation_key}};
  print HTML <<DATA;
<div class="acknowledgement">Automatically generated using text2ui script and data from <a href="https://github.com/nodrogluap/pokay">Pokay</a></div>
</body>
</html>
DATA
  close(HTML);
}

open(MAP, "$text_indir/../imagemap.txt")
  or die "Cannot open $text_indir/../imagemap.txt for reading: $!\n";
my $map_data = join("", <MAP>);
close(MAP);

open(FRONT_HTML, ">$html_outdir/$front_rel_url")
  or die "Cannot open $html_outdir/$front_rel_url for writing: $!\n";
print FRONT_HTML <<END;
<html>
<head>
  <meta name="viewport" content="width=device-width, initial-scale=1">
  <title>SARS-CoV-2 mutation functional literature references</title>
  <link rel="stylesheet" href="constellations.css">
  <link rel="stylesheet" href="https://code.jquery.com/ui/1.12.1/themes/base/jquery-ui.css">
  <script src="https://code.jquery.com/jquery-1.12.4.js"></script>
  <script src="https://code.jquery.com/ui/1.12.1/jquery-ui.js"></script>
  <script src="imageMapResizer.min.js"></script>
  <script>
  \$( function() {
    \$( "#gene_accordion" ).accordion({
      collapsible: true,
      heightStyle: "content"
    });
  } );
 \$( function() {
    \$( "#type_accordion" ).accordion({
      collapsible: true,
      heightStyle: "content"
    });
  } );
 \$(document).ready(function() {
    \$('map').imageMapResize();
  } );

  </script>
</head>
<body>
<h1>SARS-CoV-2 mutation functional literature references</h1>
<div class="tip">Have a sequence to analyze for functional mutations?<br/>Perform quality control, clade assignment, and link to the mutation constellation literature references on this site using our fork of<br/><span style="font-size: 150%"><a href="./nextclade/" target="_blank"><span style="color:rgb(67, 119, 205)">N</span><span style="color:rgb(80, 151, 186);">e</span><span style="color:rgb(99, 172, 154);">x</span><span style="color:rgb(124, 184, 121);">t</span><span style="color:rgb(185, 188, 74);">c</span><span style="color:rgb(212, 177, 63);">l</span><span style="color:rgb(228, 153, 56);">a</span><span style="color:rgb(230, 112, 48);">d</span><span style="color:rgb(222, 60, 38);">e</span></a></span></div>
<img class="genome_schematic" src="genome.png" width="100%" alt="SARS-CoV-2 genome schematic" usemap="#map" /> 
<map name="map" id="map">
$map_data
</map>
<hr/>

<table width="100%">
<tbody>
<tr><td align="center"><em>Functional literature for Variants of Concern/Interest/Under Monitoring lineages</em></td></tr>
<tr><td align="center"> 
END

my @lineage_names_html;
for my $lineage_name (sort keys %lineage_mutations){
	push @lineage_names_html, "<a href=\"lineage-$lineage_name.html\">$lineage_name". 
             (exists $lineage_nickname{$lineage_name} ? " (".$lineage_nickname{$lineage_name}.(exists $lineage2epi_id{$lineage_name}?"*":"").")" : ""). 
             "</a>";
}
print FRONT_HTML join(" &mdash;\n", @lineage_names_html);

print FRONT_HTML <<END;
</td></tr>
</tbody>
</table>
<hr/>
<table>
<tbody>
<tr><td>
<h2>Functional literature by gene</h2>
<div id="gene_accordion">
END

# Printing the genes in genomic order
for (my $i = 0; $i <= $#gene_pretty; $i+=2){
  my $gene = $gene_pretty[$i];
  my $gene_desc = $gene_pretty[$i+1];

  my $total_num_constellations_this_gene = 0;
  my $gene_rel_url = "$gene.html";
  open(GENE_HTML, ">$html_outdir/$gene_rel_url")
    or die "Cannot open $html_outdir/$gene_rel_url for writing: $!\n";
  print GENE_HTML "<html><head><link rel=\"stylesheet\" href=\"constellations.css\"></head><body><h1>Functional mutation literature for <a href=\"$front_rel_url\">SARS-CoV-2</a> $gene &mdash;<br/>$gene_desc</h1><ul class=\"category_list\">\n";

  print FRONT_HTML "<h3>$gene &mdash; $gene_desc</h3>\n<div><ul class=\"category_list\">\n";
  for my $category (sort keys %{$gene2category2constellations{$gene}}){
    my $category_pretty = $category;
    $category_pretty =~ tr/_/ /;

    open(GENE_CAT_HTML, ">$html_outdir/$gene-$category.html")
      or die "Cannot open $html_outdir/$gene-$category.html for writing: $!\n";
    print GENE_CAT_HTML <<END;
<html>
<head>
  <title>SARS-CoV-2 gene $gene literature references for $category_pretty</title>
  <link rel="stylesheet" href="constellations.css">
</head>
<body>
<h1><a href="$front_rel_url">SARS-CoV-2</a> gene <a href="$gene_rel_url">$gene</a> literature references for $category_pretty</h1>
<ul class=\"constellation_list\">
END

    my $num_constellations = 0;
    $category_pretty =~ s/^([a-z])/uc($1)/e;
    for my $constellation_key (sort @{$gene2category2constellations{$gene}->{$category}}){
      my $simplified_constellation_key = hgvsify($constellation_key);
      $simplified_constellation_key =~ s/(\S{1,24};)/$1 /g;
      print GENE_CAT_HTML "  <li class=\"constellation_list_item\"><a href=\"./$constellation_key.html#$category\"><span class=\"constellation\">$simplified_constellation_key</span></a></li>\n";
      $num_constellations++;
    }

    print GENE_CAT_HTML <<END;
</ul>
<div class="acknowledgement">Automatically generated using text2ui script and data from <a href="https://github.com/nodrogluap/pokay">Pokay</a></div>
</body>
</html>
END

    close(GENE_CAT_HTML);
    $total_num_constellations_this_gene += $num_constellations;

    print FRONT_HTML "  <li class=\"category_list_item\"><a href=\"$gene-$category.html\" class=\"gene_category_page\">$category_pretty &mdash; ", 
                     $num_constellations, " mutation constellation", ($num_constellations > 1 ? "s" : ""), "</a></li>\n";
    print GENE_HTML "  <li class=\"category_list_item\"><a href=\"$gene-$category.html\" class=\"gene_category_page\">$category_pretty &mdash; ", 
                     $num_constellations, " mutation constellation", ($num_constellations > 1 ? "s" : ""), "</a></li>\n";
  }
  print FRONT_HTML "</ul></div>\n";

  print GENE_HTML "</ul>\n";
  if($total_num_constellations_this_gene == 0){
    print GENE_HTML "<div class=\"tip\">Sorry, there are currently no functionally characterized mutations in the literature for this gene</div>\n";
  }
  print GENE_HTML "<div class=\"acknowledgement\">Automatically generated using text2ui script and data from <a href=\"https://github.com/nodrogluap/pokay\">Pokay</a></div></body></html>\n";
  close(GENE_HTML);
}

print FRONT_HTML "</div></td><td>\n<h2>Functional literature by evidence type</h2>\n<div id=\"type_accordion\">\n";

# Print the categories alphabetically
for my $category (sort keys %category2gene2constellations){
  my $category_pretty = $category;
  $category_pretty =~ s/(^[a-z])/uc($1)/e;
  $category_pretty =~ tr/_/ /;

  print FRONT_HTML "<h3>$category_pretty</h3><div><ul class=\"gene_list\">";

  for my $gene (sort keys %{$category2gene2constellations{$category}}){

     open(CAT_GENE_HTML, ">$html_outdir/$category-$gene.html")
       or die "Cannot open $html_outdir/$category-$gene.html for writing: $!\n";
     print CAT_GENE_HTML <<END;
<html>
<head>
  <title>SARS-CoV-2 $category_pretty literature references for gene $gene</title>
  <link rel="stylesheet" href="constellations.css">
</head>
<body>
<h1><a href="$front_rel_url">SARS-CoV-2</a> $category_pretty literature references for gene <a href="$gene.html">$gene</a></h1>
<ul class=\"constellation_list\">
END

    my $num_constellations = 0;
    for my $constellation_key (sort @{$gene2category2constellations{$gene}->{$category}}){
      my $simplified_constellation_key = hgvsify($constellation_key);
      # Some constellation keys are getting pretty long, let's insert breakable space in them so we don't get overlapping text in the columns list layout used
      $simplified_constellation_key =~ s/(\S{1,24};)/$1 /g;
      print CAT_GENE_HTML "  <li class=\"constellation_list_item\"><a href=\"./$constellation_key.html#$category\"><span class=\"constellation\">$simplified_constellation_key</span></a></li>\n";
      $num_constellations++;
    }

   print CAT_GENE_HTML <<END;
</ul>
<div class="acknowledgement">Automatically generated using text2ui script and data from <a href="https://github.com/nodrogluap/pokay">Pokay</a></div>
</body>
</html>
END

     close(CAT_GENE_HTML);

     print FRONT_HTML "  <li class=\"gene_list_item\"><a href=\"$category-$gene.html\" class=\"category_gene_page\">$gene &mdash; ",
                      $num_constellations, " mutation constellation", ($num_constellations > 1 ? "s" : ""), "</a></li>\n";
  }
  print FRONT_HTML "</ul></div>\n";
}


my $date = `date`;
print FRONT_HTML <<END;
</div></td></tr></tbody></table>

<div class="footer">
<div class="acknowledgement">Automatically generated using text2ui script and data from <a href="https://github.com/nodrogluap/pokay">Pokay</a></div>
This site is maintained by Paul Gordon &lt;gordonp\@ucalgary.ca&gt;. Added references are welcome by submitting data directory pull requests on <a href="https://github.com/nodrogluap/pokay/tree/master/data">GitHub</a>. Updated $date.
</div>
</body>
</html>
END
close(FRONT_HTML);


# Write all the lineage-specific pages
for my $lineage_name (keys %lineage_mutations){
	my $lineage_short_name = $lineage_name;
	$lineage_short_name =~ s/\+.*//; # get rid of B.1+E484K addditional variant info that cov-lineages.org won't handle

	# Calculate the mismatch data to graph for each lineage if a sequence was provided
	my $generate_histograms = 0;
	my $epi_id = $lineage2epi_id{$lineage_name};
	my %b_mismatch_count;
	my $b_query_total = 0;
	my $epi_fasta_file = defined $epi_id ? "$text_indir/$epi_id.fasta" : undef;
	my $epi_date;
	my $iedb_date;
	if(defined $epi_id and defined $b_cell_epitope_fasta){
		if(not -e $epi_fasta_file){
			warn "Skipping epitope mismatch histograms generation for lineage $lineage_name, specified sequence ID $epi_id (in $vocs_tab_file) does not have a corresponding $epi_fasta_file\n";
		}
		elsif(not -e "$text_indir/$lineage_name.acknowledgements.pdf"){
			warn "Skipping epitope mismatch histograms generation for lineage $lineage_name, FastA file $epi_fasta_file exists, but no sequence acknowledgment PDF exists (expected at $text_indir/$lineage_name.acknowledgements.pdf)\n";
		}
		else{
 			my $time = localtime() - (-M $epi_fasta_file)*60*60*24;
			$epi_date = $time->date;
			$time = localtime() - (-M $b_cell_epitope_fasta)*60*60*24;
			$iedb_date = $time->date;
			system("./pokay $epi_fasta_file /dev/null $b_cell_epitope_fasta $html_outdir/$lineage_name.b.html $html_outdir/$lineage_name.b.histogram.txt") >> 8 
			  and die "Running Pokay epitope mismatch alignments for $epi_fasta_file exited abnormally: exit code was ".($? >> 8)."\n";
			open(HIST, "$html_outdir/$lineage_name.b.histogram.txt")
			  or die "Cannot open expected Pokay epitope search results file $html_outdir/$lineage_name.b.histogram.txt: $!\n";
			my $hist_header = <HIST>;
			($b_query_total) = $hist_header =~ /(\d+)/;
			while(<HIST>){
				chomp;
				my @F = split /\t/, $_;
				$b_mismatch_count{$F[0]} = 100*$F[1];
			}
			close(HIST);
			system("cp $text_indir/$lineage_name.acknowledgements.pdf $html_outdir");
		}
	}
	my %t_mismatch_count;
	my $t_query_total = 0;
        if(defined $epi_id and defined $t_cell_epitope_fasta){
                if(not -e $epi_fasta_file){
                        warn "Skipping epitope mismatch histograms generation for lineage $lineage_name, specified sequence ID $epi_id (in $vocs_tab_file) does not have a corresponding $epi_fasta_file\n";
                }
		elsif(not -e "$text_indir/$lineage_name.acknowledgements.pdf"){
			warn "Skipping epitope mismatch histograms generation for lineage $lineage_name, FastA file $epi_fasta_file exists, but no sequence acknowledgment PDF exists (expected at $text_indir/$lineage_name.acknowledgements.pdf)\n";
		}
                else{   
 			my $time = localtime() - (-M $epi_fasta_file)*60*60*24;
			$epi_date = $time->date;
			$time = localtime() - (-M $t_cell_epitope_fasta)*60*60*24;
			$iedb_date = $time->date;
                        system("./pokay $epi_fasta_file /dev/null $t_cell_epitope_fasta $html_outdir/$lineage_name.t.html $html_outdir/$lineage_name.t.histogram.txt") >> 8  
                          and die "Running Pokay epitope mismatch alignments for $epi_fasta_file exited abnormally: exit code was ".($? >> 8)."\n";
			open(HIST, "$html_outdir/$lineage_name.t.histogram.txt")
			  or die "Cannot open expected Pokay epitope search results file $html_outdir/$lineage_name.t.histogram.txt: $!\n";
			my $hist_header = <HIST>;
			($t_query_total) = $hist_header =~ /(\d+)/;
			while(<HIST>){
				chomp;
				my @F = split /\t/, $_;
				$t_mismatch_count{$F[0]} = 100*$F[1];
			}
			close(HIST);
			system("cp $text_indir/$lineage_name.acknowledgements.pdf $html_outdir");
                }
        }

	open(LINEAGE_HTML, ">$html_outdir/lineage-$lineage_name.html")
		or die "Cannot open $html_outdir/lineage-$lineage_name.html for writing: $!\n";
	print LINEAGE_HTML <<END;
<html>
<head>
  <meta name="viewport" content="width=device-width, initial-scale=1">
  <title>Functional mutation literature for SARS-CoV-2 lineage $lineage_name</title>
  <link rel="stylesheet" href="constellations.css">
  <link rel="stylesheet" href="https://code.jquery.com/ui/1.12.1/themes/base/jquery-ui.css">
  <script src="https://code.jquery.com/jquery-1.12.4.js"></script>
  <script src="https://code.jquery.com/ui/1.12.1/jquery-ui.js"></script>
  <script src="https://canvasjs.com/assets/script/canvasjs.min.js"></script>
  <script src="imageMapResizer.min.js"></script>
  <script>
 \$( function() {
    \$( "#accordion" ).accordion({
      collapsible: true,
      heightStyle: "content"
    });
  } );
END

	if($b_query_total or $t_query_total){
		print LINEAGE_HTML "window.onload = function () {\n";
	}
	if($b_query_total){
		print LINEAGE_HTML <<END;
window.chartB = new CanvasJS.Chart("chartContainerB",
{
	title: {
  	text: "% B-cell epitopes with given number of AA residue mismatches (of IEDB $b_query_total total)"
  },
	data: [
	{
		type: "column",
		dataPoints: [
END

		# Print the B-epitope data
		for(my $i = 1; $i <= $max_histogram_x; $i++){
			print LINEAGE_HTML "{ x: $i, y: ", (defined $b_mismatch_count{$i} ? $b_mismatch_count{$i} : 0)," }", ($i == $max_histogram_x ? "" : ",");
		}

		print LINEAGE_HTML <<END;
		]
	}
	]
});
chartB.render();
chartB.set("dataPointWidth",Math.ceil(chartB.axisX[0].bounds.width/chartB.data[0].dataPoints.length),true);
END
	} #end if($b_query_total)

	if($t_query_total){
		print LINEAGE_HTML <<END;
window.chartT = new CanvasJS.Chart("chartContainerT",
{
	title: {
  	text: "% T-cell epitopes with given number of AA residue mismatches (of IEDB $t_query_total total)"
  },
	data: [
	{
		type: "column",
		dataPoints: [
END

		# Print the T-epitope data
		for(my $i = 1; $i <= $max_histogram_x; $i++){
			print LINEAGE_HTML "{ x: $i, y: ", (defined $t_mismatch_count{$i} ? $t_mismatch_count{$i} : 0), "}", ($i == $max_histogram_x ? "" : ",");
		}

                print LINEAGE_HTML <<END;
                ]
        }
        ]
});
chartT.render();
chartT.set("dataPointWidth",Math.ceil(chartB.axisX[0].bounds.width/chartT.data[0].dataPoints.length),true);
END

	} # end if($t_query_total)


	if($b_query_total or $t_query_total){
		print LINEAGE_HTML <<END;
\$( window ).resize(function() {
        chartB.set("dataPointWidth",Math.ceil(chartB.axisX[0].bounds.width/chartB.data[0].dataPoints.length),true);
        chartT.set("dataPointWidth",Math.ceil(chartT.axisX[0].bounds.width/chartT.data[0].dataPoints.length),true);
});
}
END
	}

	print LINEAGE_HTML <<END;
  </script>
</head>
<body>
<h1>Functional mutation literature for <a href="$front_rel_url">SARS-CoV-2</a> lineage <a href="https://cov-lineages.org/lineages/lineage_$lineage_short_name.html">$lineage_name</a></h1>
END

	# If we have a page for the parent lineage (same thing with less dot suffices), provide the link. Same for sub lineages
	my @child_lineage_links;
	while($lineage_short_name =~ s/\.\d+$//){
		if(exists $lineage_mutations{$lineage_short_name}){
			print LINEAGE_HTML "<div class='crossref'>Jump to ancestral lineage report: <a href='lineage-$lineage_short_name.html'>$lineage_short_name</a></div>\n";
		}
	}
	for my $sublineage_name (sort keys %lineage_mutations){
		if($sublineage_name ne $lineage_name and $sublineage_name =~ /^\Q$lineage_name\E\b/){
			push @child_lineage_links, "<a href='lineage-$sublineage_name.html'>$sublineage_name</a>";
		}
	}
	if(@child_lineage_links){
		print LINEAGE_HTML "<div class='crossref'>Jump to sublineage report: ", join(" &mdash; ", @child_lineage_links), "</div>\n";
	}

	if($b_query_total or $t_query_total){
		print LINEAGE_HTML "<div>\n<table>\n<tr>";
		if($b_query_total){
			print LINEAGE_HTML "<td><div id=\"chartContainerB\" class=\"epitope_mismatch_histogram\"></div></td>\n";
		}
		if($t_query_total){
			print LINEAGE_HTML "<td><div id=\"chartContainerT\" class=\"epitope_mismatch_histogram\"></div></td>\n";
		}
		print LINEAGE_HTML <<END;
</tr>
<tr><td colspan="2" align="center">Calculated using representative genome <a href="$lineage_name.acknowledgements.pdf">$epi_id</a> downloaded from <a href="http://www.gisaid.org">GISAID</a> on $epi_date, 
aligned using <a href="http://github.com/nodrogluap/pokay">Pokay</a> to an <a href="http://www.iedb.org">IEDB</a> antigens file dated $iedb_date. <a href="$lineage_name.b.html">See full alignments and disclaimers 
for epitope alignments</a></td></tr>
</table>
</div>
END
	} # end if($b_query_total or $t_query_total)

	print LINEAGE_HTML "<div id='accordion'>\n";
	# Collect all the constellation evidence that fits the mutations in this lineage
	my %effect2gene2constellation2lineage_lit_refs;
	for my $constellation (keys %constellations2categories){
		# See of this constellation is in this lineage's spec for this gene 
		my ($gene, $mutations_spec) = split /:/, $constellation;

		next unless exists $lineage_mutations{$lineage_name}->{$gene};
		my $lineage_mutations_ref = $lineage_mutations{$lineage_name}->{$gene};
		my $lineage_optional_mutations_ref = $lineage_optional_mutations{$lineage_name}->{$gene};
		my $missing_mutations = 0;
		for my $mutation (split /;/, $mutations_spec){
			if(not grep {$mutation eq $_} (@$lineage_mutations_ref, @$lineage_optional_mutations_ref)){
				$missing_mutations++;
				last;
			}
		}
		next if $missing_mutations;

		for my $short_effect_desc (keys %{$constellations2categories{$constellation}}){
			my $lit_refs = $constellations2categories{$constellation}->{$short_effect_desc};
			$effect2gene2constellation2lineage_lit_refs{$short_effect_desc}->{$gene} //= {};
			$effect2gene2constellation2lineage_lit_refs{$short_effect_desc}->{$gene}->{$constellation} //= [];
			for my $lit_ref (@$lit_refs){	
				push @{$effect2gene2constellation2lineage_lit_refs{$short_effect_desc}->{$gene}->{$constellation}}, $lit_ref;
			}
		}
	}

	for my $short_effect_desc (sort keys %effect2gene2constellation2lineage_lit_refs){
		my $short_effect_desc_pretty = prettify_effect_text($short_effect_desc);
		my $num_references = 0;
		for my $gene (@gene_pretty){ # keep in genome order
			next unless exists $effect2gene2constellation2lineage_lit_refs{$short_effect_desc}->{$gene};
			for my $constellation (sort keys %{$effect2gene2constellation2lineage_lit_refs{$short_effect_desc}->{$gene}}){
				for my $lit_ref (@{$effect2gene2constellation2lineage_lit_refs{$short_effect_desc}->{$gene}->{$constellation}}){
					$num_references++;
				}
			}
		}	

   		print LINEAGE_HTML "<h3 id=\"$short_effect_desc\">$short_effect_desc_pretty - $num_references literature reference", ($num_references > 1 ? "s" : ""), "</h3><div>\n";
		for my $gene (@gene_pretty){
			next unless exists $effect2gene2constellation2lineage_lit_refs{$short_effect_desc}->{$gene};
			print LINEAGE_HTML "<h3>Gene $gene &mdash; $gene_pretty{$gene}</h3><ul class=\"lit_list\">\n";

    			for my $constellation (sort {&constellation_cmp($a,$b)} keys %{$effect2gene2constellation2lineage_lit_refs{$short_effect_desc}->{$gene}}){
				print LINEAGE_HTML hgvsify($constellation)."<br/><ul class=\"lit_list\">\n";
    				for my $lit_ref (@{$effect2gene2constellation2lineage_lit_refs{$short_effect_desc}->{$gene}->{$constellation}}){
      					print LINEAGE_HTML htmlify_lit_ref($lit_ref);
				}
				print LINEAGE_HTML "</ul>\n";
    			}
			print LINEAGE_HTML "</ul>\n";
		}

		print LINEAGE_HTML "</div>\n";
	}
	print LINEAGE_HTML "</div><div class=\"acknowledgement\">Automatically generated using text2ui script and data from <a href=\"https://github.com/nodrogluap/pokay\">Pokay</a></div></body>\n</html>\n";
	close(LINEAGE_HTML);
}

copy("$text_indir/../constellations.css", "$html_outdir/constellations.css")
  or die "Cannot copy image file to $html_outdir/constellations.css: $!\n";
copy("$text_indir/../genome.png", "$html_outdir/genome.png")
  or die "Cannot copy image file to $html_outdir/genome.png: $!\n";
copy("$text_indir/../imageMapResizer.min.js", "$html_outdir/imageMapResizer.min.js")
  or die "Cannot copy script file to $html_outdir/imageMapResizer.min.js: $!\n";

# Sort constellations by position
sub constellation_cmp{
  my ($a, $b) = @_;
  my ($pos_a) = $a =~ /^.*?(\d+)/;
  my ($pos_b) = $b =~ /^.*?(\d+)/;
  return (($pos_a <=> $pos_b) or ($a cmp $b)); # position-based or fall back to a lexical if same position 
}

sub hgvsify{ 
  my $simplified_constellation_key = $_[0];
  while($simplified_constellation_key =~ /[A-Z](\d+)del;[A-Z](\d+)del/g){
    next if $1 != $2-1; # merge only neighbouring residue deletions
    my $last_aa = $2;
    $simplified_constellation_key =~ s/(\d+)del;([A-Z])$last_aa del/"$1_$2".$last_aa."del"/xe;
  }
  $simplified_constellation_key =~ s/(\d+)_(?:\d+_)+(\d+)del/$1."_".$2."del"/eg;
  $simplified_constellation_key =~ s/:/:p./;
  if($simplified_constellation_key =~ /,/){
    $simplified_constellation_key =~ tr/,/;/;
    $simplified_constellation_key =~ s/p\.(.*)/p.[$1]/;
  }
  return $simplified_constellation_key;
}

sub prettify_effect_text{
    	my $short_effect_desc_pretty = $_[0];
    	$short_effect_desc_pretty =~ s/^[a-z]/uc($&)/e; # Upper case first word in title
    	$short_effect_desc_pretty =~ tr/_/ /; # underscores to spaces
	return $short_effect_desc_pretty;
}

sub htmlify_lit_ref{
	if($_[0]->[2] eq "unpublished"){
		warn "Got an unpublished reference: ", @{$_[0]}, "\n";
	}
	return "<li class=\"lit_item\"><span class=\"lit_summary\">".
                                $_[0]->[0]. " (". # text desc
				($_[0]->[2] eq "unpublished" ? $_[0]->[1]. ", unpublished" : 
                                "<a href=\"".$_[0]->[2]."\" class=\"lit_link\">".$_[0]->[1]."</a>").
				")</span></li>\n";
}

sub textify_lit_ref {
    my ($ref) = @_;

    if ($ref->[2] eq "unpublished") {
        warn "Got an unpublished reference: ", @$ref, "\n";
    }

    my $text_summary = $ref->[0] . " (" .
                       ($ref->[2] eq "unpublished" ? $ref->[1] . ", unpublished" :
                       $ref->[1] . " - link available at: " . $ref->[2]) .
                       ")";
    
    # Clean the summary for JSON
    $text_summary =~ s/\"/\\\"/g;  # Escape double quotes
    $text_summary =~ s/<br\s*\/?>/ /g;  # Replace <br/> with newline
    $text_summary =~ s/<[^>]+>//g;  # Remove any other HTML tags

    return $text_summary;
}