Capstone - Sequence Toolkit

The finished program is a sequence analysis toolkit built on the standard library alone. It parses FASTA and FASTQ records by hand, decodes Phred qualities, trims each read from both ends with a sliding window, measures composition over the bases that are actually bases, searches all six reading frames for open reading frames, translates the best one through a codon table it builds itself, and writes a per-record report.

Every step here is one that Biopython performs for you in a single call. Writing them out is the exercise: the library's decisions about ambiguous bases, about where a reading frame starts, and about what counts as an ORF are decisions, and you cannot audit a result that came out of a call you did not understand.

The data

FASTA. The sandbox mounts sample.fasta:

>NM_000546.6 TP53 transcript variant 1 (excerpt)
ATGGAGGAGCCGCAGTCAGATCCTAGCGTCGAGCCCCCTCTGAGTCAGGAAACATTTTCA
GACCTATGGAAACTACTTCCTGAAAACAACGTTCTGTCCCCCTTGCCGTCCCAAGCAATG
>NM_007294.4 BRCA1 transcript variant 1 (excerpt)
...

A header line begins with >; the identifier is the text after it up to the first whitespace, and the rest is a description you do not need. The sequence follows across one or more lines and must be concatenated, as line breaks represent file formatting conventions rather than biological features.

FASTQ. The sandbox mounts reads.fastq:

@read001 length=60
ATGGAGGAGCCGCAGTCAGATCCTAGCGTCGAGCCCCCTCTGAGTCAGGAAACATTTTCA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII99755
@read002 length=60
GACCTATGGAAACTACTTCCTGAAAACAACGTTCTGTCCCCCTTGCCGTCCCAAGCAATG
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIHHHHHGGGGFFFEEDDCCBA@?>=<

Four lines per record: a header beginning @, the sequence, a separator beginning +, and one quality character per base. The qualities are Phred+33, so a character's score is its code point minus 33 (e.g. I corresponds to 40, 9 to 24, 5 to 20, and ! to 0). The mounted file is well formed: all three records carry one quality character per base. Your program must not assume that of its input. A sequence and quality line of different lengths is the most common way a truncated download announces itself, which is why requirement 2 compares the two lengths before downstream processing; zipping unequal lines would erroneously truncate quality score associations for trailing bases. The falling qualities at the end of read002 are the other thing to notice, and they are what the trimming step exists for.

Sequences may hold characters other than A, C, G and T; N is the usual one, and the toolkit has to decide what it means at every stage rather than once.

Requirements

1. Read the trimming window W, the quality threshold Q and the minimum ORF length M in amino acids as parameters, and read the FASTQ records four lines at a time.
2. Treat a record as malformed, count it, and drop it when the header does not begin with @ or carries no identifier after it, when the separator line does not begin with +, or when the sequence and quality lines differ in length. A trailing group of fewer than four lines is one malformed record.
3. Decode the qualities as Phred+33. Trim from the 5' end one base at a time while the read is at least W long and the total of its first W scores is below Q * W, then trim the same way from the 3' end. Compare totals rather than means so the decision stays in integer arithmetic and does not turn on a rounding.
4. Count a read left shorter than W as having failed quality control, and drop it.
5. Report the trimmed read's length, its GC fraction over the unambiguous bases alone (or 0.0 when none are present), alongside the count of bases deviating from A, C, G, or T. An ambiguous base is in the length and out of the GC denominator; the two counts answer different questions.
6. Search all six reading frames: the three offsets of the read and the three of its reverse complement, where the complement of an ambiguous base is N. Within a frame an ORF runs from an ATG codon to the next in-frame stop codon TAA, TAG or TGA, and the codon after that stop may begin another one.
7. Keep only ORFs of at least M amino acids, counting the codons before the stop, and report the single longest ORF of the record. Break a tie by frame in the order +1, +2, +3, -1, -2, -3, and then by the earlier start.
8. Build a codon table for the standard genetic code and translate the chosen ORF with it, giving X to any codon containing an ambiguous base and stopping before the terminator.
9. Emit the record, malformed and failed counts, one line per surviving record in input order (using ., 0, and - placeholders for records lacking qualifying ORFs), followed by a summary line detailing surviving record count, ORF-positive count, and mean trimmed length.

Milestones

Run these in order. They give you the two parsers and the codon table; the trimming policy, the ORF search and the report are yours.

Parse FASTA by hand, joining the wrapped lines back together.

NM_000546.6 120 0.525
NM_007294.4 120 0.408
ENST00000646891.2 120 0.792

The final yield after the loop is the part that is easy to leave out, and leaving it out silently drops the last record of every file you ever read.

Decode the qualities and see what the sliding window does to each read.

read001 60 -> 56 min Q 20 max Q 40
read002 60 -> 52 min Q 27 max Q 40
read003 60 -> 60 min Q 40 max Q 40

Every record passes the length test. The first gives up four bases and the second eight, both from falling tails, while the third never drops below I and keeps all sixty. Delete one quality character from the file and the record it belongs to should be rejected immediately, ensuring that corrupted records are not improperly evaluated against misaligned quality strings.

Lower the threshold to 30 and read002 gives up a single base instead of eight. The window and the threshold are not incidental settings; they are the policy, and the report is only as reproducible as your record of what they were.

Build the codon table and translate the three forward frames of one sequence.

table size: 64
ATG -> M | TAA -> * | GGN -> X
reverse complement: GGGTTACCCGGGTTTCAT
+1 MKPG*P
+2 *NPGN
+3 ETRVT

Iterating TCAG three deep produces the 64 codons in exactly the order of the NCBI standard table, so the amino acid string above lines up with it position for position. That is the whole table in three lines, and it is also why the order of those four letters is not arbitrary and must not be sorted.

What to watch for

• Reverse complementing without reversing, or reversing without complementing. Both produce a sequence of the right length in the wrong strand, and the ORF search will happily report a peptide from it.
• Taking the reverse complement of the untrimmed read. The three reverse frames are offsets into the reverse complement of what survived trimming; computing them from the original shifts every one of them.
• Trimming from the 3' end first. The two passes are not commutative, because each one changes the length the other's window test is bounded by.
• Letting the start position survive a stop codon. An ATG inside an open reading frame is a methionine, not a second start, but an ATG after the stop begins a new ORF; tracking must clear active start coordinates upon reaching valid stop codons.
• Counting the stop codon in the ORF length, or translating it into a * in the peptide. The requirement is the codons before it.
• Putting ambiguous bases in the GC denominator. A read that is half N then reports a GC fraction near 0.25 whatever its real composition.
• Assuming four lines per record without checking. A truncated file ends mid-record, and slicing four lines out of a list that has three left is not an error in Python, but returns a truncated list of three lines.

Exercise

Write the whole toolkit in one function, with the standard library alone. Biopython is not installed here, so the FASTQ parsing, the Phred decoding, the reverse complement, the codon table and the translation are all yours; that is the point of the exercise, because it is exactly the work the library was doing.

Write toolkit(text), returning a list of output lines. The first input line holds three whitespace-separated integers: the trimming window W, the quality threshold Q, and the minimum ORF length M in amino acids. Ignore blank lines in the rest, then read what remains as FASTQ records of four lines: header, sequence, separator, qualities. A record is malformed, and is counted and dropped, when its header does not begin with @ or carries no identifier after it, when its separator line does not begin with +, or when its sequence and quality lines differ in length. A trailing group of fewer than four lines counts as one malformed record. The identifier is the header text after the @ up to the first whitespace.

Decode the qualities as Phred+33, so a character's score is its code point minus 33. Trim the read from the 5' end while it is at least W long and the total of its first W scores is below Q * W, dropping one base at a time, then trim it the same way from the 3' end. A read left shorter than W fails quality control: count it and drop it. Uppercase the surviving sequence.

Report the composition of the trimmed read: its length, its GC fraction over the unambiguous bases alone -- 0.0 when there are none -- and the number of bases that are not A, C, G or T. Then search all six reading frames, the three offsets of the read and the three of its reverse complement, taking the complement of an ambiguous base to be N. Within a frame an ORF runs from an ATG codon to the next in-frame stop codon TAA, TAG or TGA, and the codon after that stop may start another. Keep only ORFs of at least M amino acids, counting the codons before the stop, and report the single longest ORF of the record, breaking a tie by frame in the order +1, +2, +3, -1, -2, -3 and then by the earlier start. Translate it with a codon table you build from the standard genetic code, giving X to any codon containing an ambiguous base and stopping before the terminator.

Print records=<n>, then malformed=<n>, then failed_qc=<n>. Follow those with one line per surviving record, in input order, holding the identifier, the trimmed length, the GC fraction to exactly three decimal places, the ambiguous base count, the frame, the ORF length in amino acids and the peptide, separated by tabs; a record with no qualifying ORF takes . for the frame, 0 for the length and - for the peptide. If no record survives, print no reads passed in place of the table. Close with summary: <n> reads, <k> with an ORF, mean trimmed length <x>, where n is the number of surviving records, k how many carried an ORF, and x their mean trimmed length to exactly one decimal place, or 0.0 when none survived.