///|
fn sum_feature_lengths(features : Array[Feature]) -> Int {
let mut total = 0
for feature in features {
total += feature.length()
}
total
}
///|
fn intron_bases_from_exons(exons : Array[Feature]) -> Int {
if exons.length() < 2 {
0
} else {
let mut total = 0
for i in 0..<(exons.length() - 1) {
let gap = exons[i + 1].start - exons[i].end - 1
if gap > 0 {
total += gap
}
}
total
}
}
///|
pub fn TranscriptModel::metrics(self : TranscriptModel) -> TranscriptMetrics {
let exon_bases = sum_feature_lengths(self.exons)
let cds_bases = sum_feature_lengths(self.cds)
let utr_bases = sum_feature_lengths(self.utrs)
{
transcript_id: self.id,
exon_count: self.exons.length(),
cds_count: self.cds.length(),
utr_count: self.utrs.length(),
exon_bases,
cds_bases,
utr_bases,
intron_count: if self.exons.length() > 0 {
self.exons.length() - 1
} else {
0
},
intron_bases: intron_bases_from_exons(self.exons),
span_bases: self.feature.length(),
}
}
///|
fn choose_longest_transcript(
transcripts : Array[TranscriptModel],
) -> (String, Int) {
if transcripts.length() == 0 {
("", 0)
} else {
let first_metrics = transcripts[0].metrics()
let mut best_id = first_metrics.transcript_id
let mut best_bases = first_metrics.exon_bases
for tx in transcripts {
let metrics = tx.metrics()
if metrics.exon_bases > best_bases {
best_id = metrics.transcript_id
best_bases = metrics.exon_bases
}
}
(best_id, best_bases)
}
}
///|
pub fn GeneModel::metrics(self : GeneModel) -> GeneMetrics {
let mut total_exon_bases = 0
let mut total_cds_bases = 0
for tx in self.transcripts {
let metrics = tx.metrics()
total_exon_bases += metrics.exon_bases
total_cds_bases += metrics.cds_bases
}
let (longest_transcript_id, longest_transcript_bases) = choose_longest_transcript(
self.transcripts,
)
{
gene_id: self.id,
transcript_count: self.transcripts.length(),
total_exon_bases,
total_cds_bases,
longest_transcript_id,
longest_transcript_bases,
gene_span_bases: self.feature.length(),
}
}
///|
pub fn Annotation::metrics(self : Annotation) -> AnnotationMetrics {
{
feature_count: self.features.length(),
gene_count: self.filter_by_type("gene").length(),
transcript_count: self.filter_by_type("transcript").length() +
self.filter_by_type("mRNA").length(),
exon_count: self.filter_by_type("exon").length(),
cds_count: self.filter_by_type("CDS").length(),
seqid_count: self.build_interval_index().seqids().length(),
total_feature_bases: sum_feature_lengths(self.features),
}
}
///|
pub fn Annotation::type_report(self : Annotation) -> Array[TypeReportRow] {
let rows : Array[TypeReportRow] = []
for count in self.count_by_type() {
let features = self.filter_by_type(count.feature_type)
rows.push({
feature_type: count.feature_type,
count: count.count,
total_bases: sum_feature_lengths(features),
})
}
rows
}
///|
fn has_attribute_key(keys : Array[String], key : String) -> Bool {
for item in keys {
if item == key {
return true
}
}
false
}
///|
fn count_attribute_key(features : Array[Feature], key : String) -> Int {
let mut count = 0
for feature in features {
for attr in feature.attributes {
if attr.key == key {
count += 1
}
}
}
count
}
///|
pub fn Annotation::attribute_frequencies(
self : Annotation,
) -> Array[AttributeFrequency] {
let keys : Array[String] = []
let result : Array[AttributeFrequency] = []
for feature in self.features {
for attr in feature.attributes {
if !has_attribute_key(keys, attr.key) {
keys.push(attr.key)
result.push({
key: attr.key,
count: count_attribute_key(self.features, attr.key),
})
}
}
}
result
}
///|
pub fn TranscriptMetrics::coding_ratio(self : TranscriptMetrics) -> Double {
if self.exon_bases == 0 {
0.0
} else {
self.cds_bases.to_double() / self.exon_bases.to_double()
}
}
///|
pub fn TranscriptMetrics::has_cds(self : TranscriptMetrics) -> Bool {
self.cds_bases > 0
}
///|
pub fn TranscriptMetrics::has_introns(self : TranscriptMetrics) -> Bool {
self.intron_count > 0
}
///|
pub fn GeneMetrics::is_multi_transcript(self : GeneMetrics) -> Bool {
self.transcript_count > 1
}
///|
pub fn AnnotationMetrics::is_empty(self : AnnotationMetrics) -> Bool {
self.feature_count == 0
}
///|
pub fn AnnotationMetrics::mean_feature_length(
self : AnnotationMetrics,
) -> Double {
if self.feature_count == 0 {
0.0
} else {
self.total_feature_bases.to_double() / self.feature_count.to_double()
}
}