CG - cancer genomics pipeline
This pipeline forms the basis of the computational work used in
“Cancer-associated mutations in the ribosomal protein L5 gene
dysregulate the HDM2/p53-mediated
ribosome biogenesis checkpoint”,
Oncogene (2020). https://doi.org/10.1038/s41388-020-1231-6
Please cite if you find it useful in your own work.
For the statistical interpretation of the clustering of
cancer-related mutations on the protein structure,
used in the same work, see the specs repo and
CodeOcean capsule
0541628.
CG is a set of scripts for extracting information about (localized) somatic
mutations in protein coding regions of cancer genomes from
TCGA
and
ICGC
databases. The icgc branch also contains tools necessary to
merge the info from the two sources. The scripts add up to a loosely connected pipeline
(or, rather, two pipelines). They re-organize the data and store it in the local MySQL
database. While the back end of each pipeline is rather generic, the front end is
geared toward answering particular questions for which they were originally written.
CG is not an out-of-the box solution. Rather, it is a starter kit, in case you would like to
do some cancer genomics data analysis on your own. Installing CG database(s) may take a
day or two if you are willing to go with the pipeline as-is. The (approximate) timings
for a default installation of the ICGC branch can be found in timing.txt.
With tweaks, a week is
not an unreasonable time estimate.
Why bother with the local copy of the data? I do not know a general answer to that question.
You should check the homepage of both databases - maybe the information you
are looking for is already available.
Here the original motive was to study the co-occurrence of mutations in different genes.
These types of questions are not readily answerable using data portals - portals tend to
agglomerate data on per-gene basis, in order to protect the privacy of sample donors.
In addition to TCGA and ICGC databases themselves, CG relies on
The tcga branch of the icgc pipeline got obsoleted before coming to production stage.
It contains various blind alleys and wrong turns. Its current use is as a prep
step for merging with the icgc branch. It has its own README page.
A general note: throughout the pipeline, you will find scripts disabled by having exit() right on the top of the file.
These are the scripts that drop tables and/or store without checking. Enable them by commenting the exit line.
(The advice is to put the comment back in once the script is done.)
A general note no. 2: all scripts that take an argument issue their usage statement.
You can set some recurring constants - such as data directories or mysql conf file(s) -
in the config.py file.
Just like the tcga branch, this branch of the pipeline starts by downloading the data from
the source, ICGC in this case. Note however that here you will need the access token.
Obtaining one is a lengthy process (as in several weeks to several months) which you can
start investigating here.
One you have the access token, place it in the environmental variable called ICGC_TOKEN, to make
the download scripts work.
Note in particular that we are grouping some cancers under the same head-group.
See 06_group_cancers.py. This is because different depositors may
use different shorthands for the same cancer (e.g. LICA == ‘Liver Cancer’,
LINC == ‘Liver Cancer - NCC’, LIRI == ‘Liver Cancer - RIKEN’), though in some cases it
might not be clear which cancer the depositors refer to. Feel free to change in you version of the code
the grouping defined in 06_group_cancers.py, or to skip it altogether.
Make sure you have
the mysql conf file, and set its path in config.py (or arrange some other way to
access the local database, by modifying connect_to_mysql() in icgc_utils/mysql.py).
The last I checked, python’s MySQLdb package did not work with
the encripted cnf files, so the only alternative is using
client field in generic mysql option file,
like this:
[client]user = blahhost = localhostpassword = "somepasswd"
In MySQL shell (or however you talk to your MySQL server) create the database icgc and the user blah,
and give blah the permissions to write to and read from icgc:
create database icgc;create user 'blah'@'localhost' identified by 'somepasswd';grant all privileges on icgc.* to 'blah'@'localhost';flush privileges;
05_find_max_field_length.py and
06_make_tables.py:
Make sure that the fields in the mysql tables are big enough
for each entry and create mysql tables.
05_find_max_field_length.py should
give you an idea about the longest entries found.
07_write_mutations_tsv.py through
10_make_indices.py.
For large tables, rather than loading them through python,
it turns out to be faster to create tsvs and then load them from mysql shell
(as in 09_load_mysql.py; alternative: use mysqlimport manually)
to read them in wholesale. These scripts take care of that part, plus some index creating on the newly loaded tables.
Make sure to run 10_make_indices_on_temp_tables.py,
10_reorganize_variants.py in 20_local_db_reorganization
pretty much does not work without it at all.
All index making is slow here (see timing.txt) - run overnight.
Some checks are thrown in here that may inform the rest of the pipeline.
15_assembly_check.py confirms that, as of v27, all ICGC entries
refer to GRCh37, and 16_donor_check.py highlights the fact that
some donor ids have no somatic mutations in ICGC. This is somewhat mysterious, because some refer to
TCGA donors with somatic mutation data available from TCGA archive.
In 20_hgnc_name_resolution_table.py and
22_ensembl_id_table.py we make and fill some tables we will use later for name resolution
(translating between gene and protein names used in different contexts).
The annotation across different submitters to TCGA/ICGC is not uniform. In particular, for the missense mutations
sometimes it is not clear which splice they refer to. Alternatively, when the reference splice(s) is listed, it
may not be clear which splice is the canonical splice. To remedy that, here we do some basic annotation on our own. For that
we will need the coding sequence of canonical transcripts.
The canonical transcript id is not readily available from Ensembl Mart. For our purposes
here you can find this info in the tables called ensembl_gene2trans_stable.tsv.bz2
and in ensembl_deprecated2new_id.tsv.bz2 the
hacks directory. Decompress them (bzip2 -d) and put someplace where
20_hgnc_name_resolution_table.py can find them.
From the same place (Ensembl Mart) all human coding sequences can be downloaded in fasta format.
You can reduce that file to sequences of canonical transcript only. This is left as na exercise for the reader
(hint: use ensembl_gene2trans_stable.tsv and
blastdbcmd
tool; blastdbcmd has batch mode (blastdbcmd -h)). When you are happy with your fasta file, put is somewhere
where 23_ensembl_coding_seqs.py can find it.
Note that this is optional: if you are happy without knowing the position of the missense mutant on
the canonical translation, you can move on without this step.
ICGC database does not have a complete consensus on location annotation, so we will be doing it ourselves.
As a prep, we download gene coordiantes from UCSC. (The coordinates are actually from Ensembl, but UCSC
keeps it in a format that is more readily usable.) The script is
24_ucsc_gene_coords_table.py. To download coordinates
from their MySQl server you will need an internet connection, mysql client, and another conf file, like this:
[client]skip-auto-rehashuser = genomehost = genome-mysql.soe.ucsc.eduport = 3306
Again, the path to that file is expected to be defined in the config.py file.
This is where we depart from ICGC original database architecture - which is pretty much
nonexistent and consists of massive duplication of annotation for each occurrence of a mutation
and for each of its interpretations within various transcripts.
So instead we reorganize the database into something like this
where _specimen, \_donor, and *_simplesomatic tables exist for each cancer type, and mutations\*
and locations_* tables exist for each chromosome.
In 06_consequence_vocab.py
we inspect the ‘consequence’ vocabulary employed by ICGC. There seems to
some confusion there about the location vs. the consequence of a mutation. This info is used
in 13_reorganize_mutations.py to
come up with the pathogenicity estimate, to be stored in the eponymous field in the
mutations* tables.
New tables are created in 08_check_mut_tables_and_make_new_ones.py.
Note 1: In 10_reorganize_variants.py,
11_delete_variants_from_normal.py,
13_reorganize_mutations.py, and
14_reorganize_locations.py
you can choose to run in parallel (the number of ‘chunks’ in main()).
Note 2: In 11_delete_variants_from_normal.py,
we are dropping variants from normal samples. This is something you might not want to do if you are trying to
annotate variants yourself. Though in that case you might want to go back to without
08_check_mut_tables_and_make_new_ones.py
and keep the ‘matched_icgc_sample_id’ field.
Note 3: doing things carefully leads to some interesting results. NACA (Nasopharyngeal cancer) set, for example,
consists of normal tissue samples only.
14_reorganize_locations.py script uses
UCSC gene annotation to check chromosome addresses. The only information we are looking for here is the
possibility that the location falls within the splice region just outside of an exon. Mutations at these positions
are annotated as (possibly) pathogenic
in 15_update_pathogenicity_in_mutation_tables.py.
In 16_pathg_from_mutations_to_variants.py we are adding this
same info (it is a boolean flag, no big elaboration on the table) to the variantis (*_simple_somatic) tables.
(Do not forget to create indices using 10_make_indices_on_temp_tables.py.)
We add a couple of values to each row to later make the search for meaningful entries faster
we are adding mutant_allele_read_count/total_read_count ratio and pathogenicity estimate (boolean)
to simple_somatic tables. These scripts can be run in the order indicated in the name, rather than at this point.
In the following script,
21_add_realiability_annotation_to_somatic.py,
we combine these two columns into a reliability estimate: a somatic mutation in individual patient is considered reliable if mutant_allele_read_count>=10
and mut_to_total_read_count_ratio>=0.2. Information about the mutation in general (mutations_chromosome tables;
27_reliability_from_variants_to_mutations.py)
is considered reliable if there is at least one patient for which it was reliably established.
ICGC is rife with data duplication, coming from various sources. Some seem to be bookkeeping mistakes with the
same patient data finding its way into the dataset through various depositors; some are the results of the re-sampling
of the same tumor, while some are completely obscure, with all identifiers being identical everywhere.
All identifiers identical.
18_cleanup_duplicate_entries.py:
Some mutations have identical tuple
of identifiers (icgc_mutation_id, icgc_donor_id, icgc_specimen_id, icgc_sample_id). Note that this
is after we have reorganized the database so that the mutation and location info sit in
different tables from the donor info. Not sure what this is about (the same sample analyzed independently multiple
times?), but when a duplicate is found, this script chooses the entry with the greatest depth reported.
See the script for the full resolution strategy and make sure to run
17_make_jumbo_index
beforehand.
Same submitted_sample_id, different donor ids.
There might be further problems: See for example, mutation MU2003689, which,
so the ICGC page claims can be found in two distinct donors.
The closer inspection of the two donors shows however that their submitter ID is the same, as is the age
of the ‘two’ women. (The tumour subtype has different description, reflecting, apparently, the
curator’s preference.) Indeed, donors table for BRCA, at this point in the pipeline has
1976 distinct ICGC donor ids, and 1928 distinct submitter IDs. BRCA does turn out to be the biggest offender here,
followed by LICA with 8 duplicated donors. It is not clear whether these duplicates refer to the same
tumor at the same stage because even the submitter sample ids might be different
(see 19_cleanup_multiple_donor_for_the_same_submitted_id.py).
Same donor with differing specimen and sample ids. Apparently ICGC refers to biological replicates - i.e.
samples taken from different sites - as specimens, and
to technical replicates as samples. Perhaps they might have a role when answering different
types of questions than what we have in mind. Here, however we do not want to have these results mistaken for recurring mutations,
thus we remove them in 22_cleanup_duplicate_specimens.py
and 23_cleanup_multiple_sample_ids.py, but not before checking
which of the samples produced more reliable reads (see below).
In this version of the pipeline we keep only the sample annotated as ‘Primary tumour - solid tissue.’ Out of
these, if multiple refer to the same submitter id, we keep the ones with the largest reported number of
somatic mutations. The investigation of the source of this duplication is again outside of our zone of interest.
Different donor and submitter ids with overlapping variants.
Unfortunately, the duplicates do not stop here. See for example DO224621
and DO230968 that have different donor and submitter
ids, and mysteriously have 1703 identical variants.
Both donors are males diagnosed with lung squamous carcinoma. One sample turns out to be WES, while the other is WGS.
25_cleanup_duplicate_donors_by_variant_overlap.py
attempts to detect such cases by looking for suspiciously high overlap in reported variants. When such case is detected,
the sample with the higher average depth of sampling is retained.
See the script for the criteria used.
29_index_on_mutation_tables.py
and 30_mutation2gene_maps.py
together build a table that maps icgc_mutation_id to HUGO gene
symbol, and vice versa.
It may be helpful to change the annotation of missense mutations
(34_reannotate_missense_mutations.py)
to include the currently accepted canonical transcript, according to Ensembl.
After running reference 34_reannotate_missense_mutations.py,
the reference will be lost to other (non-canonical) transcripts, which can be retrieved from the locations table.
The prerequisite is 32_load_ensembl_coord_patches.py,
a hacky solution to include the latest gene coordinates from Ensembl. The coordinate tables can be found in
hacks/coord_patches.tar.bz2 (rightclick to download)
A step toward independent annotation.
It is possible to stop here and move to production scripts, if there is no interest in including TCGA.
It might be advisable to backup the newly-created tables, by storing them as an sqldump, for example.
There is a couple of scripts in hacks directory to help along.
somatic_tables_dump.pl will dump them out, just make sure you move to the storage direcotry,
and somatic_tables_load.pl will load them back in if needs be.
If you choose to do the full database dump,
somatic_tables_from_dump.pl can extract only *_simple_somatic tables,
the process if slow, however.
The scripts 28_add_tables_for_tcga_variants.py
through 40_mutation2gene_maps.py
concern themselves with merging TCGA info created in TCGA branch with the ICGC.
This sub-pipe runs from creating tables, to checking for duplicates in the newly merged database.
Note thus, that some scripts are actually softlinks to the scripts we used in cleaning up the ICGC version.
If everything is ok, 35_database_stats.py should report
no donor or specimen duplicates in any of the tables. Samples might be duplicated, we take these are
bona fide from the same biopsy. The potential duplicates should have been resolved on the level of variants -
the entries from different samples are ok unless they report the exact same variant (location and replacement).
Assorted stats and checks, including the late-point decision to mark
variants from hypermutating cancers (by def, more than 1000 pathogenic mutations in coding regions)
as unreliable, 42_mark_hypermutators_as_unreliable.py.
The fun starts here, and you probably might want to do something else for fun. Here is what we looked into:
42_gene_stats_generic.py. Given an HGNC symbol for a protein coding gene,
it produces a table in tsv format of all pathogenic mutations found in the database. It also reports the mutation status
of a background gene (default is TP53). The table is created in the directory where the script is invoked.
44_mutation_freq_stats.py. Given an HGNC symbol (or two) for a protein coding gene,
it produces a table in tsv format containing the information about
the overall stats of the cancers in which the query gene appears mutated. This should give you an idea
of how extraordinary (or not) is the number of mutations from the previous script. Thus for example, if a gene is
mutated in 5% of all patients in the database, it might not be very significant if the same is true for 7,000 other
genes, as might be the case if the cancer at hand is mutation-happy. To the contrary, 2% of the patients might be meaningful if it holds
for a small percentage of the genome. The table is created in the directory where the script is invoked.
45_two_gene_co-ocurrence.py. Given a background gene and a list
of other gene names, the script reports how often the members of this list, taken as a group, appear mutated in the
same tumor sample as the background gene. The script also attempts to give a statistical evaluation of how (un)likely
would be the size of the overlap under random selection of tumors (or patients; in our database there is only
one specimen per patient).
Both simple Fisher test and mutation-propensity-corrected test results are reported.
To use the latter you need to compile the C code in c-utils/random_binsize_fisher.
A simple make will do, no need to move the executable, no further dependencies. Alternatively you can comment it out
in the python script (the invocation point is called ‘size_corrected_pvals_C’).
The C code is neutral on the issue of
how to measure the mutation propensity (actually it does not know anything about the mutation, it only samples
from a non-uniform distribution). Trial and error seems to indicate that assigning a measure that is directly
proportional to the number of mutations tends to significantly over-estimate the expected co-occurrence of mutations.
As of this writing, the question remains open.
48_co-ocurrence_of_mutations_and_lack_thereof.py.
This script will take a single gene and output some prep files for calculating the co-occurrence propensity of mutations
in this gene against all other protein coding genes in the human genome. It is supposed to work in pair with the following script.
50_co-ocurrence_postprocess.py. What
45_two_gene_co-ocurrence.py script does for a single gene,
this script repeats for one-against-all case: it reads the output from
48_co-ocurrence_of_mutations_and_lack_thereof.py
and adds Fisher and mutation-propensity-corrected stats for the co-occurrence of mutations.
This might take several hours for a highly mutated gene such as TP53.
56_silent_to_nonsilent_ratio.py.
Silent-to-nonsilent ratio of mutations for all protein coding genes in the database.
58_silent_to_nonsilent_plot.py.
Plot of the previous, using Matplotlib.
60_clustering_on_the_structure.py.
Given gene symbol, a protein structure in the PDB format and the name of the “background” gene,
the script produces statistical estimate of the (non)randomness of the mutation positions seen in cancers when mapped
on the protein structure. The script also produces input for PyMol.
The script relies on a piece of C-code, to be found in c-utils/clustering.
A simple make will do, no need to move the executable, no further dependencies.
Three pairs of output files are produced, one for mutations seen in the wt background gene, one for mutated background
gene, and one for all mutations, irrespective of the background status: output files that list the membership
of the clusters on the structure, and pml files, the PyMol visualization scripts.
If you are not familiar with PyMol, click here
for some pointers.
Note that you are on your own to procure the PDB structure
and assure that the residues there are not mutated and
their numbering corresponds to the canonical transcript of the gene.
Do not use locks in MySQL. Locks are evil. Good luck with the rest.