function tripal_get_feature_sequences
| 2.x tripal_feature.api.inc | tripal_get_feature_sequences($feature, $options) |
| 3.x tripal_chado.module.DEPRECATED.api.inc | tripal_get_feature_sequences($feature, $options) |
Retrieves the sequences for a given feature.
If a feature has multiple alignments or multiple relationships then multiple sequences will be returned.
Parameters
$feature: An associative array describing the feature. Valid keys include:
- feature_id: The feature_id of the feature for which the sequence will be retrieved
- name: The feature name. This will appear on the FASTA definition line
- parent_id: (optional) only retrieve a sequence if 'derive_from_parent' is true and the parent matches this ID.
- featureloc_id: (optional) only retrieve a sequence if 'derive_from_parent' is true and the alignment is defined with this featureloc_id
$options: An associative array of options. Valid keys include:
- width: Indicate the number of bases to use per line. A new line will be added after the specified number of bases on each line.
- is_html: Set to '1' if the sequence is meant to be displayed on a web page. This will cause a <br> tag to separate lines of the FASTA sequence.
- derive_from_parent: Set to '1' if the sequence should be obtained from the parent to which this feature is aligned.
- aggregate: Set to '1' if the sequence should only contain sub features, excluding intro sub feature sequence. For example, set this option to obtain just the coding sequence of an mRNA.
- upstream: An integer specifing the number of upstream bases to include in the output
- downstream: An integer specifying the number of downstream bases to include in the output.
- sub_feature_types: Only include sub features (or child features) of the types provided in the array
- relationship_type: If a relationship name is provided (e.g. sequence_of) then any sequences that are in relationships of this type with matched sequences are also included
- relationship_part: If a relationship is provided in the preceeding argument then the rel_part must be either 'object' or 'subject' to indicate which side of the relationship the matched features belong
Return value
an array of matching sequence in the following keys for each sequence: 'types' => an array of feature types that were used to derive the sequence (e.g. from an aggregated sequence) 'upstream' => the number of upstream bases included in the sequence 'downstream' => the number of downstream bases included in the sequence 'defline' => the definintion line used to create a FASTA sequence 'residues' => the residues 'featureloc_id' => the featureloc_id if the sequences is from an alignment
Related topics
2 calls to tripal_get_feature_sequences()
- tripal_feature_sequence.tpl.php in tripal_feature/
theme/ templates/ tripal_feature_sequence.tpl.php - tripal_get_bulk_feature_sequences in tripal_feature/
api/ tripal_feature.api.inc
File
- tripal_feature/
api/ tripal_feature.api.inc, line 100 - Provides an application programming interface (API) for working with features
Code
function tripal_get_feature_sequences($feature, $options) {
// Default values for finding the feature.
$feature_id = array_key_exists('feature_id', $feature) ? $feature['feature_id'] : 0;
$parent_id = array_key_exists('parent_id', $feature) ? $feature['parent_id'] : 0;
$featureloc_id = array_key_exists('featureloc_id', $feature) ? $feature['featureloc_id'] : 0;
$feature_name = array_key_exists('name', $feature) ? $feature['name'] : '';
// Default values for building the sequence.
$num_bases_per_line = array_key_exists('width', $options) ? $options['width'] : 50;
$derive_from_parent = array_key_exists('derive_from_parent', $options) ? $options['derive_from_parent'] : 0;
$aggregate = array_key_exists('aggregate', $options) ? $options['aggregate'] : 0;
$upstream = array_key_exists('upstream', $options) ? $options['upstream'] : 0;
$downstream = array_key_exists('downstream', $options) ? $options['downstream'] : 0;
$sub_features = array_key_exists('sub_feature_types', $options) ? $options['sub_feature_types'] : array();
$relationship = array_key_exists('relationship_type', $options) ? $options['relationship_type'] : '';
$rel_part = array_key_exists('relationship_part', $options) ? $options['relationship_part'] : '';
$is_html = array_key_exists('is_html', $options) ? $options['is_html'] : 0;
if (!$upstream) {
$upstream = 0;
}
if (!$downstream) {
$downstream = 0;
}
// Make sure the sub_features variable is an array.
if (!is_array($sub_features)) {
tripal_report_error('tripal_feature', TRIPAL_ERROR,
"'sub_features' option must be an array for function tripal_get_feature_sequences().",
array()
);
return array();
}
// If a relationship was specified then retreive and the sequences that
// have the given relationship and the recurse to extract the appropriate
// sequence.
if ($rel_part == "object" or $rel_part == "subject") {
if ($rel_part == "subject") {
$sql = '
SELECT FO.feature_id, FO.name, FO.uniquename, CVTO.name as feature_type, O.genus, O.species
FROM {feature} FS
INNER JOIN {feature_relationship} FR ON FR.subject_id = FS.feature_id
INNER JOIN {cvterm} CVTFR ON CVTFR.cvterm_id = FR.type_id
INNER JOIN {feature} FO ON FO.feature_id = FR.object_id
INNER JOIN {cvterm} CVTO ON CVTO.cvterm_id = FO.type_id
INNER JOIN {organism} O ON O.organism_id = FO.organism_id
WHERE
FS.feature_id = :feature_id AND
CVTFR.name = :relationship
';
$features = chado_query($sql, array(':feature_id' => $feature_id, ':relationship' => $relationship));
}
if ($rel_part == "object") {
$sql = '
SELECT FS.feature_id, FS.name, FS.uniquename, CVTO.name as feature_type, O.genus, O.species
FROM {feature} FO
INNER JOIN {feature_relationship} FR ON FR.object_id = FO.feature_id
INNER JOIN {cvterm} CVTFR ON CVTFR.cvterm_id = FR.type_id
INNER JOIN {feature} FS ON FS.feature_id = FR.subject_id
INNER JOIN {cvterm} CVTO ON CVTO.cvterm_id = FS.type_id
INNER JOIN {organism} O ON O.organism_id = FS.organism_id
WHERE
FO.feature_id = :feature_id AND
CVTFR.name = :relationship
';
$features = chado_query($sql, array(':feature_id' => $feature_id, ':relationship' => $relationship));
}
$sequences = '';
while ($feature = $features->fetchObject()) {
// Recurse and get the sequences for these in the relationship.
if ($rel_part == "subject") {
$defline = "$feature_name, $relationship, $feature->uniquename $feature->feature_type ($feature->genus $feature->species)";
}
if ($rel_part == "object") {
$defline = "$feature->uniquename $feature->feature_type ($feature->genus $feature->species), $relationship, $feature_name";
}
return tripal_get_feature_sequences(
array(
'feature_id' => $feature->feature_id,
'name' => $defline,
'parent_id' => $parent_id,
),
array(
'width' => $num_bases_per_line,
'derive_from_parent' => $derive_from_parent,
'aggregate' => $aggregate,
'upstream' => $upstream,
'downstream' => $downstream,
'sub_features' => $sub_features,
)
);
}
}
// Prepare the queries we're going to use later during the render phase
// This SQL statement uses conditionals in the select clause to handle
// cases cases where the alignment is in the reverse direction and when
// the upstream and downstream extensions go beyond the lenght of the
// parent sequence.
$parent_sql = '
SELECT featureloc_id, srcname, srcfeature_id, strand, srctypename, typename,
fmin, fmax, upstream, downstream, adjfmin, adjfmax,
substring(residues from (cast(adjfmin as int4) + 1) for cast((upstream + (fmax - fmin) + downstream) as int4)) as residues,
genus, species
FROM (
SELECT
FL.featureloc_id, OF.name srcname, FL.srcfeature_id, FL.strand,
OCVT.name as srctypename, SCVT.name as typename,
FL.fmin, FL.fmax, OO.genus, OO.species,
CASE
WHEN FL.strand >= 0 THEN
CASE
WHEN FL.fmin - :upstream <= 0 THEN 0
ELSE FL.fmin - :upstream
END
WHEN FL.strand < 0 THEN
CASE
WHEN FL.fmin - :downstream <= 0 THEN 0
ELSE FL.fmin - :downstream
END
END as adjfmin,
CASE
WHEN FL.strand >= 0 THEN
CASE
WHEN FL.fmax + :downstream > OF.seqlen THEN OF.seqlen
ELSE FL.fmax + :downstream
END
WHEN FL.strand < 0 THEN
CASE
WHEN FL.fmax + :upstream > OF.seqlen THEN OF.seqlen
ELSE FL.fmax + :upstream
END
END as adjfmax,
CASE
WHEN FL.strand >= 0 THEN
CASE
WHEN FL.fmin - :upstream <= 0 THEN FL.fmin
ELSE :upstream
END
ELSE
CASE
WHEN FL.fmax + :upstream > OF.seqlen THEN OF.seqlen - FL.fmax
ELSE :upstream
END
END as upstream,
CASE
WHEN FL.strand >= 0 THEN
CASE
WHEN FL.fmax + :downstream > OF.seqlen THEN OF.seqlen - FL.fmax
ELSE :downstream
END
ELSE
CASE
WHEN FL.fmin - :downstream <= 0 THEN FL.fmin
ELSE :downstream
END
END as downstream,
OF.residues
FROM {featureloc} FL
INNER JOIN {feature} SF on FL.feature_id = SF.feature_id
INNER JOIN {cvterm} SCVT on SF.type_id = SCVT.cvterm_id
INNER JOIN {feature} OF on FL.srcfeature_id = OF.feature_id
INNER JOIN {cvterm} OCVT on OF.type_id = OCVT.cvterm_id
INNER JOIN {organism} OO on OF.organism_id = OO.organism_id
WHERE SF.feature_id = :feature_id and NOT (OF.residues = \'\' or OF.residues IS NULL)) as tbl1
';
// This query is meant to get all of the sub features of any given
// feature (arg #1) and order them as they appear on the reference
// feature (arg #2).
$sfsql = '
SELECT SF.feature_id, CVT.name as type_name, SF.type_id
FROM {feature_relationship} FR
INNER JOIN {feature} SF ON SF.feature_id = FR.subject_id
INNER JOIN {cvterm} CVT ON CVT.cvterm_id = SF.type_id
INNER JOIN {featureloc} FL ON FL.feature_id = FR.subject_id
INNER JOIN {feature} PF ON PF.feature_id = FL.srcfeature_id
WHERE FR.object_id = :feature_id and PF.feature_id = :srcfeature_id
ORDER BY FL.fmin ASC
';
// For counting the number of children.
$fsql = '
SELECT count(*) as num_children
FROM {feature_relationship} FR
INNER JOIN {feature} SF ON SF.feature_id = FR.subject_id
INNER JOIN {cvterm} CVT ON CVT.cvterm_id = SF.type_id
INNER JOIN {featureloc} FL ON FL.feature_id = FR.subject_id
INNER JOIN {feature} PF ON PF.feature_id = FL.srcfeature_id
WHERE FR.object_id = :feature_id and PF.feature_id = :srcfeature_id
';
// The array to be returned.
$sequences = array();
// If we need to get the sequence from the parent then do so now.
if ($derive_from_parent) {
// Execute the query to get the sequence from the parent.
$parents = chado_query($parent_sql, array(':upstream' => $upstream, ':downstream' => $downstream, ':feature_id' => $feature_id));
while ($parent = $parents->fetchObject()) {
// If the user specified a particular parent and this one doesn't
// match then skip it.
if ($parent_id and $parent_id != $parent->srcfeature_id) {
continue;
}
// if the user specified a particular featureloc_id and this one
// doesn't match then skip it.
if ($featureloc_id and $featureloc_id != $parent->featureloc_id) {
continue;
}
// Initialize the sequence for each parent.
$seq = '';
$notes = '';
$types = array();
// if we are to aggregate then we will ignore the feature returned
// by the query above and rebuild it using the sub features
if ($aggregate) {
// now get the sub features that are located on the parent.
$children = chado_query($sfsql, array(':feature_id' => $feature_id, ':srcfeature_id' => $parent->srcfeature_id));
$num_children = chado_query($fsql, array(':feature_id' => $feature_id, ':srcfeature_id' => $parent->srcfeature_id))->fetchField();
// Iterate through the sub features and concat their sequences. They
// should already be in order.
$i = 0;
while ($child = $children->fetchObject()) {
// If the callee has specified that only certain sub features should be
// included then continue if this child is not one of those allowed
// subfeatures.
if (count($sub_features) > 0 and !in_array($child->type_name, $sub_features)) {
$i++;
continue;
}
// keep up with the types
if (!in_array($child->type_name, $types)) {
$types[] = $child->type_name;
}
// if the first sub feature we need to include the upstream bases. first check if
// the feature is in the foward direction or the reverse.
if ($i == 0 and $parent->strand >= 0) { // forward direction
// -------------------------- ref
// ....----> ---->
// up 1 2
$q = chado_query($parent_sql, array(':upstream' => $upstream, ':downstream' => 0, ':feature_id' => $child->feature_id));
}
elseif ($i == 0 and $parent->strand < 0) { // reverse direction
// -------------------------- ref
// ....<---- <----
// down 1 2
$q = chado_query($parent_sql, array(':upstream' => 0, ':downstream' => $downstream, ':feature_id' => $child->feature_id));
}
// Next, if the last sub feature we need to include the downstream bases. first check if
// the feature is in teh forward direction or the reverse
elseif ($i == $num_children - 1 and $parent->strand >= 0) { // forward direction
// -------------------------- ref
// ----> ---->....
// 1 2 down
$q = chado_query($parent_sql, array(':upstream' => 0, ':downstream' => $downstream, ':feature_id' => $child->feature_id));
}
elseif ($i == $num_children - 1 and $parent->strand < 0) { // reverse direction
// -------------------------- ref
// <---- <----....
// 1 2 up
$q = chado_query($parent_sql, array(':upstream' => $upstream, ':downstream' => 0, ':feature_id' => $child->feature_id));
}
// for internal sub features we don't want upstream or downstream bases
else {
$q = chado_query($parent_sql, array(':upstream' => 0, ':downstream' => 0, ':feature_id' => $child->feature_id));
}
while ($subseq = $q->fetchObject()) {
// concatenate the sequences of all the sub features
if ($subseq->srcfeature_id == $parent->srcfeature_id) {
$seq .= $subseq->residues;
}
if ($subseq->upstream > 0) {
$notes .= "Includes " . $subseq->upstream . " bases upstream. ";
}
if ($subseq->downstream > 0) {
$notes .= "Includes " . $subseq->downstream . " bases downstream. ";
}
}
$i++;
}
}
// if this isn't an aggregate then use the parent residues
else {
$seq = $parent->residues;
if ($parent->upstream > 0) {
$notes .= "Includes " . $parent->upstream . " bases upstream. ";
}
if ($parent->downstream > 0) {
$notes .= "Includes " . $parent->downstream . " bases downstream. ";
}
}
// get the reverse compliment if feature is on the reverse strand
$dir = 'forward';
$length = strlen($seq);
if ($parent->strand < 0) {
$seq = tripal_reverse_compliment_sequence($seq);
$dir = 'reverse';
}
// now format for display
if ($is_html) {
$seq = wordwrap($seq, $num_bases_per_line, "<br>", TRUE);
}
else {
$seq = wordwrap($seq, $num_bases_per_line, "\r\n", TRUE);
}
if (!$seq) {
$notes .= "No sequence available.";
}
if (count($types) > 0) {
$notes .= "Excludes all bases but those of type(s): " . implode(', ', $types) . ". ";
}
// Construct the definition line for this feature. To construct the
// defline we need a featureloc record, so we'll create one using
// the information we have.
$featureloc = new stdClass;
$featureloc->feature_id = $feature;
$featureloc->fmin = $parent->adjfmin;
$featureloc->fmax = $parent->adjfmax;
$featureloc->strand = $parent->strand;
$featureloc->srcfeature_id = new stdClass;
$featureloc->srcfeature_id->name = $parent->srcname;
$featureloc->srcfeature_id->type_id = $parent->srctypename;
$featureloc->srcfeature_id->organism_id = new stdClass;
$featureloc->srcfeature_id->organism_id->genus = $parent->genus;
$featureloc->srcfeature_id->organism_id->species = $parent->species;
// Get a proper feature object.
$f = chado_generate_var('feature', array('feature_id' => $feature_id));
$defline = tripal_get_fasta_defline($f, $notes, $featureloc, '', $length);
$sequences[] = array(
'types' => $types,
'upstream' => $parent->upstream,
'downstream' => $parent->downstream,
'defline' => $defline,
'residues' => $seq,
'featureloc_id' => $parent->featureloc_id,
'length' => $length,
);
}
}
// If we are not getting the sequence from the parent sequence then
// use what comes through from the feature record.
else {
$f = chado_generate_var('feature', array('feature_id' => $feature_id));
$f = chado_expand_var($f, 'field', 'feature.residues');
$residues = $f->residues;
$length = strlen($residues);
if ($is_html) {
$residues = wordwrap($residues, $num_bases_per_line, "<br>", TRUE);
}
else {
$residues = wordwrap($residues, $num_bases_per_line, "\r\n", TRUE);
}
// get the definintion line for this feature
$defline = tripal_get_fasta_defline($f, '', NULL, '', $length);
// add to the sequence array
$sequences[] = array(
'types' => $f->type_id->name,
'upstream' => 0,
'downstream' => 0,
'defline' => $defline,
'residues' => $residues,
'length' => $length,
);
}
return $sequences;
}