Fastq

5 problems from Rosalind — Bioinformatics Armory. Press Run on any block to execute it in your browser.

# Rosalind: TFSQ — FASTQ Format Introduction
# https://rosalind.info/problems/tfsq/
#
# Given: FASTQ entries.
# Return: Corresponding FASTA records (strip quality, change @ to >).

# Sample FASTQ data
let fastq_id = "SEQ_ID"
let fastq_seq = "GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT"

# Convert to FASTA format
let fasta_text = ">" + fastq_id + "\n" + fastq_seq

println("Result:")
println(fasta_text)
println("")
println("Expected:")
println(">SEQ_ID")
println("GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT")

let ok = fasta_text == ">SEQ_ID\nGATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT"
println("\nMatch: " + str(ok))

fn test_tfsq_fastq_to_fasta() {
    assert ok, "TFSQ: FASTA conversion did not match the expected record"
}

# Rosalind: PHRE — Read Quality Distribution
# https://rosalind.info/problems/phre/
#
# Given: A quality threshold and FASTQ entries.
# Return: Number of reads whose average quality is below the threshold.

let threshold = 28

# Sample reads with Phred33 quality strings
let quality_strings = [
    "6.3536354;.151<211/0?::6/-2051)-*\"40/.,+%)",
    "AH@FGGGJ<GB<<9:GD=D@GG9=?A@DC=;:?>839/4856",
    "@DJEJEA?JHJ@8?F?IA3=;8@C95=;=?;>D/:;74792."
]

# Parse Phred33: each character's ASCII value minus 33 gives quality score
let below_count = 0
for qual_str in quality_strings {
    let chars = split(qual_str, "")
    let scores = chars |> filter(|c| c != "") |> map(|c| ascii(c) - 33)
    let avg = (scores |> reduce(|a, b| a + b)) / len(scores)
    if avg < threshold then
        below_count = below_count + 1
}

println("Result:   " + str(below_count))
println("Expected: 1")
println("Match:    " + str(below_count == 1))

fn test_phre_reads_below_threshold() {
    assert below_count == 1, "PHRE: expected 1, got " + str(below_count)
}

# Rosalind: FILT — Read Filtration by Quality
# https://rosalind.info/problems/filt/
#
# Given: Quality threshold q, percentage p, and FASTQ entries.
# Return: Number of reads where at least p% of bases have quality >= q.

let q = 20
let p = 90

let reads = [
    {id: "Rosalind_0049_1", seq: "GCAGAGACCAGTAGATGTGTTTGCGGACGGTCGGGCTCCATGTGACACAG",
     qual: "FD@@;C<AI?4BA:=>C<G=:AE=><A??>764A8B797@A:58:527+,"},
    {id: "Rosalind_0049_2", seq: "AATGGGGGGGGGAGACAAAATACGGCTAAGGCAGGGGTCCTTGATGTCAT",
     qual: "1<<65:793967<4:92568-34:.>1;2752)24')*15;1,.3*3+*!"},
    {id: "Rosalind_0049_3", seq: "ACCCCATACGGCGAGCGTCAGCATCTGATATCCTCTTTCAATCCTAGCTA",
     qual: "B:EI>JDB5=>DA?E6B@@CA?C;=;@@C:6D:3=@49;@87;::;;?8+"}
]

let passing = 0
for read in reads {
    let chars = split(read.qual, "") |> filter(|c| c != "")
    let scores = chars |> map(|c| ascii(c) - 33)
    let good = scores |> filter(|s| s >= q) |> count()
    let pct = (good * 100) / len(scores)
    if pct >= p then
        passing = passing + 1
}

println("Result:   " + str(passing))
println("Expected: 2")
println("Match:    " + str(passing == 2))

fn test_filt_reads_passing() {
    assert passing == 2, "FILT: expected 2, got " + str(passing)
}

# Rosalind: BPHR — Base Quality Distribution
# https://rosalind.info/problems/bphr/
#
# Given: FASTQ file and quality threshold q.
# Return: Number of positions where mean base quality falls below q.

let threshold = 26

let quality_strings = [
    ">?F?@6<C<HF?<85486B;85:8488/2/",
    "@J@H@>B9:B;<D==:<;:,<::?463-,,",
    "=88;99637@5,4664-65)/?4-2+)$)$",
    "<@BGE@8C9=B9:B<>>>7?B>7:02+33."
]

# All reads should be the same length
let read_len = len(split(quality_strings[0], "") |> filter(|c| c != ""))

# Calculate mean quality at each position across all reads
let below_count = 0
let i = 0
while i < read_len {
    let total = 0
    for qual_str in quality_strings {
        let chars = split(qual_str, "") |> filter(|c| c != "")
        total = total + ascii(chars[i]) - 33
    }
    let mean = total / len(quality_strings)
    if mean < threshold then
        below_count = below_count + 1
    i = i + 1
}

println("Result:   " + str(below_count))
println("Expected: 17")
println("Match:    " + str(below_count == 17))

fn test_bphr_positions_below_threshold() {
    assert below_count == 17, "BPHR: expected 17, got " + str(below_count)
}

# Rosalind: BFIL — Base Filtration by Quality
# https://rosalind.info/problems/bfil/
#
# Given: FASTQ file and quality threshold q (Phred33).
# Return: FASTQ with leading and trailing low-quality bases trimmed.

let q = 20

let reads = [
    {id: "Rosalind_0049",
     seq:  "GCAGAGACCAGTAGATGTGTTTGCGGACGGTCGGGCTCCATGTGACACAG",
     qual: "FD@@;C<AI?4BA:=>C<G=:AE=><A??>764A8B797@A:58:527+,"},
    {id: "Rosalind_0049",
     seq:  "AATGGGGGGGGGAGACAAAATACGGCTAAGGCAGGGGTCCTTGATGTCAT",
     qual: "1<<65:793967<4:92568-34:.>1;2752)24')*15;1,.3*3+*!"},
    {id: "Rosalind_0049",
     seq:  "ACCCCATACGGCGAGCGTCAGCATCTGATATCCTCTTTCAATCCTAGCTA",
     qual: "B:EI>JDB5=>DA?E6B@@CA?C;=;@@C:6D:3=@49;@87;::;;?8+"}
]

let trimmed_reads = []

println("Result:")
for read in reads {
    let chars_seq = split(read.seq, "") |> filter(|c| c != "")
    let chars_qual = split(read.qual, "") |> filter(|c| c != "")
    let scores = chars_qual |> map(|c| ascii(c) - 33)

    # Find first position from left with quality >= q
    let start = 0
    while start < len(scores) && scores[start] < q {
        start = start + 1
    }
    # Find last position from right with quality >= q
    let end_pos = len(scores) - 1
    while end_pos >= 0 && scores[end_pos] < q {
        end_pos = end_pos - 1
    }

    let trimmed_seq = chars_seq |> slice(start, end_pos + 1) |> join("")
    let trimmed_qual = chars_qual |> slice(start, end_pos + 1) |> join("")
    trimmed_reads = push(trimmed_reads, trimmed_seq)

    println("@" + read.id)
    println(trimmed_seq)
    println("+")
    println(trimmed_qual)
}

println("\nExpected:")
println("@Rosalind_0049")
println("GCAGAGACCAGTAGATGTGTTTGCGGACGGTCGGGCTCCATGTGACAC")
println("+")
println("FD@@;C<AI?4BA:=>C<G=:AE=><A??>764A8B797@A:58:527")
println("@Rosalind_0049")
println("ATGGGGGGGGGAGACAAAATACGGCTAAGGCAGGGGTCCT")
println("+")
println("<<65:793967<4:92568-34:.>1;2752)24')*15;")
println("@Rosalind_0049")
println("ACCCCATACGGCGAGCGTCAGCATCTGATATCCTCTTTCAATCCTAGCT")
println("+")
println("B:EI>JDB5=>DA?E6B@@CA?C;=;@@C:6D:3=@49;@87;::;;?8")

fn test_bfil_trims_low_quality_ends() {
    let expected = [
        "GCAGAGACCAGTAGATGTGTTTGCGGACGGTCGGGCTCCATGTGACAC",
        "ATGGGGGGGGGAGACAAAATACGGCTAAGGCAGGGGTCCT",
        "ACCCCATACGGCGAGCGTCAGCATCTGATATCCTCTTTCAATCCTAGCT"
    ]
    assert len(trimmed_reads) == 3, "BFIL: expected 3 reads, got " + str(len(trimmed_reads))
    let i = 0
    while i < 3 {
        assert trimmed_reads[i] == expected[i], "BFIL: read " + str(i + 1) + " trimmed to " + trimmed_reads[i]
        i = i + 1
    }
}