neksa / mutagene Goto Github PK
View Code? Open in Web Editor NEWPython library and package for mutational analysis
Home Page: https://www.ncbi.nlm.nih.gov/projects/mutagene/
License: Other
Python library and package for mutational analysis
Home Page: https://www.ncbi.nlm.nih.gov/projects/mutagene/
License: Other
the webserver and package produce different results:
(env_mutagene) Annas-MBP:mutagene panch$ mutagene signature -i TSVC_variants_IonXpress_019_(BCMDS-116).vcf -g hg38 -f VCF -s 10 -U
sample signature exposure mutations
VCF 2 0.0535 3
VCF 3 0.1164 7
VCF 4 0.1134 7
VCF 5 0.1244 7
VCF 6 0.3912 23
VCF 7 0.1603 9
VCF 8 0.0109 1
VCF 9 0.0274 2
(env_mutagene) Annas-MBP:mutagene panch$ mutagene signature -i TSVC_variants_IonXpress_019_(BCMDS-116).vcf -g hg38 -f VCF -s 10
sample signature exposure mutations
VCF 1 0.0144 1
VCF 3 0.1123 7
VCF 4 0.0514 3
VCF 5 0.1074 6
VCF 6 0.1604 9
VCF 7 0.0341 2
VCF 9 0.0328 2
$ mutagene -v rank -g hg38 -i TCGA.COAD.mutect.somatic.maf --profile TCGA.COAD.mutect.somatic.maf
INFO Loaded 230248 mutations
INFO Loaded 190240 protein mutations in 399 samples skipped 74543 protein mutations
INFO Cohort size: 399
INFO Profile overridden
INFO THRESHOLD_DRIVER: 8.030647e-05
INFO THRESHOLD_PASSENGER: 0.003440945
Traceback (most recent call last):
File "/Library/Developer/CommandLineTools/Library/Frameworks/Python3.framework/Versions/3.7/lib/python3.7/runpy.py", line 193, in _run_module_as_main
"main", mod_spec)
File "/Library/Developer/CommandLineTools/Library/Frameworks/Python3.framework/Versions/3.7/lib/python3.7/runpy.py", line 85, in _run_code
exec(code, run_globals)
File "/Users/panch/mutagene/env_mutagene/lib/python3.7/site-packages/mutagene/main.py", line 113, in
MutaGeneApp()
File "/Users/panch/mutagene/env_mutagene/lib/python3.7/site-packages/mutagene/main.py", line 104, in init
parser_class.callback(args)
File "/Users/panch/mutagene/env_mutagene/lib/python3.7/site-packages/mutagene/cli/rank_menu.py", line 114, in callback
rank(protein_mutations, args.outfile, profile, cohort_aa_mutations, cohort_size, args.threshold_driver, args.threshold_passenger)
File "/Users/panch/mutagene/env_mutagene/lib/python3.7/site-packages/mutagene/mutability/mutability.py", line 138, in rank
mutation_model = calculate_base_substitution_mutability(profile, cohort_size)
File "/Users/panch/mutagene/env_mutagene/lib/python3.7/site-packages/mutagene/mutability/mutability.py", line 111, in calculate_base_substitution_mutability
assert len(counts_profile) == 96
AssertionError
mutagene currently can't process VCFs due to updates in MAF reader
steps to reproduce error:
run: mutagene motif search -g hg19 -i sample2.vcf
caroline's draft code to fix:
in motif_menu.py:
if ".vcf" not in args.infile.name:
mutations, mutations_with_context, processing_stats = read_MAF_with_context_window(args.infile, args.genome, args.window_size)
if len(mutations_with_context) == 0:
logger.warning("No mutations loaded")
matching_motifs = identify_motifs(mutations_with_context, custom_motif,
args.strand) if mutations_with_context is not None else []
else:
mutations, mutations_with_context, processing_stats = read_VCF_with_context_window(args.infile, args.genome, args.window_size)
if len(mutations_with_context) == 0:
logger.warning("No mutations loaded")
matching_motifs = identify_motifs(mutations_with_context, custom_motif,
args.strand) if mutations_with_context is not None else []
in motifs/init.py:
try:
for sample, mutations in samples_mutations.items():
if mutations is not None and len(mutations) > 0:
first_mut_seq_with_coords = mutations[0][-1]
window_size = (len(first_mut_seq_with_coords) - 1) // 2
for m in search_motifs:
for s in strand:
# print("IDENTIFYING MOTIF: ", m['name'])
result = get_enrichment(mutations, m['motif'], m['position'], m['ref'], m['alt'], window_size, s)
debug_data = {'sample': sample, 'motif': m['logo'], 'strand': s}
debug_data.update(result)
debug_string = pprint.pformat(debug_data, indent=4)
logger.debug(debug_string)
if result['mutation_load'] == 0:
continue
motif_dict = {
'sample': sample,
'name': m['name'],
'motif': m['logo'],
'strand': s,
'enrichment': result['enrichment'],
'pvalue': result['pvalue_fisher'],
'mutations_low_est': result['mutation_load'],
'mutations_high_est': result['bases_mutated_in_motif'],
# 'mut_motif': result['bases_mutated_in_motif'],
# 'mut_not_in_motif': result['bases_mutated_not_in_motif'],
# 'stat_count': result['bases_not_mutated_in_motif'],
# 'ref_count': result['bases_not_mutated_not_in_motif']
}
motif_matches.append(motif_dict)
except AttributeError:
# print(sample, len(mutations))
if samples_mutations is not None and len(samples_mutations) > 0:
first_mut_seq_with_coords = samples_mutations[0][-1]
window_size = (len(first_mut_seq_with_coords) - 1) // 2
for m in search_motifs:
for s in strand:
# print("IDENTIFYING MOTIF: ", m['name'])
result = get_enrichment(samples_mutations, m['motif'], m['position'], m['ref'], m['alt'], window_size, s)
debug_data = {'sample': "VCF", 'motif': m['logo'], 'strand': s}
debug_data.update(result)
debug_string = pprint.pformat(debug_data, indent=4)
logger.debug(debug_string)
if result['mutation_load'] == 0:
continue
motif_dict = {
'sample': "VCF",
'name': m['name'],
'motif': m['logo'],
'strand': s,
'enrichment': result['enrichment'],
'pvalue': result['pvalue_fisher'],
'mutations_low_est': result['mutation_load'],
'mutations_high_est': result['bases_mutated_in_motif'],
# 'mut_motif': result['bases_mutated_in_motif'],
# 'mut_not_in_motif': result['bases_mutated_not_in_motif'],
# 'stat_count': result['bases_not_mutated_in_motif'],
# 'ref_count': result['bases_not_mutated_not_in_motif']
}
motif_matches.append(motif_dict)
return motif_matches
in context_window.py:
def read_VCF_with_context_window(muts, asm, window_size):
cn = complementary_nucleotide
mutations = defaultdict(float)
N_skipped = 0
# N_skipped_indels = 0
raw_mutations = []
for line in muts.read().split("\n"):
if line.startswith("#"):
continue
if len(line) < 10:
continue
col_list = line.split()
if len(col_list) < 4:
continue
# ID = col_list[2]
# chromosome is expected to be one or two number or one letter
chrom = col_list[0] # VCF CHROM
if chrom.lower().startswith("chr"):
chrom = chrom[3:]
# if len(chrom) == 2 and chrom[1] not in "0123456789":
# chrom = chrom[0]
pos = int(col_list[1]) # VCF POS
x = col_list[3] # VCF REF
y = col_list[4] # VCF ALT
# if multiple REF or ALT alleles are given, ignore mutation entry (could mean seq error, could mean deletion)
if len(x) != 1:
N_skipped += 1
continue
if len(y) != 1:
N_skipped += 1
continue
transcript_strand = '+' # bad assumption about transcript strand
raw_mutations.append((chrom, pos, transcript_strand, x, y))
# print("RAW", raw_mutations)
# print("INDELS", N_skipped)
mutations_with_context = []
if len(raw_mutations) > 0:
contexts = get_context_twobit_window(raw_mutations, asm, window_size)
if contexts is None or len(contexts) == 0:
return None, None
for (chrom, pos, transcript_strand, x, y) in raw_mutations:
(p5, p3), seq_with_coords = contexts.get((chrom, pos), (("N", "N"), []))
# print("RESULT: {} {}".format(p5, p3))
if len(set([p5, x, y, p3]) - set(nucleotides)) > 0:
# print(chrom, pos, p5, p3, x)
# print("Skipping invalid nucleotides")
N_skipped += 1
continue
if x in "CT":
mutations[p5 + p3 + x + y] += 1.0
else:
# complementary mutation
mutations[cn[p3] + cn[p5] + cn[x] + cn[y]] += 1.0
mutations_with_context.append((chrom, pos, transcript_strand, x, y, seq_with_coords))
N_loaded = int(sum(mutations.values()))
processing_stats = {'loaded': N_loaded, 'skipped': N_skipped, 'format': 'VCF'}
return mutations, mutations_with_context, processing_stats
not an informative warning message if no signature set is specified:
WARNING Set of signatures required. Use 5 and 10 for MUTAGENE-5 and MUTAGENE-10. Use 30 for COSMIC-30
Replace this message into - WARNING Use "-s" option to specify a set of signatures
No option to list available genomes to fetch like how cohorts can be listed.
mutagene -v rank -i test.maf -g hg38 -c gcb_lymphomas -p test.profile
WARNING Parsing MAF failed. Check that the input file is in MAF format or specify a different format using option -f
ValueError('Unexpected value of transcript_strand in MAF file')
mutagene report error when transcript_strand is not available in some lines of input MAF files.
Explore CircleCI integration with GitHub for use in CI testing.
Specifying a wrong format type (VCF, not MAF) produces an uninformative message:
(env_mutagene) Annas-MBP:mutagene panch$ mutagene signature identify -i data_mutations_extended.txt -g hg19 -s 49 -U -b -f VCF
Traceback (most recent call last):
File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/runpy.py", line 193, in _run_module_as_main
"main", mod_spec)
File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/runpy.py", line 85, in _run_code
exec(code, run_globals)
File "/Users/panch/mutagene/env_mutagene/lib/python3.7/site-packages/mutagene/main.py", line 113, in
MutaGeneApp()
File "/Users/panch/mutagene/env_mutagene/lib/python3.7/site-packages/mutagene/main.py", line 104, in init
parser_class.callback(args)
File "/Users/panch/mutagene/env_mutagene/lib/python3.7/site-packages/mutagene/cli/signature_menu.py", line 117, in callback
self.identify(args)
File "/Users/panch/mutagene/env_mutagene/lib/python3.7/site-packages/mutagene/cli/signature_menu.py", line 94, in identify
mutations, processing_stats = read_auto_profile(args.infile, fmt=args.input_format, asm=args.genome)
File "/Users/panch/mutagene/env_mutagene/lib/python3.7/site-packages/mutagene/io/mutations_profile.py", line 115, in read_auto_profile
mutations, processing_stats = read_VCF_profile(mutations_lines, asm)
File "/Users/panch/mutagene/env_mutagene/lib/python3.7/site-packages/mutagene/io/mutations_profile.py", line 234, in read_VCF_profile
pos = int(col_list[1]) # VCF POS
ValueError: invalid literal for int() with base 10: 'Entrez_Gene_Id'
Error parsing MAF file
Reproduce:
$ mutagene rank -i test.maf -g hg38 -o ranking.txt -c pancancer
WARNING MAF format not recognized
WARNING No mutations to rank
Sample:
Hugo_Symbol Entrez_Gene_Id Center NCBI_Build Chromosome Start_Position End_Position Strand Variant_Classification tumor_VAF Variant_Type Reference_Allele Tumor_Seq_Allele1 Tumor_Seq_Allele2 dbSNP_RS Tumor_Sample_Barcode Mutation_Status AAChange Transcript_Id TxChange GeneDetail.refGene
SAMD11 NCC hg38 chr1 930314 930314 + Missense_Mutation 0.5 SNP C C T NSLA_475_exonic_anno.txt Somatic p.H78Y NM_152486 c.C232T
...
MutaGene currently only reports 1 p-value - should functionality be added to let user choose which p-values are reported?
auto-generate package CLI and API documentation
Use readthedocs to host it
some mutation low estimates & all high estimates end in 0.5, should be whole number
error can be reproduced by running mutagens motif search on a sample
please copy the new documentation from "signature_doc*" file
specifically, PDF is generated but preview cannot open b/c "file may be damaged or use a file format that Preview doesn’t recognize", other pdf editors also will not open
issue will be reproduced if user runs -o file_name.pdf and tries to open PDF on their machine
formats that work: txt, csv, doc, tex, etc. seems to be a pdf issue
idea to fix:
It is not clear which motif is used by default. If "--motif" argument is missing, the program is hanging for long time and then returns an error. It should return a warning: "motif is not specified" in the beginning.
mutagene motif -i sample.maf
Traceback (most recent call last):
File "/cvmfs/soft.computecanada.ca/easybuild/software/2017/Core/python/3.7.4/lib/python3.7/runpy.py", line 193, in _run_module_as_main
"main", mod_spec)
File "/cvmfs/soft.computecanada.ca/easybuild/software/2017/Core/python/3.7.4/lib/python3.7/runpy.py", line 85, in _run_code
exec(code, run_globals)
File "/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/main.py", line 113, in
MutaGeneApp()
File "/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/main.py", line 104, in init
parser_class.callback(args)
File "/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/cli/motif_menu.py", line 149, in callback
MotifMenu.search(args)
File "/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/cli/motif_menu.py", line 127, in search
dump_matches=args.save_motif_matches) if mutations_with_context is not None else []
File "/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/motifs/init.py", line 46, in identify_motifs
motifs = get_known_motifs()
File "/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/io/motifs.py", line 21, in get_known_motifs
with open(fname, 'r') as f:
FileNotFoundError: [Errno 2] No such file or directory: '/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/io/../data/motifs/motifs.json'
(venv) [hpc4498@caclogin02 mutagene]$
(venv) [hpc4498@caclogin02 mutagene]$ mutagene motif -i TCGA.COAD.mutect.somatic.maf -w 250 -t 0.1 -o tmp1 -g hg38
Traceback (most recent call last):
File "/cvmfs/soft.computecanada.ca/easybuild/software/2017/Core/python/3.7.4/lib/python3.7/runpy.py", line 193, in _run_module_as_main
"main", mod_spec)
File "/cvmfs/soft.computecanada.ca/easybuild/software/2017/Core/python/3.7.4/lib/python3.7/runpy.py", line 85, in _run_code
exec(code, run_globals)
File "/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/main.py", line 113, in
MutaGeneApp()
File "/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/main.py", line 104, in init
parser_class.callback(args)
File "/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/cli/motif_menu.py", line 149, in callback
MotifMenu.search(args)
File "/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/cli/motif_menu.py", line 127, in search
dump_matches=args.save_motif_matches) if mutations_with_context is not None else []
File "/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/motifs/init.py", line 46, in identify_motifs
motifs = get_known_motifs()
File "/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/io/motifs.py", line 21, in get_known_motifs
with open(fname, 'r') as f:
FileNotFoundError: [Errno 2] No such file or directory: '/global/home/hpc4498/venv/lib/python3.7/site-packages/mutagene/io/../data/motifs/motifs.json'
When running mutagene signature identify, the signatures are reported as 1, 2, 3, etc., instead of Mutagene-1, Mutagene-2, ... or COSMIC-1, COSMIC-2, etc. This could cause confusion for user.
I performed ranking for about 100 MAF files but had the following Key errors for several of them. Do you have any ideas about this error? Thanks!
mutagene rank -i maf -g hg38 -o ranking.txt -p sample.profile
Traceback (most recent call last):
File "/usr/local/Anaconda/envs/py3.7/lib/python3.7/runpy.py", line 193, in _run_module_as_main
"main", mod_spec)
File "/usr/local/Anaconda/envs/py3.7/lib/python3.7/runpy.py", line 85, in _run_code
exec(code, run_globals)
File "/home/pengy10/.local/lib/python3.7/site-packages/mutagene/main.py", line 113, in
MutaGeneApp()
File "/home/pengy10/.local/lib/python3.7/site-packages/mutagene/main.py", line 104, in init
parser_class.callback(args)
File "/home/pengy10/.local/lib/python3.7/site-packages/mutagene/cli/rank_menu.py", line 75, in callback
protein_mutations, processing_stats = read_protein_mutations_MAF(args.infile, args.genome)
File "/home/pengy10/.local/lib/python3.7/site-packages/mutagene/io/protein_mutations_MAF.py", line 239, in read_protein_mutations_MAF
flat_mutations[(gene, protein_mutation)]['seq5'][props['seq5']] += 1
KeyError: 'GCTGA'
$ mutagene fetch cohorts --list
WARNING Looks like the file has been downloaded already: cohorts.tar.gz
There is no description in help files what does it mean "pre-identified motifs".
Specifying the wrong genome assembly returns an error about the maf format:
(env_mutagene) Annas-MacBook-Pro:mutagene panch$ mutagene -v rank -g hg18 -i TCGA.UCEC.mutect.somatic.maf
WARNING Parsing MAF failed. Check that the input file is in MAF format or specify a different format using option -f
FileNotFoundError(2, 'No such file or directory')
"mutagene rank" does not report an error if specifying genome assembly which does not match the one from maf file
currently requires user to manually input genome
ideas to fix:
Currently thresholds "-t" are not included in "advanced arguments". We should also include thresholds on individual p-values.
I don’t think it is clear what “TN”, “NA”… mean;
I would call these options “transcribed”,”non-transcribed”, “any” or “both” (“TN”?)
--save-motif-matches
“save mutations in matching motifs to a BED file”
There appear to be syntactic issues with motif_doc.rst along with formatting, numbering, and other issues. In addition, the documentation should be reviewed for accuracy and completeness.
need some sort of separation/clarification for optional arguments
user can't see distinction between these groups
Simplify rank sub package to calculate everything from one input file
"Identify" returns an error message for a small one sample file:
"Error processing the sample. The file could be too large for web processing - please try again. If error persists we recommend using stand-alone mutagene package. In some cases input format could not be parsed correctly."
Currently rank module does not support VCF
[Method] argument is not clearly described, there is some mix-up between target functions and optimization methods. Some "methods" which account or not account for context-independent processes contradict "--no-unexplained-variance" option.
Too many methods which are not necessary.
Log file must be generated for every program run (opt-in mode)
Log file should contain all necessary information to reproduce the results, in particular input parameters
It should be possible to read log file back as input configuration. Suggested to use YAML
Moved the last point from #27 "Include in the output a file with information on the arguments provided when running the command"
cohort args overlap, need clearer documentation + possible visuals/expl.
In mutagene rank module, we should specify in the output what was used as a profile, cohort size and observed mutations' count source.
if motif following -m is not in string format, program works unless user uses > to represent mutation
is user uses >, motif truncates beforehand
problem: having to manually put in string format is confusing to user, esp. b/c works sometimes w/o being a str
ways to fix:
Add a description of how to handle gene transcript in Ranking. Choose one longest transcript per gene.
Propose options for standardizing the output of the identify function. Specifically:
Related:
#15
example in "motif search --help" is not working:
mutagene motif --infile --input-format sample1.maf -f MAF --genome hg19 --motif 'C[A>T]' --strand A
It should be:
mutagene motif --infile sample1.maf -f MAF --genome hg19 --motif 'C[A>T]' --strand A
We should have an example with the significant motif matches.
mutagene signature identify -i TSVC_variants_IonXpress_019_(BCMDS-116).vcf -g hg38 -s 10 -f VCF
Produces different output with and without "-b" option, this issue is not observed for other VCF files
.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.