")
+
+ exon_lines.append(f"Exon region: {exon_region} ")
+ exon_lines.append(f"Nucleotide percent identity: {pid} | BLOSUM: {blosum} ")
+ exon_lines.append(f"Intersects assembly gaps: {gaps} ")
+ exon_lines.append(f"Exon alignment class: {ali_class} ")
+ exon_lines.append(f"Detected within expected region ({exp_reg}): {in_exp_region} ")
+ exon_lines.append(f" ")
+ exon_lines.append(f"Sequence alignment between reference and query exon: ")
+ exon_lines.append(f"{ali} ")
+ exon_line = "".join(exon_lines)
+ transcript_lines.append(exon_line)
+ transcript_line = "".join(transcript_lines)
+ trans_to_exon_line[transcript] = transcript_line
+ return trans_to_exon_line
+
+
+def __get_ref_transcript(q_trans):
+ return ".".join(q_trans.split(".")[:-1])
+
+
+def get_ref_ts_to_link(p_dir):
+ """Get ref_transcript: link dict."""
+ ret = {}
+ ref_2bit_link = os.path.join(p_dir, "t2bit.link")
+ ref_2bit_path = os.readlink(ref_2bit_link)
+ ref_2bit_basename = os.path.basename(ref_2bit_path)
+ ref_genome_gname = ".".join(ref_2bit_basename.split(".")[:-1])
+ ref_toga_dir = f"/projects/hillerlab/genome/gbdb-HL/{ref_genome_gname}/TOGA/"
+ ref_ts_to_link_path = os.path.join(ref_toga_dir, "toga.transcript2Link.txt")
+ # dirname = os.path.dirname(ref_2bit_path)
+ # chrom_sizes_path = os.path.join(dirname, "chrom.sizes")
+ if not os.path.isfile(ref_ts_to_link_path):
+ return ret, ref_genome_gname
+ f = open(ref_ts_to_link_path, "r")
+ for line in f:
+ ld = line.rstrip().split("\t")
+ ret[ld[0]] = ld[1]
+ f.close()
+ return ret, ref_genome_gname
+
+
+def get_que_chrom_sizes(p_dir):
+ que_2bit_link = os.path.join(p_dir, "q2bit.link")
+ que_2bit_path = os.readlink(que_2bit_link)
+ dirname = os.path.dirname(que_2bit_path)
+ chrom_sizes_path = os.path.join(dirname, "chrom.sizes")
+ que_2bit_basename = os.path.basename(que_2bit_path)
+ que_genome_gname = ".".join(que_2bit_basename.split(".")[:-1])
+ return chrom_sizes_path, que_genome_gname
+
+
+def merge_all_tables(project_dir, r_ts_to_link):
+ """Merge all tabs in ready-to-merge bed file."""
+ print("Joining tsv for bigbed")
+ tabs_dir = os.path.join(project_dir, "tabs")
+ sorted_bed = os.path.join(tabs_dir, "query_annot_sorted.bed")
+ toga_data_file = os.path.join(tabs_dir, "togaData.tab")
+ toga_inact_file = os.path.join(tabs_dir, "togaInactMut.tab")
+ toga_nucl_file = os.path.join(tabs_dir, "togaNucl.tab")
+ transcripts_ordered, t_to_bed = __get_ts_from_bed(sorted_bed)
+ ts_to_data_track = __get_ts_to_toga_data(toga_data_file)
+ ts_to_inact_data = __get_ts_to_toga_inact(toga_inact_file)
+ ts_to_exon_data = __get_ts_to_toga_exons(toga_nucl_file)
+
+
+ bigbed_material = os.path.join(tabs_dir, "query_annot_for_big.tsv")
+ f = open(bigbed_material, "w")
+ for transcript in transcripts_ordered:
+ ref_trans = __get_ref_transcript(transcript)
+ bed_fields = t_to_bed[transcript]
+ data_track = ts_to_data_track[transcript]
+ inact_line = ts_to_inact_data.get(transcript, " ")
+ exon_line = ts_to_exon_data[transcript]
+ link = r_ts_to_link.get(ref_trans, ref_trans) # just reference transcript if no link
+ combined_datum_lst = bed_fields + data_track + [link] + [inact_line] + [exon_line]
+ combined_datum = "\t".join(combined_datum_lst)
+ f.write(f"{combined_datum}\n")
+ f.close()
+
+
+def sort_and_make_bb(project_dir, chrom_sizes):
+ """bed12+22 as the result
+ sort -k1,1 -k2,2n query_annot_for_big.tsv > query_annot_for_big_sorta.tsv
+ bedToBigBed -type=bed12+22 query_annot_for_big_sorta.tsv
+ /projects/hillerlab/genome/gbdb-HL/mm10/chrom.sizes query_annot.bb -tab
+ """
+ print("Making bigbed...")
+ tabs_dir = os.path.join(project_dir, "tabs")
+ bigbed_material_not_sorted = os.path.join(tabs_dir, "query_annot_for_big.tsv")
+ bigbed_material = os.path.join(tabs_dir, "query_annot_for_big__sorted.tsv")
+ sort_cmd = f"sort -k1,1 -k2,2n {bigbed_material_not_sorted} > {bigbed_material}"
+ rc = subprocess.call(sort_cmd, shell=True)
+ if rc != 0:
+ sys.exit(f"Error! Command:\n{sort_cmd}\ncrashed")
+ output = os.path.join(tabs_dir, "query_annotation.bb")
+ # now transform into bigbed
+ print("Calling:")
+ bb_cmd = f"bedToBigBed -type=bed12+22 {bigbed_material} {chrom_sizes} {output} -tab -extraIndex=name -as={SCHEMA_LOCATION}"
+ print(bb_cmd)
+ rc = subprocess.call(bb_cmd, shell=True)
+ if rc != 0:
+ sys.exit(f"Error! Command:\n{bb_cmd}\ncrashed")
+
+ # make ix.txt
+ query_annot = os.path.join(project_dir, "query_annotation.bed")
+ make_ix_script = os.path.join(LOCATION, "get_names_from_bed.py")
+ ix_txt = os.path.join(tabs_dir, "query_annotation.ix.txt")
+ ix_cmd = f"{make_ix_script} {query_annot} | sort -u > {ix_txt}"
+ rc = subprocess.call(ix_cmd, shell=True)
+ if rc != 0:
+ sys.exit(f"Error! Command:\n{ix_cmd}\ncrashed")
+
+ # make .ix and .ixx indexes
+ ix_path = os.path.join(tabs_dir, "query_annotation.bb.ix")
+ ixx_path = os.path.join(tabs_dir, "query_annotation.bb.ixx")
+ ixixx_cmd = f"ixIxx {ix_txt} {ix_path} {ixx_path}"
+ rc = subprocess.call(ixixx_cmd, shell=True)
+ if rc != 0:
+ sys.exit(f"Error! Command:\n{ixixx_cmd}\ncrashed")
+
+
+def cleanup(project_dir, dont):
+ """Clean project dir up."""
+ if dont: # do not cleanup
+ return
+ tabs_dir = os.path.join(project_dir, "tabs")
+ bigbed_material_not_sorted = os.path.join(tabs_dir, "query_annot_for_big.tsv")
+ bigbed_material = os.path.join(tabs_dir, "query_annot_for_big__sorted.tsv")
+ toga_info_file = os.path.join(tabs_dir, "togaInfo.tab")
+ toga_inact_feat_file = os.path.join(tabs_dir, "togaInactFeat.tab")
+ toga_plot_file = os.path.join(tabs_dir, "togaPlot.tab")
+ toga_prot_file = os.path.join(tabs_dir, "togaProt.tab")
+ sorted_bed = os.path.join(tabs_dir, "query_annot_sorted.bed")
+ toga_data_file = os.path.join(tabs_dir, "togaData.tab")
+ toga_inact_file = os.path.join(tabs_dir, "togaInactMut.tab")
+ toga_nucl_file = os.path.join(tabs_dir, "togaNucl.tab")
+ to_rm = [bigbed_material_not_sorted, bigbed_material, toga_info_file,
+ toga_inact_feat_file, toga_plot_file, toga_prot_file,
+ sorted_bed, toga_data_file, toga_inact_file, toga_nucl_file]
+ for path in to_rm:
+ os.remove(path) if os.path.isfile(path) else None
+
+
+def __ssh_mkdir(uname, dirname):
+ cmd = f"ssh {uname}@genome mkdir -p {dirname}"
+ print(f"Calling {cmd}")
+ subprocess.call(cmd, shell=True)
+
+
+def __ssh_ln(uname, src, dest):
+ """Call ln on delta"""
+ cmd = f"ssh {uname}@genome ln -sf {src} {dest}"
+ print(f"Calling {cmd}")
+ subprocess.call(cmd, shell=True)
+
+
+def __get_ccase_refname(refname):
+ if refname == "hg38":
+ return "Hg38"
+ elif refname == "mm10":
+ return "Mm10"
+ elif refname == "mm39":
+ return "Mm39"
+ else:
+ return refname
+
+
+def load_bb_track(project_dir, dont, ref, que, bb_ver):
+ if dont: # do not load anything
+ return
+ tabs_dir = os.path.join(project_dir, "tabs")
+ bb_file = os.path.join(tabs_dir, "query_annotation.bb")
+ uname = getpass.getuser()
+ # scp query_annotation.bb genome:/genome/gbdb-HL/$ref/TOGA/vs_$query/HLTOGAannotVs$refv1.bb
+ ref_dir__genome = f"/genome/gbdb-HL/{ref}"
+ toga_ref_dir__genome = f"{ref_dir__genome}/TOGA"
+ # make if needed
+ __ssh_mkdir(uname, toga_ref_dir__genome)
+ vs_que_dir__genome = f"{toga_ref_dir__genome}/vs_{que}"
+ __ssh_mkdir(uname, vs_que_dir__genome)
+ # transfer data
+ bb_file = os.path.join(tabs_dir, "query_annotation.bb")
+ ix_file = os.path.join(tabs_dir, "query_annotation.bb.ix")
+ ixx_file = os.path.join(tabs_dir, "query_annotation.bb.ixx")
+
+ ref_camel_case = __get_ccase_refname(ref)
+ bb_filename__genome = f"HLTOGAannotVs{ref_camel_case}{bb_ver}.bb"
+ ix_filename__genome = f"HLTOGAannotVs{ref_camel_case}{bb_ver}.ix"
+ ixx_filename__genome = f"HLTOGAannotVs{ref_camel_case}{bb_ver}.ixx"
+
+ bb_path_genome = f"{vs_que_dir__genome}/{bb_filename__genome}"
+ ix_path_genome = f"{vs_que_dir__genome}/{ix_filename__genome}"
+ ixx_path_genome = f"{vs_que_dir__genome}/{ixx_filename__genome}"
+
+ rsync_cmd = f"rsync -av {bb_file} {uname}@genome:{bb_path_genome}"
+ subprocess.call(rsync_cmd, shell=True)
+
+ rsync_cmd = f"rsync -av {ix_file} {uname}@genome:{ix_path_genome}"
+ subprocess.call(rsync_cmd, shell=True)
+
+ rsync_cmd = f"rsync -av {ixx_file} {uname}@genome:{ixx_path_genome}"
+ subprocess.call(rsync_cmd, shell=True)
+
+ # make link in /var/www
+ www_data__genome = "/var/www/data"
+ www_que__genome = f"{www_data__genome}/{que}"
+ __ssh_mkdir(uname, www_que__genome)
+
+ www_bb_link_dest = f"{www_que__genome}/{bb_filename__genome}"
+ www_ix_link_dest = f"{www_que__genome}/{ix_filename__genome}"
+ www_ixx_link_dest = f"{www_que__genome}/{ixx_filename__genome}"
+
+ __ssh_ln(uname, bb_path_genome, www_bb_link_dest)
+ __ssh_ln(uname, ix_path_genome, www_ix_link_dest)
+ __ssh_ln(uname, ixx_path_genome, www_ixx_link_dest)
+
+
+def main():
+ """Entry point."""
+ args = parse_args()
+ r_ts_to_link, ref_name = get_ref_ts_to_link(args.project_dir)
+ chrom_sizes, que_name = get_que_chrom_sizes(args.project_dir)
+ # call generate_tab_files.py and make_togaPlot.py
+ call_gen_tab_files(args.project_dir)
+ call_gen_plot_files(args.project_dir, no_plots=args.no_plots)
+ # merge 4 transcriptID-related tables into a single one
+ make_toga_data(args.project_dir)
+
+ # OK, now making bigBed
+ # togaNucl.tab and togaInactMut.tab
+ # -> merge into HTML_formatted strings
+ # -> add to togaData
+ # save to as.tab for bedToBigBed
+ make_sorted_bed(args.project_dir)
+
+ merge_all_tables(args.project_dir, r_ts_to_link)
+ sort_and_make_bb(args.project_dir, chrom_sizes)
+
+ cleanup(args.project_dir, args.do_not_cleanup)
+
+
+if __name__ == "__main__":
+ main()
diff --git a/ucsc_browser_visualisation/make_sql_data.py b/ucsc_browser_visualisation/make_sql_data.py
deleted file mode 100755
index 41eba17..0000000
--- a/ucsc_browser_visualisation/make_sql_data.py
+++ /dev/null
@@ -1,131 +0,0 @@
-#!/usr/bin/env python3
-"""Perform all operations to create SQL tables for UCSC browser."""
-import os
-import sys
-import argparse
-import subprocess
-
-GENERATE_TAB_FILES = "generate_tab_files.py"
-GENERATE_PLOTS = "make_togaPlot.py"
-LOCATION = os.path.dirname(__file__)
-
-IFEAT_PLACEHOLDER = ["0.0" for _ in range(6)]
-PROT_PLACEHOLDER = ["NO_DATA"]
-UNDEF = "UNDEFINED"
-SVG_PLACEHOLDER = """
-"""
-PLOT_PLACEHOLDER = [SVG_PLACEHOLDER, ]
-
-
-def parse_args():
- """Parse args."""
- app = argparse.ArgumentParser()
- app.add_argument("project_dir", help="Directory containing TOGA output")
- if len(sys.argv) < 2:
- app.print_help()
- sys.exit(0)
- args = app.parse_args()
- return args
-
-
-def call_gen_tab_files(project_dir):
- """Call Generate tab files subprocess."""
- print(f"Calling {GENERATE_TAB_FILES}")
- out_dir = os.path.join(project_dir, "tabs")
- os.mkdir(out_dir) if not os.path.isdir(out_dir) else None
- exe_ = os.path.join(LOCATION, GENERATE_TAB_FILES)
- cmd = f"{exe_} {project_dir} {out_dir}"
- rc = subprocess.call(cmd, shell=True)
- if rc != 0:
- raise ValueError(f"Command {cmd} died - cannot generate tab files")
- print("5 of 6 tab files generated")
-
-
-def call_gen_plot_files(project_dir):
- """Generate svg plots."""
- print(f"Calling {GENERATE_PLOTS}")
- out_file = os.path.join(project_dir, "tabs", "togaPlot.tab")
- exe_ = os.path.join(LOCATION, GENERATE_PLOTS)
- cmd = f"{exe_} {project_dir} {out_file}"
- rc = subprocess.call(cmd, shell=True)
- if rc != 0:
- raise ValueError(f"Command {cmd} died - cannot generate plots")
- print("All tab files generated")
-
-
-def read_tab_file(tab_file):
- """Make trans_id: data track."""
- trans_to_dat = {}
- f = open(tab_file, "r")
- for line in f:
- line_data = line.rstrip().split("\t")
- trans = line_data[0]
- data = line_data[1:]
- trans_to_dat[trans] = data
- f.close()
- return trans_to_dat
-
-
-
-def __gen_info_placeholder(p_name):
- """Generate info placeholder."""
- t_name = ".".join(p_name.split(".")[:-1])
- placeholder = [t_name, UNDEF, UNDEF, "0.0", "0", "0.0", "0.0", "0.0", "0.0", "0.0", UNDEF]
- return placeholder
-
-
-def make_toga_data(project_dir):
- """Merge togaPlot, togaProt, togaInfo and togaInactFeat tables into togaData."""
- tabs_dir = os.path.join(project_dir, "tabs")
- toga_info_file = os.path.join(tabs_dir, "togaInfo.tab")
- toga_inact_feat_file = os.path.join(tabs_dir, "togaInactFeat.tab")
- toga_plot_file = os.path.join(tabs_dir, "togaPlot.tab")
- toga_prot_file = os.path.join(tabs_dir, "togaProt.tab")
- out_file = os.path.join(tabs_dir, "togaData.tab")
-
- print("Reading tab files to merge")
- trans_to_info = read_tab_file(toga_info_file)
- trans_to_ifeat = read_tab_file(toga_inact_feat_file)
- trans_to_plot = read_tab_file(toga_plot_file)
- trans_to_prot = read_tab_file(toga_prot_file)
-
- transcripts = trans_to_info.keys()
- print(f"Got data for {len(transcripts)} transcripts")
-
- f = open(out_file, "w")
- for t in transcripts:
- info = trans_to_info.get(t)
- if info is None:
- info = __gen_info_placeholder(t)
- ifeat = trans_to_ifeat.get(t, IFEAT_PLACEHOLDER.copy())
- prot = trans_to_prot.get(t, PROT_PLACEHOLDER.copy())
- plot = trans_to_plot.get(t, PLOT_PLACEHOLDER.copy())
- data_lst = [t] + info + ifeat + prot + plot
- tab_str = "\t".join(data_lst)
- f.write(tab_str)
- f.write("\n")
- f.close()
- # do cleanup
- os.remove(toga_info_file)
- os.remove(toga_inact_feat_file)
- os.remove(toga_plot_file)
- os.remove(toga_prot_file)
- print("Done")
-
-
-def main():
- """Entry point."""
- args = parse_args()
- # call generate_tab_files.py and make_togaPlot.py
- call_gen_tab_files(args.project_dir)
- call_gen_plot_files(args.project_dir)
- # merge 4 transcriptID-related tables into a single one
- make_toga_data(args.project_dir)
-
-
-if __name__ == "__main__":
- main()
diff --git a/ucsc_browser_visualisation/readme.md b/ucsc_browser_visualisation/readme.md
index 4241c7e..c2f53d2 100644
--- a/ucsc_browser_visualisation/readme.md
+++ b/ucsc_browser_visualisation/readme.md
@@ -76,13 +76,22 @@ We will add another condition to handle TOGA tracks.
Put another condition as shown here:
+
```c
-else if (startsWith("HLTOGA", table) && hTableExists(database, "TOGAData"))
-{
- doHillerLabTOGAGene(tdb, item);
-}
+else if (startsWith("HLTOGAannot", trackHubSkipHubName(table)))
+ {
+ doHillerLabTOGAGene(database, tdb, item, table);
+ }
```
+
+[//]: # (```c)
+[//]: # (else if (startsWith("HLTOGA", table) && hTableExists(database, "TOGAData")))
+[//]: # ({)
+[//]: # ( doHillerLabTOGAGene(database, tdb, item, table);)
+[//]: # (}
+[//]: # (```)
+
Then re-build your browser.
### Loading TOGA tables
@@ -92,26 +101,18 @@ Then re-build your browser.
${project_dir}: directory containing TOGA results and intermediate data.
```shell
-./ucsc_browser_visualisation/make_sql_data.py ${project_dir}```
+./ucsc_browser_visualisation/make_bigbed_data_public.py ${project_dir}```
# Wait, it will take a few minutes, most likely less than an hour
```
#### Load tab files to browser database
-${tab_files_dir} = ${project_dir}/tabs
-${query} - annotated genome identifier.
-Use *.sql files located in the ucsc_browser_visualisation directory.
-Call the following commands:
+Transfer created *.bb, *.ix, and *.ixx files to the machine hosting your instance of UCSC genome browser, if need be.
+Create a table schema for bigBed tracks, using bigDataUrl field to specify the bigBed URL, and searchTrix field to specify the *.ix file URL (no need to specify *.ixx file separately, it just should be located in the same directory with the *.ix file.)
-```shell
-hgLoadSqlTab ${query} TOGAData togaData.sql ${tab_files_dir}/togaInfo.tab
-hgLoadSqlTab ${query} TOGANucl togaNucl.sql ${tab_files_dir}/togaNucl.tab
-hgLoadSqlTab ${query} TOGAInactMut togaInactMut.sql ${tab_files_dir}/togaInactMut.tab
+Please also specify the following fields:
```
-
-Also load bed file using hgLoadBed command.
-Create a track starting with HLTOGA, for instance HLTOGAannotation.
-
-```shell script
-hgLoadBed ${query} HLTOGAannotation ${project_dir}/query_annotation.bed
+type bigBed 12 +
+labelFields name
+searchIndex name
```
diff --git a/ucsc_browser_visualisation/togaClick.c b/ucsc_browser_visualisation/togaClick.c
index e2f10e6..e3d53c8 100644
--- a/ucsc_browser_visualisation/togaClick.c
+++ b/ucsc_browser_visualisation/togaClick.c
@@ -3,6 +3,44 @@
#include "hgc.h"
#include "togaClick.h"
#include "string.h"
+#include "htmshell.h"
+#include "chromAlias.h"
+
+
+struct togaDataBB *togaDataBBLoad(char **row)
+/* Load a togaData from row fetched with select * from togaData
+ * from database. Dispose of this with togaDataFree(). */
+{
+ struct togaDataBB *ret;
+ AllocVar(ret);
+ ret->projection = cloneString(row[0]);
+ ret->ref_trans_id = cloneString(row[1]);
+ ret->ref_region = cloneString(row[2]);
+ ret->query_region = cloneString(row[3]);
+ ret->chain_score = cloneString(row[4]);
+
+ ret->chain_synteny = cloneString(row[5]);
+ ret->chain_flank = cloneString(row[6]);
+ ret->chain_gl_cds_fract = cloneString(row[7]);
+ ret->chain_loc_cds_fract = cloneString(row[8]);
+ ret->chain_exon_cov = cloneString(row[9]);
+
+ ret->chain_intron_cov = cloneString(row[10]);
+ ret->status = cloneString(row[11]);
+ ret->perc_intact_ign_M = cloneString(row[12]);
+ ret->perc_intact_int_M = cloneString(row[13]);
+ ret->intact_codon_prop = cloneString(row[14]);
+
+ ret->ouf_prop = cloneString(row[15]);
+ ret->mid_intact = cloneString(row[16]);
+ ret->mid_pres = cloneString(row[17]);
+ ret->prot_alignment = cloneString(row[18]);
+ ret->svg_line = cloneString(row[19]);
+ ret->ref_link = cloneString(row[20]);
+ ret->inact_mut_html_table = cloneString(row[21]);
+ ret->exon_ali_html = cloneString(row[22]);
+ return ret;
+}
struct togaData *togaDataLoad(char **row)
@@ -38,6 +76,43 @@ struct togaData *togaDataLoad(char **row)
}
+void togaDataBBFree(struct togaDataBB **pEl)
+/* Free a single dynamically allocated togaDatasuch as created
+ * with togaDataLoad(). */
+{
+ struct togaDataBB *el;
+
+ if ((el = *pEl) == NULL) return;
+ freeMem(el->projection);
+ freeMem(el->ref_trans_id);
+ freeMem(el->ref_region);
+ freeMem(el->query_region);
+ freeMem(el->chain_score);
+
+ freeMem(el->chain_synteny);
+ freeMem(el->chain_flank);
+ freeMem(el->chain_gl_cds_fract);
+ freeMem(el->chain_loc_cds_fract);
+ freeMem(el->chain_exon_cov);
+
+ freeMem(el->chain_intron_cov);
+ freeMem(el->status);
+ freeMem(el->perc_intact_ign_M);
+ freeMem(el->perc_intact_int_M);
+ freeMem(el->intact_codon_prop);
+
+ freeMem(el->ouf_prop);
+ freeMem(el->mid_intact);
+ freeMem(el->mid_pres);
+ freeMem(el->prot_alignment);
+ freeMem(el->svg_line);
+ freeMem(el->ref_link);
+ freeMem(el->inact_mut_html_table);
+ freeMem(el->exon_ali_html);
+ freez(pEl);
+}
+
+
void togaDataFree(struct togaData **pEl)
/* Free a single dynamically allocated togaDatasuch as created
* with togaDataLoad(). */
@@ -168,10 +243,207 @@ Prefix must be HLTOGAannot */
}
-void doHillerLabTOGAGene(struct trackDb *tdb, char *item, char *table_name)
+void HLprintQueryProtSeqForAli(char *proteinAlignment) {
+ // take protein sequence alignment
+ // print only the query sequence
+ char *str = proteinAlignment;
+ int printed_char_num = 0;
+ while ((str = strstr(str, "que:")) != NULL)
+ {
+ str += 10;
+ char ch;
+ while ((ch = *str++) != '<') {
+ if (ch != '-') {
+ putchar(ch);
+ ++printed_char_num;
+ }
+ if (printed_char_num == 80) {
+ printed_char_num = 0;
+ printf(" ");
+ }
+ }
+ }
+}
+
+
+
+void doHillerLabTOGAGeneBig(char *database, struct trackDb *tdb, char *item, char *table_name)
+/* Put up TOGA Gene track info. */
+// To think about -> put into a single bigBed
+// string: HTML formatted inact mut
+// string: HTML formatted exon ali section
+{
+int start = cartInt(cart, "o");
+int end = cartInt(cart, "t");
+char *chrom = cartString(cart, "c");
+char *fileName = bbiNameFromSettingOrTable(tdb, NULL, tdb->table);
+struct bbiFile *bbi = bigBedFileOpenAlias(hReplaceGbdb(fileName), chromAliasFindAliases);
+struct lm *lm = lmInit(0);
+struct bigBedInterval *bbList = bigBedIntervalQuery(bbi, chrom, start, end, 0, lm);
+struct bigBedInterval *bb;
+char *fields[bbi->fieldCount];
+for (bb = bbList; bb != NULL; bb = bb->next)
+ {
+ if (!(bb->start == start && bb->end == end))
+ continue;
+
+ // our names are unique
+ char *name = cloneFirstWordByDelimiterNoSkip(bb->rest, '\t');
+ boolean match = (isEmpty(name) && isEmpty(item)) || sameOk(name, item);
+ if (!match)
+ continue;
+
+ char startBuf[16], endBuf[16];
+ bigBedIntervalToRow(bb, chrom, startBuf, endBuf, fields, bbi->fieldCount);
+ break;
+ }
+
+printf("
Projection %s
\n", item);
+struct togaDataBB *info = togaDataBBLoad(&fields[11]); // Bogdan: why 11? 0-11 are bed-like fields likely
+
+printf("Reference transcript: %s ", info->ref_link);
+printf("Genomic locus in reference: %s \n", info->ref_region);
+printf("Genomic locus in query: %s \n", info->query_region);
+
+printf("Projection classification: %s \n", info->status);
+printf("Probability that query locus is orthologous: %s \n", info->chain_score);
+// list of chain features (for orthology classification)
+printf("Show features used for ortholog probability\n");
+printf("
\n");
+printf("
\n");
+printf("
Synteny (log10 value): %s
\n", info->chain_synteny);
+printf("
Global CDS fraction: %s
\n", info->chain_gl_cds_fract);
+printf("
Local CDS fraction: %s
\n", info->chain_loc_cds_fract);
+printf("
Local intron fraction: %s
\n", info->chain_intron_cov);
+printf("
Local CDS coverage: %s
\n", info->chain_exon_cov);
+printf("
Flank fraction: %s
\n", info->chain_flank);
+printf("
\n");
+
+printf(" \nFeature description:\n");
+printf("For each projection (one reference transcript and one overlapping chain),\n");
+printf("TOGA computes the following features by intersecting the reference coordinates of aligning\n");
+printf("blocks in the chain with different gene parts (coding exons, UTR (untranslated region) exons, introns)\n");
+printf("and the respective intergenic regions.\n \n");
+
+printf("We define the following variables:\n
\n");
+printf("
c: number of reference bases in the intersection between chain blocks and coding exons of the gene under consideration.
\n");
+printf("
C: number of reference bases in the intersection between chain blocks and coding exons of all genes.
\n");
+printf("
a: number of reference bases in the intersection between chain blocks and coding exons and introns of the gene under consideration.
\n");
+printf("
A: number of reference bases in the intersection between chain blocks and coding exons and introns of all genes and the intersection\n");
+printf("between chain blocks and intergenic regions (excludes UTRs).
\n");
+printf("
f: number of reference bases in chain blocks overlapping the 10 kb flanks of the gene under consideration.\n");
+printf("Alignment blocks overlapping exons of another gene that is located in these 10 kb flanks are ignored.
\n");
+printf("
i: number of reference bases in the intersection between chain blocks and introns of the gene under consideration.
\n");
+printf("
CDS (coding sequence): length of the coding region of the gene under consideration.
\n");
+printf("
I: sum of all intron lengths of the gene under consideration.
\n");
+printf("
\n");
+printf("Using these variables, TOGA computes the following features:\n");
+printf("
\n");
+printf("
“global CDS fraction” as C / A. Chains with a high value have alignments that largely overlap coding exons,");
+printf("which is a hallmark of paralogous or processed pseudogene chains. In contrast, chains with a low value also align many ");
+printf("intronic and intergenic regions, which is a hallmark of orthologous chains.
\n");
+printf("
“local CDS fraction” as c / a. Orthologous chains tend to have a lower value, as intronic ");
+printf("regions partially align. This feature is not computed for single-exon genes.
\n");
+printf("
“local intron fraction” as i / I. Orthologous chains tend to have a higher value.");
+printf("This feature is not computed for single-exon genes.
\n");
+printf("
“flank fraction” as f / 20,000. Orthologous chains tend to have higher values,");
+printf("as flanking intergenic regions partially align. This feature is important to detect orthologous loci of single-exon genes.
\n");
+printf("
“synteny” as log10 of the number of genes, whose coding exons overlap by at least one base aligning");
+printf("blocks of this chain. Orthologous chains tend to cover several genes located in a conserved order, resulting in higher synteny values.
\n");
+printf("
“local CDS coverage” as c / CDS, which is only used for single-exon genes.
\n");
+printf("
\n");
+
+
+printf("\n
\n \n");
+htmlHorizontalLine();
+
+// show inact mut plot
+printf("
Visualization of inactivating mutations on exon-intron structure
\n");
+printf("%s \n", info->svg_line);
+printf(" Exons shown in grey are missing (often overlap assembly gaps).\nExons shown in");
+printf(" red or blue are deleted or do not align at all.\nRed indicates that the exon deletion ");
+printf("shifts the reading frame, while blue indicates that exon deletion(s) are framepreserving. \n");
+
+// GLP features
+printf("Show features used for transcript classification\n");
+printf("
\n");
+printf("
\n");
+printf("
Percent intact, ignoring missing sequence: %s
\n", info->perc_intact_ign_M);
+printf("
Percent intact, treating missing as intact sequence: %s
\n", info->perc_intact_int_M);
+printf("
Proportion of intact codons: %s
\n", info->intact_codon_prop);
+printf("
Percent of CDS not covered by this chain (0 unless the chain covers only a part of the gene): %s
\n");
+// printf("{protein seq of the query without dashes or other things. Should end with *}\n");
+printf("");
+HLprintQueryProtSeqForAli(info->prot_alignment);
+printf("\n \n\n
\n");
+
+// and show protein sequence
+htmlHorizontalLine();
+printf("
\n");
+// printf("%s\n", info->exon_ali_string);
+printf("%s\n", info->exon_ali_html);
+
+htmlHorizontalLine();
+
+// TODO: check whether I need this
+printf("%s", hgTracksPathAndSettings());
+hPrintf("");
+hPrintf("");
+hPrintf("");
+
+
+printTrackHtml(tdb); // and do I need this?
+}
+
+
+void doHillerLabTOGAGene(char *database, struct trackDb *tdb, char *item, char *table_name)
/* Put up TOGA Gene track info. */
{
- int start = cartInt(cart, "o");
+ //int start = cartInt(cart, "o");
char headerTitle[512];
char suffix[512];
strcpy(suffix, table_name);
@@ -180,6 +452,13 @@ void doHillerLabTOGAGene(struct trackDb *tdb, char *item, char *table_name)
genericHeader(tdb, headerTitle);
printf("
TOGA gene annotation
\n");
// htmlHorizontalLine();
+
+ if (startsWith("bigBed", tdb->type))
+ {
+ doHillerLabTOGAGeneBig(database, tdb, item, table_name);
+ return;
+ }
+
struct sqlConnection *conn = hAllocConn(database);
// define TOGA table names: initate with pre-defined prefixes
@@ -209,57 +488,103 @@ void doHillerLabTOGAGene(struct trackDb *tdb, char *item, char *table_name)
if ((row = sqlNextRow(sr)) != NULL) {
info = togaDataLoad(row); // parse sql output
// fill HTML template:
- printf("Projected via: %s ",
+ printf("Reference transcript: %s ",
info->ref_trans_id, info->ref_trans_id);
- printf("Region in reference: %s \n", info->ref_region);
- printf("Region in query: %s \n", info->query_region);
+ printf("Genomic locus in reference: %s \n", info->ref_region);
+ printf("Genomic locus in query: %s \n", info->query_region);
- printf("Projection class: %s \n", info->status);
- printf("Chain score: %s \n", info->chain_score);
+ printf("Projection classification: %s \n", info->status);
+ printf("Probability that query locus is orthologous: %s \n", info->chain_score);
// list of chain features (for orthology classification)
- printf("Show chain features for classification\n");
+ printf("Show features used for ortholog probability\n");
printf("
\n");
printf("
\n");
- printf("
Synteny: %s
\n", info->chain_synteny);
+ printf("
Synteny (log10 value): %s
\n", info->chain_synteny);
printf("
Global CDS fraction: %s
\n", info->chain_gl_cds_fract);
printf("
Local CDS fraction: %s
\n", info->chain_loc_cds_fract);
printf("
Local intron fraction: %s
\n", info->chain_intron_cov);
printf("
Local CDS coverage: %s
\n", info->chain_exon_cov);
printf("
Flank fraction: %s
\n", info->chain_flank);
+ printf("
\n");
+
+ printf(" \nFeature description:\n");
+ printf("For each projection (one reference transcript and one overlapping chain),\n");
+ printf("TOGA computes the following features by intersecting the reference coordinates of aligning\n");
+ printf("blocks in the chain with different gene parts (coding exons, UTR (untranslated region) exons, introns)\n");
+ printf("and the respective intergenic regions.\n \n");
+
+ printf("We define the following variables:\n
\n");
+ printf("
c: number of reference bases in the intersection between chain blocks and coding exons of the gene under consideration.
\n");
+ printf("
C: number of reference bases in the intersection between chain blocks and coding exons of all genes.
\n");
+ printf("
a: number of reference bases in the intersection between chain blocks and coding exons and introns of the gene under consideration.
\n");
+ printf("
A: number of reference bases in the intersection between chain blocks and coding exons and introns of all genes and the intersection\n");
+ printf("between chain blocks and intergenic regions (excludes UTRs).
\n");
+ printf("
f: number of reference bases in chain blocks overlapping the 10 kb flanks of the gene under consideration.\n");
+ printf("Alignment blocks overlapping exons of another gene that is located in these 10 kb flanks are ignored.
\n");
+ printf("
i: number of reference bases in the intersection between chain blocks and introns of the gene under consideration.
\n");
+ printf("
CDS (coding sequence): length of the coding region of the gene under consideration.
\n");
+ printf("
I: sum of all intron lengths of the gene under consideration.
\n");
+ printf("
\n");
+ printf("Using these variables, TOGA computes the following features:\n");
+ printf("
\n");
+ printf("
“global CDS fraction” as C / A. Chains with a high value have alignments that largely overlap coding exons,");
+ printf("which is a hallmark of paralogous or processed pseudogene chains. In contrast, chains with a low value also align many ");
+ printf("intronic and intergenic regions, which is a hallmark of orthologous chains.
\n");
+ printf("
“local CDS fraction” as c / a. Orthologous chains tend to have a lower value, as intronic ");
+ printf("regions partially align. This feature is not computed for single-exon genes.
\n");
+ printf("
“local intron fraction” as i / I. Orthologous chains tend to have a higher value.");
+ printf("This feature is not computed for single-exon genes.
\n");
+ printf("
“flank fraction” as f / 20,000. Orthologous chains tend to have higher values,");
+ printf("as flanking intergenic regions partially align. This feature is important to detect orthologous loci of single-exon genes.
\n");
+ printf("
“synteny” as log10 of the number of genes, whose coding exons overlap by at least one base aligning");
+ printf("blocks of this chain. Orthologous chains tend to cover several genes located in a conserved order, resulting in higher synteny values.
\n");
+ printf("
“local CDS coverage” as c / CDS, which is only used for single-exon genes.
\n");
+ printf("
\n");
+
+
printf("\n
\n \n");
htmlHorizontalLine();
// show inact mut plot
- printf("
Inactivating mutations plot
\n");
+ printf("
Visualization of inactivating mutations on exon-intron structure
\n");
printf("%s \n", info->svg_line);
+ printf(" Exons shown in grey are missing (often overlap assembly gaps).\nExons shown in");
+ printf(" red or blue are deleted or do not align at all.\nRed indicates that the exon deletion ");
+ printf("shifts the reading frame, while blue indicates that exon deletion(s) are framepreserving. \n");
// GLP features
- printf("Show GLP features\n");
+ printf("Show features used for transcript classification\n");
printf("
\n");
printf("
\n");
- printf("
Percent intact ignoring missing seq: %s
\n", info->perc_intact_ign_M);
- printf("
Percent intact (miss == intact): %s
\n", info->perc_intact_int_M);
- printf("
Intact codon proportion %s
\n", info->intact_codon_prop);
- printf("
Out of chain proportion: %s
\n", info->ouf_prop);
+ printf("
Percent intact, ignoring missing sequence: %s
\n", info->perc_intact_ign_M);
+ printf("
Percent intact, treating missing as intact sequence: %s
\n", info->perc_intact_int_M);
+ printf("
Proportion of intact codons: %s
\n", info->intact_codon_prop);
+ printf("
Percent of CDS not covered by this chain (0 unless the chain covers only a part of the gene): %s
\n", info->ouf_prop);
if (sameWord(info->mid_intact, ONE_))
{
- printf("
Middle 80 percent intact: %s
\n", YES_);
+ printf("
Middle 80 percent of CDS intact: %s
\n", YES_);
} else {
- printf("
Middle 80 percent intact: %s
\n", NO_);
+ printf("
Middle 80 percent of CDS intact: %s
\n", NO_);
}
if (sameWord(info->mid_pres, ONE_))
{
- printf("