"""
Protein sequence alignment creation protocols/workflows.
Authors:
Thomas A. Hopf
Anna G. Green - complex protocol, hmm_build_and_search
Chan Kang - hmm_build_and_search
"""
from collections import OrderedDict, Iterable
import re
from shutil import copy
import os
import numpy as np
import pandas as pd
from evcouplings.align import tools as at
from evcouplings.align.alignment import (
detect_format, parse_header, read_fasta,
write_fasta, Alignment
)
from evcouplings.couplings.mapping import Segment
from evcouplings.utils import BailoutException
from evcouplings.utils.config import (
check_required, InvalidParameterError, MissingParameterError,
read_config_file, write_config_file
)
from evcouplings.utils.system import (
create_prefix_folders, get, valid_file,
verify_resources, ResourceError
)
from evcouplings.align.ena import (
extract_embl_annotation,
extract_cds_ids,
add_full_header
)
def _verify_sequence_id(sequence_id):
"""
Verify if a target sequence identifier is in proper
format for the pipeline to run without errors
(not none, and contains no whitespace)
Parameters
----------
id : str
Target sequence identifier to verify
Raises
------
InvalidParameterError
If sequence identifier is not valid
"""
if sequence_id is None:
raise InvalidParameterError(
"Target sequence identifier (sequence_id) must be defined and "
"cannot be None/null."
)
try:
if len(sequence_id.split()) != 1 or len(sequence_id) != len(sequence_id.strip()):
raise InvalidParameterError(
"Target sequence identifier (sequence_id) may not contain any "
"whitespace (spaces, tabs, ...)"
)
except AttributeError:
raise InvalidParameterError(
"Target sequence identifier (sequence_id) must be a string"
)
def _make_hmmsearch_raw_fasta(alignment_result, prefix):
"""
HMMsearch results do not contain the query sequence
so we must construct a raw_fasta file with the query
sequence as the first hit, to ensure proper numbering.
The search result is filtered to only contain the columns with
match states to the HMM, which has a one to one mapping to the
query sequence.
Paramters
---------
alignment_result : dict
Alignment result dictionary, output by run_hmmsearch
prefix : str
Prefix for file creation
Returns
-------
str
path to raw focus alignment file
"""
def _add_gaps_to_query(query_sequence_ali, ali):
# get the index of columns that do not contain match states (indicated by an x)
gap_index = [
i for i, x in enumerate(ali.annotation["GC"]["RF"]) if x != "x"
]
# get the index of columns that contain match states (indicated by an x)
match_index = [
i for i, x in enumerate(ali.annotation["GC"]["RF"]) if x == "x"
]
# ensure that the length of the match states
# match the length of the sequence
if len(match_index) != query_sequence_ali.L:
raise ValueError(
"HMMsearch result {} does not have a one-to-one"
" mapping to the query sequence columns".format(
alignment_result["raw_alignment_file"]
)
)
gapped_query_sequence = ""
seq = list(query_sequence_ali.matrix[0, :])
# loop through every position in the HMMsearch hits
for i in range(len(ali.annotation["GC"]["RF"])):
# if that position should be a gap, add a gap
if i in gap_index:
gapped_query_sequence += "-"
# if that position should be a letter, pop the next
# letter in the query sequence
else:
gapped_query_sequence += seq.pop(0)
new_sequence_ali = Alignment.from_dict({
query_sequence_ali.ids[0]: gapped_query_sequence
})
return new_sequence_ali
# open the sequence file
with open(alignment_result["target_sequence_file"]) as a:
query_sequence_ali = Alignment.from_file(a, format="fasta")
# if the provided alignment is empty, just return the target sequence
raw_focus_alignment_file = prefix + "_raw.fasta"
if not valid_file(alignment_result["raw_alignment_file"]):
# write the query sequence to a fasta file
with open(raw_focus_alignment_file, "w") as of:
query_sequence_ali.write(of)
# return as an alignment object
return raw_focus_alignment_file
# else, open the HMM search result
with open(alignment_result["raw_alignment_file"]) as a:
ali = Alignment.from_file(a, format="stockholm")
# make sure that the stockholm alignment contains the match annotation
if not ("GC" in ali.annotation and "RF" in ali.annotation["GC"]):
raise ValueError(
"Stockholm alignment {} missing RF"
" annotation of match states".format(alignment_result["raw_alignment_file"])
)
# add insertions to the query sequence in order to preserve correct
# numbering of match sequences
gapped_sequence_ali = _add_gaps_to_query(query_sequence_ali, ali)
# write a new alignment file with the query sequence as
# the first entry
with open(raw_focus_alignment_file, "w") as of:
gapped_sequence_ali.write(of)
ali.write(of)
return raw_focus_alignment_file
[docs]def fetch_sequence(sequence_id, sequence_file,
sequence_download_url, out_file):
"""
Fetch sequence either from database based on identifier, or from
input sequence file.
Parameters
----------
sequence_id : str
Identifier of sequence that should be retrieved
sequence_file : str
File containing sequence. If None, sqeuence will
be downloaded from sequence_download_url
sequence_download_url : str
URL from which to download missing sequence. Must
contain "{}" at the position where sequence ID will
be inserted into download URL (using str.format).
out_file : str
Output file in which sequence will be stored, if
sequence_file is not existing.
Returns
-------
str
Path of file with stored sequence (can be sequence_file
or out_file)
tuple (str, str)
Identifier of sequence as stored in file, and sequence
"""
if sequence_file is None:
get(
sequence_download_url.format(sequence_id),
out_file,
allow_redirects=True
)
else:
# if we have sequence file, try to copy it
try:
copy(sequence_file, out_file)
except FileNotFoundError:
raise ResourceError(
"sequence_file does not exist: {}".format(
sequence_file
)
)
# also make sure input file has something in it
verify_resources(
"Input sequence missing", out_file
)
with open(out_file) as f:
seq = next(read_fasta(f))
return out_file, seq
[docs]def cut_sequence(sequence, sequence_id, region=None, first_index=None, out_file=None):
"""
Cut a given sequence to sub-range and save it in a file
Parameters
----------
sequence : str
Full sequence that will be cut
sequence_id : str
Identifier of sequence, used to construct header
in output file
region : tuple(int, int), optional (default: None)
Region that will be cut out of full sequence.
If None, full sequence will be returned.
first_index : int, optional (default: None)
Define index of first position in sequence.
Will be set to 1 if None.
out_file : str, optional (default: None)
Save sequence in a FASTA file (header:
>sequence_id/start_region-end_region)
Returns
------
str
Subsequence contained in region
tuple(int, int)
Region. If no input region is given, this will be
(1, len(sequence)); otherwise, the input region is
returned.
Raises
------
InvalidParameterError
Upon invalid region specification (violating boundaries
of sequence)
"""
cut_seq = None
# (not using 1 as default value to allow parameter
# to be unspecified in config file)
if first_index is None:
first_index = 1
# last index is *inclusive*!
if region is None:
region = (first_index, first_index + len(sequence) - 1)
cut_seq = sequence
else:
start, end = region
str_start = start - first_index
str_end = end - first_index + 1
cut_seq = sequence[str_start:str_end]
# make sure bounds are valid given the sequence that we have
if str_start < 0 or str_end > len(sequence):
raise InvalidParameterError(
"Invalid sequence range: "
"region={} first_index={} len(sequence)={}".format(
region,
first_index,
len(sequence)
)
)
# save sequence to file
if out_file is not None:
with open(out_file, "w") as f:
header = "{}/{}-{}".format(sequence_id, *region)
write_fasta([(header, cut_seq)], f)
return region, cut_seq
[docs]def search_thresholds(use_bitscores, seq_threshold, domain_threshold, seq_len):
"""
Set homology search inclusion parameters.
HMMER hits get included in the HMM according to a two-step rule
1. sequence passes sequence-level treshold
2. domain passes domain-level threshold
Therefore, search thresholds are set based on the following logic:
1. If only sequence threshold is given, a MissingParameterException is raised
2. If only bitscore threshold is given, sequence threshold is set to the same
3. If both thresholds are given, they are according to defined values
Valid inputs for bitscore thresholds:
1. int or str: taken as absolute score threshold
2. float: taken as relative threshold (absolute threshold derived by
multiplication with domain length)
Valid inputs for integer thresholds:
1. int: Used as negative exponent, threshold will be set to 1E-<exponent>
2. float or str: Interpreted literally
Parameters
----------
use_bitscores : bool
Use bitscore threshold instead of E-value threshold
domain_threshold : str or int or float
Domain-level threshold. See rules above.
seq_threshold : str or int or float
Sequence-level threshold. See rules above.
seq_len : int
Length of sequence. Used to calculate absolute bitscore
threshold for relative bitscore thresholds.
Returns
-------
tuple (str, str)
Sequence- and domain-level thresholds ready to be fed into HMMER
"""
def transform_bitscore(x):
if isinstance(x, float):
# float: interpret as relative fraction of length
return "{:.1f}".format(x * seq_len)
else:
# otherwise interpret as absolute score
return str(x)
def transform_evalue(x):
if isinstance(x, int):
# if integer, interpret as negative exponent
return "1E{}".format(-x)
else:
# otherwise interpret literally
# (mantissa-exponent string or float)
return str(x).upper()
if domain_threshold is None:
raise MissingParameterError(
"domain_threshold must be explicitly defined "
"and may not be None/empty"
)
if use_bitscores:
transform = transform_bitscore
else:
transform = transform_evalue
if seq_threshold is not None:
seq_threshold = transform(seq_threshold)
if domain_threshold is not None:
domain_threshold = transform(domain_threshold)
# set "outer" sequence threshold so that it matches domain threshold
if domain_threshold is not None and seq_threshold is None:
seq_threshold = domain_threshold
return seq_threshold, domain_threshold
[docs]def describe_seq_identities(alignment, target_seq_index=0):
"""
Calculate sequence identities of any sequence
to target sequence and create result dataframe.
Parameters
----------
alignment : Alignment
Alignment for which description statistics
will be calculated
Returns
-------
pandas.DataFrame
Table giving the identity to target sequence
for each sequence in alignment (in order of
occurrence)
"""
id_to_query = alignment.identities_to(
alignment[target_seq_index]
)
return pd.DataFrame(
{"id": alignment.ids, "identity_to_query": id_to_query}
)
[docs]def describe_frequencies(alignment, first_index, target_seq_index=None):
"""
Get parameters of alignment such as gaps, coverage,
conservation and summarize.
Parameters
----------
alignment : Alignment
Alignment for which description statistics
will be calculated
first_index : int
Sequence index of first residue in target sequence
target_seq_index : int, optional (default: None)
If given, will add the symbol in the target sequence
into a separate column of the output table
Returns
-------
pandas.DataFrame
Table detailing conservation and symbol frequencies
for all positions in the alignment
"""
fi = alignment.frequencies
conservation = alignment.conservation()
fi_cols = {c: fi[:, i] for c, i in alignment.alphabet_map.items()}
if target_seq_index is not None:
target_seq = alignment[target_seq_index]
else:
target_seq = np.full((alignment.L), np.nan)
info = pd.DataFrame(
{
"i": range(first_index, first_index + alignment.L),
"A_i": target_seq,
"conservation": conservation,
**fi_cols
}
)
# reorder columns
info = info.loc[:, ["i", "A_i", "conservation"] + list(alignment.alphabet)]
return info
[docs]def describe_coverage(alignment, prefix, first_index, minimum_column_coverage):
"""
Produce "classical" buildali coverage statistics, i.e.
number of sequences, how many residues have too many gaps, etc.
Only to be applied to alignments focused around the
target sequence.
Parameters
----------
alignment : Alignment
Alignment for which coverage statistics will be calculated
prefix : str
Prefix of alignment file that will be stored as identifier in table
first_index : int
Sequence index of first position of target sequence
minimum_column_coverage : Iterable(float) or float
Minimum column coverage threshold(s) that will be tested
(creating one row for each threshold in output table).
.. note::
``int`` values given to this function instead of a float will be divided by 100 to create the corresponding
floating point representation. This parameter is 1.0 - maximum fraction of gaps per column.
Returns
-------
pd.DataFrame
Table with coverage statistics for different gap thresholds
"""
res = []
NO_MEFF = np.nan
if not isinstance(minimum_column_coverage, Iterable):
minimum_column_coverage = [minimum_column_coverage]
pos = np.arange(first_index, first_index + alignment.L)
f_gap = alignment.frequencies[:, alignment.alphabet_map[alignment._match_gap]]
for threshold in minimum_column_coverage:
if isinstance(threshold, int):
threshold /= 100
# all positions that have enough sequence information (i.e. little gaps),
# and their indeces
uppercase = f_gap <= 1 - threshold
uppercase_idx = np.nonzero(uppercase)[0]
# where does coverage of sequence by good alignment start and end?
cov_first_idx, cov_last_idx = uppercase_idx[0], uppercase_idx[-1]
# calculate indeces in sequence numbering space
first, last = pos[cov_first_idx], pos[cov_last_idx]
# how many lowercase positions in covered region?
num_lc_cov = np.sum(~uppercase[cov_first_idx:cov_last_idx + 1])
# total number of upper- and lowercase positions,
# and relative percentage
num_cov = uppercase.sum()
num_lc = (~uppercase).sum()
perc_cov = num_cov / len(uppercase)
res.append(
(prefix, threshold, alignment.N, alignment.L,
num_cov, num_lc, perc_cov, first, last,
last - first + 1, num_lc_cov, NO_MEFF)
)
df = pd.DataFrame(
res, columns=[
"prefix", "minimum_column_coverage", "num_seqs",
"seqlen", "num_cov", "num_lc", "perc_cov",
"1st_uc", "last_uc", "len_cov",
"num_lc_cov", "N_eff",
]
)
return df
[docs]def existing(**kwargs):
"""
Protocol:
Use external sequence alignment and extract all relevant
information from there (e.g. sequence, region, etc.),
then apply gap & fragment filtering as usual
Parameters
----------
Mandatory kwargs arguments:
See list below in code where calling check_required
Returns
-------
outcfg : dict
Output configuration of the pipeline, including
the following fields:
* sequence_id (passed through from input)
* alignment_file
* raw_focus_alignment_file
* statistics_file
* sequence_file
* first_index
* target_sequence_file
* annotation_file (None)
* frequencies_file
* identities_file
* focus_mode
* focus_sequence
* segments
"""
check_required(
kwargs,
[
"prefix", "input_alignment",
"sequence_id", "first_index",
"extract_annotation"
]
)
prefix = kwargs["prefix"]
# make sure output directory exists
create_prefix_folders(prefix)
# this file is starting point of pipeline;
# check if input alignment actually exists
input_alignment = kwargs["input_alignment"]
verify_resources(
"Input alignment does not exist",
input_alignment
)
# first try to autodetect format of alignment
with open(input_alignment) as f:
format = detect_format(f)
if format is None:
raise InvalidParameterError(
"Format of input alignment {} could not be "
"automatically detected.".format(
input_alignment
)
)
with open(input_alignment) as f:
ali_raw = Alignment.from_file(f, format)
# save annotation in sequence headers (species etc.)
annotation_file = None
if kwargs["extract_annotation"]:
annotation_file = prefix + "_annotation.csv"
from_anno_line = (format == "stockholm")
annotation = extract_header_annotation(
ali_raw, from_annotation=from_anno_line
)
annotation.to_csv(annotation_file, index=False)
# Target sequence of alignment
sequence_id = kwargs["sequence_id"]
# check if sequence identifier is valid
_verify_sequence_id(sequence_id)
# First, find focus sequence in alignment
focus_index = None
for i, id_ in enumerate(ali_raw.ids):
if id_.startswith(sequence_id):
focus_index = i
break
# if we didn't find it, cannot continue
if focus_index is None:
raise InvalidParameterError(
"Target sequence {} could not be found in alignment"
.format(sequence_id)
)
# identify what columns (non-gap) to keep for focus
focus_seq = ali_raw[focus_index]
focus_cols = np.array(
[c not in [ali_raw._match_gap, ali_raw._insert_gap] for c in focus_seq]
)
# extract focus alignment
focus_ali = ali_raw.select(columns=focus_cols)
focus_seq_nogap = "".join(focus_ali[focus_index])
# determine region of sequence. If first_index is given,
# use that in any case, otherwise try to autodetect
full_focus_header = ali_raw.ids[focus_index]
focus_id = full_focus_header.split()[0]
# try to extract region from sequence header
id_, region_start, region_end = parse_header(focus_id)
# override with first_index if given
if kwargs["first_index"] is not None:
region_start = kwargs["first_index"]
region_end = region_start + len(focus_seq_nogap) - 1
if region_start is None or region_end is None:
raise InvalidParameterError(
"Could not extract region information " +
"from sequence header {} ".format(full_focus_header) +
"and first_index parameter is not given."
)
# resubstitute full sequence ID from identifier
# and region information
header = "{}/{}-{}".format(
id_, region_start, region_end
)
focus_ali.ids[focus_index] = header
# write target sequence to file
target_sequence_file = prefix + ".fa"
with open(target_sequence_file, "w") as f:
write_fasta(
[(header, focus_seq_nogap)], f
)
# apply sequence identity and fragment filters,
# and gap threshold
mod_outcfg, ali = modify_alignment(
focus_ali, focus_index, id_, region_start, **kwargs
)
# generate output configuration of protocol
outcfg = {
**mod_outcfg,
"sequence_id": sequence_id,
"sequence_file": target_sequence_file,
"first_index": region_start,
"target_sequence_file": target_sequence_file,
"focus_sequence": header,
"focus_mode": True,
}
if annotation_file is not None:
outcfg["annotation_file"] = annotation_file
# dump config to YAML file for debugging/logging
write_config_file(prefix + ".align_existing.outcfg", outcfg)
# return results of protocol
return outcfg
[docs]def modify_alignment(focus_ali, target_seq_index, target_seq_id, region_start, **kwargs):
"""
Apply pairwise identity filtering, fragment filtering, and exclusion
of columns with too many gaps to a sequence alignment. Also generates
files describing properties of the alignment such as frequency distributions,
conservation, and "old-style" alignment statistics files.
.. note::
assumes focus alignment (otherwise unprocessed) as input.
.. todo::
come up with something more clever to filter fragments than fixed width
(e.g. use 95% quantile of length distribution as reference point)
Parameters
----------
focus_ali : Alignment
Focus-mode input alignment
target_seq_index : int
Index of target sequence in alignment
target_seq_id : str
Identifier of target sequence (without range)
region_start : int
Index of first sequence position in target sequence
kwargs : See required arguments in source code
Returns
-------
outcfg : Dict
File products generated by the function:
* alignment_file
* statistics_file
* frequencies_file
* identities_file
* raw_focus_alignment_file
ali : Alignment
Final processed alignment
"""
check_required(
kwargs,
[
"prefix", "seqid_filter", "hhfilter",
"minimum_sequence_coverage", "minimum_column_coverage",
"compute_num_effective_seqs", "theta",
]
)
prefix = kwargs["prefix"]
create_prefix_folders(prefix)
focus_fasta_file = prefix + "_raw_focus.fasta"
outcfg = {
"alignment_file": prefix + ".a2m",
"statistics_file": prefix + "_alignment_statistics.csv",
"frequencies_file": prefix + "_frequencies.csv",
"identities_file": prefix + "_identities.csv",
"raw_focus_alignment_file": focus_fasta_file,
}
# swap target sequence to first position if it is not
# the first sequence in alignment;
# this is particularly important for hhfilter run
# because target sequence might otherwise be filtered out
if target_seq_index != 0:
indices = np.arange(0, len(focus_ali))
indices[0] = target_seq_index
indices[target_seq_index] = 0
target_seq_index = 0
focus_ali = focus_ali.select(sequences=indices)
with open(focus_fasta_file, "w") as f:
focus_ali.write(f, "fasta")
# apply pairwise identity filter (using hhfilter)
if kwargs["seqid_filter"] is not None:
filtered_file = prefix + "_filtered.a3m"
at.run_hhfilter(
focus_fasta_file, filtered_file,
threshold=kwargs["seqid_filter"],
columns="first", binary=kwargs["hhfilter"]
)
with open(filtered_file) as f:
focus_ali = Alignment.from_file(f, "a3m")
# final FASTA alignment before applying A2M format modifications
filtered_fasta_file = prefix + "_raw_focus_filtered.fasta"
with open(filtered_fasta_file, "w") as f:
focus_ali.write(f, "fasta")
ali = focus_ali
# filter fragments
# come up with something more clever here than fixed width
# (e.g. use 95% quantile of length distribution as reference point)
min_cov = kwargs["minimum_sequence_coverage"]
if min_cov is not None:
if isinstance(min_cov, int):
min_cov /= 100
keep_seqs = (1 - ali.count("-", axis="seq")) >= min_cov
ali = ali.select(sequences=keep_seqs)
# Calculate frequencies, conservation and identity to query
# on final alignment (except for lowercase modification)
# Note: running hhfilter might cause a loss of the target seque
# if it is not the first sequence in the file! To be sure that
# nothing goes wrong, target_seq_index should always be 0.
describe_seq_identities(
ali, target_seq_index=target_seq_index
).to_csv(
outcfg["identities_file"], float_format="%.3f", index=False
)
describe_frequencies(
ali, region_start, target_seq_index=target_seq_index
).to_csv(
outcfg["frequencies_file"], float_format="%.3f", index=False
)
coverage_stats = describe_coverage(
ali, prefix, region_start, kwargs["minimum_column_coverage"]
)
# keep list of uppercase sequence positions in alignment
pos_list = np.arange(region_start, region_start + ali.L, dtype="int32")
# Make columns with too many gaps lowercase
min_col_cov = kwargs["minimum_column_coverage"]
if min_col_cov is not None:
if isinstance(min_col_cov, int):
min_col_cov /= 100
lc_cols = ali.count(ali._match_gap, axis="pos") > 1 - min_col_cov
ali = ali.lowercase_columns(lc_cols)
# if we remove columns, we have to update list of positions
pos_list = pos_list[~lc_cols]
else:
lc_cols = None
# compute effective number of sequences
# (this is intended for cases where coupling stage is
# not run, but this number is wanted nonetheless)
if kwargs["compute_num_effective_seqs"]:
# make sure we only compute N_eff on the columns
# that would be used for model inference, dispose
# the rest
if lc_cols is None:
cut_ali = ali
else:
cut_ali = ali.select(columns=~lc_cols)
# compute sequence weights
cut_ali.set_weights(kwargs["theta"])
# N_eff := sum of all sequence weights
n_eff = float(cut_ali.weights.sum())
# patch into coverage statistics (N_eff column)
coverage_stats.loc[:, "N_eff"] = n_eff
# create table with number of cluster members (inverse sequence
# weights) for each sequence
inv_seq_weights = pd.DataFrame({
"id": cut_ali.ids,
"num_cluster_members": cut_ali.num_cluster_members
})
# save sequence weights to file and add to output config
outcfg["sequence_weights_file"] = prefix + "_inverse_sequence_weights.csv"
inv_seq_weights.to_csv(
outcfg["sequence_weights_file"], index=False
)
else:
n_eff = None
# save coverage statistics to file
coverage_stats.to_csv(
outcfg["statistics_file"], float_format="%.3f",
index=False
)
# store description of final sequence alignment in outcfg
# (note these parameters will be updated by couplings protocol)
outcfg.update(
{
"num_sites": len(pos_list),
"num_sequences": len(ali),
"effective_sequences": n_eff,
"region_start": region_start,
}
)
# create segment in outcfg
outcfg["segments"] = [
Segment(
"aa", target_seq_id, region_start, region_start + ali.L - 1, pos_list
).to_list()
]
with open(outcfg["alignment_file"], "w") as f:
ali.write(f, "fasta")
return outcfg, ali
[docs]def jackhmmer_search(**kwargs):
"""
Protocol:
Iterative jackhmmer search against a sequence database.
Parameters
----------
Mandatory kwargs arguments:
See list below in code where calling check_required
.. todo::
explain meaning of parameters in detail.
Returns
-------
outcfg : dict
Output configuration of the protocol, including
the following fields:
* sequence_id (passed through from input)
* first_index (passed through from input)
* target_sequence_file
* sequence_file
* raw_alignment_file
* hittable_file
* focus_mode
* focus_sequence
* segments
"""
check_required(
kwargs,
[
"prefix", "sequence_id", "sequence_file",
"sequence_download_url", "region", "first_index",
"use_bitscores", "domain_threshold", "sequence_threshold",
"database", "iterations", "cpu", "nobias", "reuse_alignment",
"checkpoints_hmm", "checkpoints_ali", "jackhmmer",
"extract_annotation"
]
)
prefix = kwargs["prefix"]
# check if sequence identifier is valid
_verify_sequence_id(kwargs["sequence_id"])
# make sure output directory exists
create_prefix_folders(prefix)
# store search sequence file here
target_sequence_file = prefix + ".fa"
full_sequence_file = prefix + "_full.fa"
# make sure search sequence is defined and load it
full_seq_file, (full_seq_id, full_seq) = fetch_sequence(
kwargs["sequence_id"],
kwargs["sequence_file"],
kwargs["sequence_download_url"],
full_sequence_file
)
# cut sequence to target region and save in sequence_file
# (this is the main sequence file used downstream)
(region_start, region_end), cut_seq = cut_sequence(
full_seq,
kwargs["sequence_id"],
kwargs["region"],
kwargs["first_index"],
target_sequence_file
)
# run jackhmmer... allow to reuse pre-exisiting
# Stockholm alignment file here
ali_outcfg_file = prefix + ".align_jackhmmer_search.outcfg"
# determine if to rerun, only possible if previous results
# were stored in ali_outcfg_file
if kwargs["reuse_alignment"] and valid_file(ali_outcfg_file):
ali = read_config_file(ali_outcfg_file)
# check if the alignment file itself is also there
verify_resources(
"Tried to reuse alignment, but empty or "
"does not exist",
ali["alignment"], ali["domtblout"]
)
else:
# otherwise, we have to run the alignment
# modify search thresholds to be suitable for jackhmmer
seq_threshold, domain_threshold = search_thresholds(
kwargs["use_bitscores"],
kwargs["sequence_threshold"],
kwargs["domain_threshold"],
len(cut_seq)
)
# run search process
ali = at.run_jackhmmer(
query=target_sequence_file,
database=kwargs[kwargs["database"]],
prefix=prefix,
use_bitscores=kwargs["use_bitscores"],
domain_threshold=domain_threshold,
seq_threshold=seq_threshold,
iterations=kwargs["iterations"],
nobias=kwargs["nobias"],
cpu=kwargs["cpu"],
checkpoints_hmm=kwargs["checkpoints_hmm"],
checkpoints_ali=kwargs["checkpoints_ali"],
binary=kwargs["jackhmmer"],
)
# get rid of huge stdout log file immediately
# (do not use /dev/null option of jackhmmer function
# to make no assumption about operating system)
try:
os.remove(ali.output)
except OSError:
pass
# turn namedtuple into dictionary to make
# restarting code nicer
ali = dict(ali._asdict())
# save results of search for possible restart
write_config_file(ali_outcfg_file, ali)
# prepare output dictionary with result files
outcfg = {
"sequence_id": kwargs["sequence_id"],
"target_sequence_file": target_sequence_file,
"sequence_file": full_sequence_file,
"first_index": kwargs["first_index"],
"focus_mode": True,
"raw_alignment_file": ali["alignment"],
"hittable_file": ali["domtblout"],
}
# define a single protein segment based on target sequence
outcfg["segments"] = [
Segment(
"aa", kwargs["sequence_id"],
region_start, region_end,
range(region_start, region_end + 1)
).to_list()
]
outcfg["focus_sequence"] = "{}/{}-{}".format(
kwargs["sequence_id"], region_start, region_end
)
return outcfg
[docs]def hmmbuild_and_search(**kwargs):
"""
Protocol:
Build HMM from sequence alignment using hmmbuild and
search against a sequence database using hmmsearch.
Parameters
----------
Mandatory kwargs arguments:
See list below in code where calling check_required
Returns
-------
outcfg : dict
Output configuration of the protocol, including
the following fields:
* target_sequence_file
* sequence_file
* raw_alignment_file
* hittable_file
* focus_mode
* focus_sequence
* segments
"""
def _format_alignment_for_hmmbuild(input_alignment_file, **kwargs):
# this file is starting point of pipeline;
# check if input alignment actually exists
verify_resources(
"Input alignment does not exist",
input_alignment_file
)
# first try to autodetect format of alignment
with open(input_alignment_file) as f:
format = detect_format(f)
if format is None:
raise InvalidParameterError(
"Format of input alignment {} could not be "
"automatically detected.".format(
input_alignment_file
)
)
with open(input_alignment_file) as f:
ali_raw = Alignment.from_file(f, format)
# Target sequence of alignment
sequence_id = kwargs["sequence_id"]
if sequence_id is None:
raise InvalidParameterError(
"Parameter sequence_id must be defined"
)
# First, find focus sequence in alignment
focus_index = None
for i, id_ in enumerate(ali_raw.ids):
if id_.startswith(sequence_id):
focus_index = i
break
# if we didn't find it, cannot continue
if focus_index is None:
raise InvalidParameterError(
"Target sequence {} could not be found in alignment"
.format(sequence_id)
)
# identify what columns (non-gap) to keep for focus
# this should be all columns in the raw_focus_alignment_file
# but checking anyway
focus_seq = ali_raw[focus_index]
focus_cols = np.array(
[c not in [ali_raw._match_gap, ali_raw._insert_gap] for c in focus_seq]
)
# extract focus alignment
focus_ali = ali_raw.select(columns=focus_cols)
focus_seq_nogap = "".join(focus_ali[focus_index])
# determine region of sequence. If first_index is given,
# use that in any case, otherwise try to autodetect
full_focus_header = ali_raw.ids[focus_index]
focus_id = full_focus_header.split()[0]
# try to extract region from sequence header
id_, region_start, region_end = parse_header(focus_id)
# override with first_index if given (but respect region from alignment if defined)
if kwargs["first_index"] is not None and (region_start is None or region_end is None):
region_start = kwargs["first_index"]
region_end = region_start + len(focus_seq_nogap) - 1
if region_start is None or region_end is None:
raise InvalidParameterError(
"Could not extract region information " +
"from sequence header {} ".format(full_focus_header) +
"and first_index parameter is not given."
)
# resubstitute full sequence ID from identifier
# and region information
header = "{}/{}-{}".format(
id_, region_start, region_end
)
focus_ali.ids[focus_index] = header
# write target sequence to file
target_sequence_file = prefix + ".fa"
with open(target_sequence_file, "w") as f:
write_fasta(
[(header, focus_seq_nogap)], f
)
# swap target sequence to first position if it is not
# the first sequence in alignment;
# this is particularly important for hhfilter run
# because target sequence might otherwise be filtered out
if focus_index != 0:
indices = np.arange(0, len(focus_ali))
indices[0] = focus_index
indices[focus_index] = 0
focus_index = 0
focus_ali = focus_ali.select(sequences=indices)
# write the raw focus alignment for hmmbuild
focus_fasta_file = prefix + "_raw_focus_input.fasta"
with open(focus_fasta_file, "w") as f:
focus_ali.write(f, "fasta")
return focus_fasta_file, target_sequence_file, region_start, region_end
# define the gap threshold for inclusion in HMM's build by HMMbuild.
SYMFRAC_HMMBUILD = 0.0
# check for required options
check_required(
kwargs,
[
"prefix", "sequence_id", "alignment_file",
"use_bitscores", "domain_threshold", "sequence_threshold",
"database", "cpu", "nobias", "reuse_alignment",
"hmmbuild", "hmmsearch"
]
)
prefix = kwargs["prefix"]
# check if sequence identifier is valid
_verify_sequence_id(kwargs["sequence_id"])
# make sure output directory exists
create_prefix_folders(prefix)
# prepare input alignment for hmmbuild
focus_fasta_file, target_sequence_file, region_start, region_end = \
_format_alignment_for_hmmbuild(
kwargs["alignment_file"], **kwargs
)
# run hmmbuild_and_search... allow to reuse pre-exisiting
# Stockholm alignment file here
ali_outcfg_file = prefix + ".align_hmmbuild_and_search.outcfg"
# determine if to rerun, only possible if previous results
# were stored in ali_outcfg_file
if kwargs["reuse_alignment"] and valid_file(ali_outcfg_file):
ali = read_config_file(ali_outcfg_file)
# check if the alignment file itself is also there
verify_resources(
"Tried to reuse alignment, but empty or "
"does not exist",
ali["alignment"], ali["domtblout"]
)
else:
# otherwise, we have to run the alignment
# modify search thresholds to be suitable for hmmsearch
sequence_length = region_end - region_start + 1
seq_threshold, domain_threshold = search_thresholds(
kwargs["use_bitscores"],
kwargs["sequence_threshold"],
kwargs["domain_threshold"],
sequence_length
)
# create the hmm
hmmbuild_result = at.run_hmmbuild(
alignment_file=focus_fasta_file,
prefix=prefix,
symfrac=SYMFRAC_HMMBUILD,
cpu=kwargs["cpu"],
binary=kwargs["hmmbuild"],
)
hmmfile = hmmbuild_result.hmmfile
# run the alignment from the hmm
ali = at.run_hmmsearch(
hmmfile=hmmfile,
database=kwargs[kwargs["database"]],
prefix=prefix,
use_bitscores=kwargs["use_bitscores"],
domain_threshold=domain_threshold,
seq_threshold=seq_threshold,
nobias=kwargs["nobias"],
cpu=kwargs["cpu"],
binary=kwargs["hmmsearch"],
)
# get rid of huge stdout log file immediately
try:
os.remove(ali.output)
except OSError:
pass
# turn namedtuple into dictionary to make
# restarting code nicer
ali = dict(ali._asdict())
# only item from hmmsearch_result to save is the hmmfile
ali["hmmfile"] = hmmfile
# save results of search for possible restart
write_config_file(ali_outcfg_file, ali)
# prepare output dictionary with result files
outcfg = {
"sequence_file": target_sequence_file,
"first_index": region_start,
"input_raw_focus_alignment": focus_fasta_file,
"target_sequence_file": target_sequence_file,
"focus_mode": True,
"raw_alignment_file": ali["alignment"],
"hittable_file": ali["domtblout"],
}
# convert the raw output alignment to fasta format
# and add the appropriate query sequecne
raw_focus_alignment_file = _make_hmmsearch_raw_fasta(outcfg, prefix)
outcfg["raw_focus_alignment_file"] = raw_focus_alignment_file
# define a single protein segment based on target sequence
outcfg["segments"] = [
Segment(
"aa", kwargs["sequence_id"],
region_start, region_end,
range(region_start, region_end + 1)
).to_list()
]
outcfg["focus_sequence"] = "{}/{}-{}".format(
kwargs["sequence_id"], region_start, region_end
)
return outcfg
[docs]def standard(**kwargs):
"""
Protocol:
Standard buildali4 workflow (run iterative jackhmmer
search against sequence database, than determine which
sequences and columns to include in the calculation based
on coverage and maximum gap thresholds).
Parameters
----------
Mandatory kwargs arguments:
See list below in code where calling check_required
Returns
-------
outcfg : dict
Output configuration of the pipeline, including
the following fields:
* sequence_id (passed through from input)
* first_index (passed through from input)
* alignment_file
* raw_alignment_file
* raw_focus_alignment_file
* statistics_file
* target_sequence_file
* sequence_file
* annotation_file
* frequencies_file
* identities_file
* hittable_file
* focus_mode
* focus_sequence
* segments
ali : Alignment
Final sequence alignment
"""
check_required(
kwargs,
[
"prefix", "extract_annotation",
]
)
prefix = kwargs["prefix"]
# make sure output directory exists
create_prefix_folders(prefix)
# first step of protocol is to get alignment using
# jackhmmer; initialize output configuration with
# results of this search
jackhmmer_outcfg = jackhmmer_search(**kwargs)
stockholm_file = jackhmmer_outcfg["raw_alignment_file"]
segment = Segment.from_list(jackhmmer_outcfg["segments"][0])
target_seq_id = segment.sequence_id
region_start = segment.region_start
region_end = segment.region_end
# read in stockholm format (with full annotation)
with open(stockholm_file) as a:
ali_raw = Alignment.from_file(a, "stockholm")
# and store as FASTA file first (disabled for now
# since equivalent information easily be obtained
# from Stockholm file
"""
ali_raw_fasta_file = prefix + "_raw.fasta"
with open(ali_raw_fasta_file, "w") as f:
ali_raw.write(f, "fasta")
"""
# save annotation in sequence headers (species etc.)
if kwargs["extract_annotation"]:
annotation_file = prefix + "_annotation.csv"
annotation = extract_header_annotation(ali_raw)
annotation.to_csv(annotation_file, index=False)
# center alignment around focus/search sequence
focus_cols = np.array([c != "-" for c in ali_raw[0]])
focus_ali = ali_raw.select(columns=focus_cols)
target_seq_index = 0
mod_outcfg, ali = modify_alignment(
focus_ali, target_seq_index, target_seq_id, region_start, **kwargs
)
# merge results of jackhmmer_search and modify_alignment stage
outcfg = {
**jackhmmer_outcfg,
**mod_outcfg,
"annotation_file": annotation_file
}
# dump output config to YAML file for debugging/logging
write_config_file(prefix + ".align_standard.outcfg", outcfg)
if len(ali) <= 1:
raise BailoutException("align: No sequences found")
# return results of protocol
return outcfg
[docs]def complex(**kwargs):
"""
Protocol:
Run monomer alignment protocol and postprocess it for
EVcomplex calculations
Parameters
----------
Mandatory kwargs arguments:
See list below in code where calling check_required
Returns
-------
outcfg : dict
Output configuration of the alignment protocol, and
the following additional field:
genome_location_file : path to file containing
the genomic locations for CDs's corresponding to
identifiers in the alignment.
"""
check_required(
kwargs,
[
"prefix", "alignment_protocol",
"uniprot_to_embl_table",
"ena_genome_location_table"
]
)
verify_resources(
"Uniprot to EMBL mapping table does not exist",
kwargs["uniprot_to_embl_table"]
)
verify_resources(
"ENA genome location table does not exist",
kwargs["ena_genome_location_table"]
)
prefix = kwargs["prefix"]
# make sure output directory exists
create_prefix_folders(prefix)
# run the regular alignment protocol
# (standard, existing, ...)
alignment_protocol = kwargs["alignment_protocol"]
if alignment_protocol not in PROTOCOLS:
raise InvalidParameterError(
"Invalid choice for alignment protocol: {}".format(
alignment_protocol
)
)
outcfg = PROTOCOLS[kwargs["alignment_protocol"]](**kwargs)
# if the user selected the existing alignment protocol
# they can supply an input annotation file
# which overwrites the annotation file generated by the existing protocol
if alignment_protocol == "existing":
check_required(kwargs, ["override_annotation_file"])
if kwargs["override_annotation_file"] is not None:
verify_resources(
"Override annotation file does not exist",
kwargs["override_annotation_file"]
)
outcfg["annotation_file"] = prefix + "_annotation.csv"
annotation_data = pd.read_csv(kwargs["override_annotation_file"])
annotation_data.to_csv(outcfg["annotation_file"])
# extract cds identifiers for alignment uniprot IDs
cds_ids = extract_cds_ids(
outcfg["alignment_file"],
kwargs["uniprot_to_embl_table"]
)
# extract genome location information from ENA
genome_location_filename = prefix + "_genome_location.csv"
genome_location_table = extract_embl_annotation(
cds_ids,
kwargs["ena_genome_location_table"],
genome_location_filename
)
genome_location_table = add_full_header(
genome_location_table, outcfg["alignment_file"]
)
genome_location_table.to_csv(genome_location_filename)
outcfg["genome_location_file"] = genome_location_filename
# dump output config to YAML file for debugging/logging
write_config_file(prefix + ".align_complex.outcfg", outcfg)
return outcfg
# list of available alignment protocols
PROTOCOLS = {
# standard buildali protocol (iterative hmmer search)
"standard": standard,
# build raw multiple sequence alignment using jackmmer
"jackhmmer_search": jackhmmer_search,
# start from an existing (external) alignment
"existing": existing,
# run alignment protocol and postprocess output for
# complex pipeline
"complex": complex,
}
[docs]def run(**kwargs):
"""
Run alignment protocol to generate multiple sequence
alignment from input sequence.
Parameters
----------
Mandatory kwargs arguments:
protocol: Alignment protocol to run
prefix: Output prefix for all generated files
Optional:
Returns
-------
Alignment
Dictionary with results of stage in following fields (in brackets - not returned by all protocols):
* alignment_file
* [raw_alignment_file]
* statistics_file
* target_sequence_file
* sequence_file
* [annotation_file]
* frequencies_file
* identities_file
* [hittable_file]
* focus_mode
* focus_sequence
* segments
"""
check_required(kwargs, ["protocol"])
if kwargs["protocol"] not in PROTOCOLS:
raise InvalidParameterError(
"Invalid protocol selection: " +
"{}. Valid protocols are: {}".format(
kwargs["protocol"], ", ".join(PROTOCOLS.keys())
)
)
return PROTOCOLS[kwargs["protocol"]](**kwargs)