amplicon_workflow.py · Scientific Agent Skills

1#!/usr/bin/env python3

2"""Validate paired-end amplicons and run a provenance-preserving QIIME 2 workflow."""

3from __future__ import annotations

4

5import argparse

6import csv

7import gzip

8import hashlib

9import itertools

10import json

11import math

12from pathlib import Path

13import re

14import subprocess

15

16IUPAC = dict(zip('ACGTRYSWKMBDHVN', ['A', 'C', 'G', 'T', 'AG', 'CT', 'GC', 'AT', 'GT', 'AC', 'CGT', 'AGT', 'ACT', 'ACG', 'ACGT']))

17ID_HEADERS = {'id', 'sampleid', 'sample id', 'sample-id', 'featureid', 'feature id', 'feature-id'}

18LEGACY_ID_HEADERS = {'#SampleID', '#Sample ID', '#OTUID', '#OTU ID', 'sample_name'}

19

20

21def id_header(value: str) -> bool:

22 return value.lower() in ID_HEADERS or value in LEGACY_ID_HEADERS

23

24

25def fastq(path: Path):

26 opener = gzip.open if path.suffix == '.gz' else open

27 with opener(path, 'rt', encoding='ascii') as handle:

28 while True:

29 name = handle.readline()

30 if not name:

31 break

32 sequence, plus, quality = (handle.readline().rstrip('\r\n') for _ in range(3))

33 if not name.startswith('@') or not plus.startswith('+') or not sequence or len(sequence) != len(quality):

34 raise ValueError(f'Malformed four-line FASTQ record in {path}')

35 if set(sequence.upper()) - set('ACGTN') or any(not 33 <= ord(c) <= 126 for c in quality):

36 raise ValueError(f'Invalid sequence or Phred+33 characters in {path}')

37 tokens = name[1:].split()

38 identifier = re.sub(r'/[12]$', '', tokens[0]) if tokens else ''

39 if not identifier:

40 raise ValueError(f'Empty FASTQ read identifier in {path}')

41 yield identifier, sequence.upper()

42

43

44def primer_match(sequence: str, primer: str) -> bool:

45 return len(sequence) >= len(primer) and all(base in IUPAC[p] for base, p in zip(sequence, primer))

46

47

48def read_metadata(path: Path) -> set[str]:

49 with path.open(encoding='utf-8-sig') as handle:

50 rows = [[cell.strip() for cell in row] for row in csv.reader(handle, delimiter='\t')]

51 rows = [row for row in rows if any(row) and

52 (not row[0].startswith('#') or row[0] in LEGACY_ID_HEADERS)]

53 if not rows:

54 raise ValueError('Metadata is empty')

55 if not id_header(rows[0][0]):

56 raise ValueError('Use sample-id as the first metadata header')

57 columns = rows[0][1:]

58 if any(not c or id_header(c) for c in columns) or len({c.lower() for c in columns}) != len(columns):

59 raise ValueError('Metadata has empty, reserved, or duplicate column names')

60 if any(len(row) != len(rows[0]) for row in rows[1:]):

61 raise ValueError('Metadata rows must have the same number of columns as the header')

62 samples = [row[0] for row in rows[1:]]

63 if not samples or any(not s or id_header(s) for s in samples) or len(samples) != len(set(samples)):

64 raise ValueError('Metadata has empty, reserved, or duplicate sample IDs')

65 return set(samples)

66

67

68def validate(manifest: Path, metadata: Path, primer_f: str, primer_r: str,

69 trunc_f: int, trunc_r: int, amplicon_max: int, min_overlap: int = 12) -> dict:

70 primers = (primer_f.upper(), primer_r.upper())

71 if any(not p or set(p) - set(IUPAC) for p in primers):

72 raise ValueError('Both primers must be nonempty IUPAC DNA, written 5-prime to 3-prime as sequenced')

73 if min(trunc_f, trunc_r, amplicon_max) <= 0 or min_overlap < 12:

74 raise ValueError('Positive post-primer truncation/amplicon lengths and overlap >=12 are required')

75 overlap = trunc_f + trunc_r - amplicon_max

76 if overlap < min_overlap:

77 raise ValueError(f'Only {overlap} bp overlap remains; need at least {min_overlap}')

78 with manifest.open(newline='', encoding='utf-8-sig') as handle:

79 reader = csv.DictReader(handle, delimiter='\t')

80 fields = ['sample-id', 'forward-absolute-filepath', 'reverse-absolute-filepath']

81 if reader.fieldnames != fields:

82 raise ValueError(f'PairedEndFastqManifestPhred33V2 requires tab-separated columns {fields}')

83 rows = list(reader)

84 seen, paths, results = set(), set(), []

85 for row in rows:

86 if None in row or any(row[key] is None or not row[key].strip() for key in fields):

87 raise ValueError('Manifest rows must contain exactly three nonempty tab-separated cells')

88 sample = row['sample-id'].strip()

89 if sample.startswith('#') or id_header(sample) or sample in seen:

90 raise ValueError(f'Reserved or duplicate sample ID: {sample}')

91 seen.add(sample)

92 files = [Path(row[key].strip()) for key in fields[1:]]

93 if any(not p.is_absolute() or not p.is_file() for p in files):

94 raise ValueError(f'{sample}: FASTQ paths must be absolute existing files')

95 if len({p.resolve() for p in files}) != 2 or any(p.resolve() in paths for p in files):

96 raise ValueError(f'{sample}: repeated forward/reverse FASTQ file')

97 paths.update(p.resolve() for p in files)

98 count, matched, too_short, long_pairs = 0, [0, 0], [0, 0], 0

99 for forward, reverse in itertools.zip_longest(fastq(files[0]), fastq(files[1])):

100 if forward is None or reverse is None or forward[0] != reverse[0]:

101 raise ValueError(f'{sample}: paired read IDs/order/counts do not match')

102 count += 1

103 long_pairs += all(len(read[1]) - len(primers[i]) >= (trunc_f, trunc_r)[i]

104 for i, read in enumerate((forward, reverse)))

105 for i, read in enumerate((forward, reverse)):

106 sequence = read[1]

107 if count <= 1000:

108 matched[i] += primer_match(sequence, primers[i])

109 too_short[i] += len(sequence) - len(primers[i]) < (trunc_f, trunc_r)[i]

110 if long_pairs == 0:

111 raise ValueError(f'{sample}: no read pairs meet nominal truncation lengths after primer removal')

112 results.append({'sample_id': sample, 'raw_pairs': count,

113 'nominal_length_eligible_pairs': long_pairs,

114 'exact_primer_fraction_first_1000': [v/min(count, 1000) for v in matched],

115 'short_fraction_before_quality_filter': [v/count for v in too_short]})

116 if not seen or read_metadata(metadata) != seen:

117 raise ValueError('Manifest and metadata sample IDs must match exactly')

118 warnings = [f"{s['sample_id']}: fewer than 80% exact primer-prefix matches; inspect orientation/adapters"

119 for s in results if min(s['exact_primer_fraction_first_1000']) < 0.8]

120 return {'samples': results, 'minimum_predicted_overlap': overlap, 'warnings': warnings}

121

122

123def retention(path: Path, raw_counts: dict[str, int] | None = None) -> dict:

124 if raw_counts is not None and any(type(v) is not int or v < 0 for v in raw_counts.values()):

125 raise ValueError('Raw read-pair counts must be nonnegative integers')

126 with path.open() as handle:

127 reader = csv.DictReader((line for line in handle if not line.startswith('#q2:')), delimiter='\t')

128 rows = list(reader)

129 summaries, warnings = [], []

130 for row in rows:

131 sample = row.get('sample-id', row.get('sampleid', row.get('#SampleID')))

132 if not sample or not sample.strip():

133 raise ValueError('DADA2 statistics lack sample IDs')

134 sample = sample.strip()

135 if 'concatenated' in row:

136 raise ValueError('Retention helper requires merged-only DADA2 statistics (retain_unmerged=False)')

137 counts = [float(row[key]) for key in ['input', 'filtered', 'denoised', 'merged', 'non-chimeric']]

138 if any(not math.isfinite(x) or x < 0 or not x.is_integer() for x in counts):

139 raise ValueError('Invalid DADA2 read counts')

140 if any(b > a for a, b in zip(counts, counts[1:])):

141 raise ValueError('DADA2 retained counts increase between stages')

142 if raw_counts is not None and sample not in raw_counts:

143 raise ValueError('DADA2 statistics contain samples absent from raw read counts')

144 raw = raw_counts[sample] if raw_counts is not None else int(counts[0])

145 if counts[0] > raw:

146 raise ValueError('DADA2 input exceeds raw read pairs')

147 fraction = counts[-1]/raw if raw else 0.0

148 basis = 'raw pairs' if raw_counts is not None else 'DADA2 input pairs'

149 if fraction < 0.5:

150 warnings.append(f'{sample}: retained {fraction:.1%} of {basis}; inspect trimming, quality, merging, chimeras')

151 summaries.append({'sample_id': sample, 'raw_pairs': raw if raw_counts is not None else None,

152 'retention_denominator': basis, 'dada2_input': int(counts[0]),

153 'non_chimeric': int(counts[-1]), 'retained_fraction': fraction})

154 ids = [r['sample_id'] for r in summaries]

155 if not ids or len(ids) != len(set(ids)) or (raw_counts is not None and set(ids) != set(raw_counts)):

156 raise ValueError('Missing/duplicate samples in DADA2 statistics')

157 return {'samples': summaries, 'warnings': warnings}

158

159

160def validate_runtime(info: str) -> None:

161 versions = {key.strip().replace('_', '-'): value.strip()

162 for line in info.splitlines() if ':' in line

163 for key, value in [line.split(':', 1)]}

164 required = ['rachis version', 'q2cli version', 'cutadapt', 'dada2', 'demux',

165 'feature-table', 'feature-classifier', 'taxa', 'types']

166 incompatible = [name for name in required

167 if not re.fullmatch(r'2026\.7\.\d+', versions.get(name, ''))]

168 if incompatible:

169 raise ValueError('This runner requires stable QIIME 2 2026.7 and matching plugins; '

170 'missing or incompatible: ' + ', '.join(incompatible))

171

172

173def run(args) -> dict:

174 qc = validate(args.manifest, args.metadata, args.primer_f, args.primer_r,

175 args.trunc_f, args.trunc_r, args.amplicon_max, args.min_overlap)

176 output = args.output.resolve()

177 if output.exists() and any(output.iterdir()):

178 raise ValueError('Use a fresh output directory; retain prior QIIME artifacts and provenance')

179 if not args.classifier.is_file():

180 raise ValueError('Provide a release-compatible TaxonomicClassifier .qza')

181 info = subprocess.run([args.qiime, 'info'], check=True, text=True, capture_output=True).stdout

182 validate_runtime(info)

183 output.mkdir(parents=True, exist_ok=True)

184 (output/'qiime-info.txt').write_text(info)

185 (output/'input-qc.json').write_text(json.dumps(qc, indent=2)+'\n')

186 commands = []

187 def invoke(*parts):

188 command = [args.qiime, *map(str, parts)]

189 commands.append(command)

190 (output/'commands.json').write_text(json.dumps(commands, indent=2)+'\n')

191 with (output/'workflow.log').open('a') as log:

192 subprocess.run(command, stdout=log, stderr=subprocess.STDOUT, check=True)

193 invoke('tools', 'validate', args.classifier.resolve(), '--level', 'max')

194 invoke('tools', 'import', '--type', 'SampleData[PairedEndSequencesWithQuality]',

195 '--input-format', 'PairedEndFastqManifestPhred33V2', '--input-path', args.manifest.resolve(),

196 '--output-path', output/'demux.qza')

197 invoke('cutadapt', 'trim-paired', '--i-demultiplexed-sequences', output/'demux.qza',

198 '--p-front-f', '^'+args.primer_f.upper(), '--p-front-r', '^'+args.primer_r.upper(),

199 '--p-discard-untrimmed', '--p-cores', args.threads, '--output-dir', output/'cutadapt')

200 invoke('demux', 'summarize', '--i-data', output/'cutadapt/trimmed_sequences.qza', '--o-visualization', output/'trimmed.qzv')

201 invoke('dada2', 'denoise-paired', '--i-demultiplexed-seqs', output/'cutadapt/trimmed_sequences.qza',

202 '--p-trunc-len-f', args.trunc_f, '--p-trunc-len-r', args.trunc_r,

203 '--p-trim-left-f', 0, '--p-trim-left-r', 0, '--p-min-overlap', args.min_overlap,

204 '--p-n-threads', args.threads, '--output-dir', output/'denoise')

205 invoke('feature-table', 'summarize', '--i-table', output/'denoise/table.qza',

206 '--m-metadata-file', args.metadata.resolve(), '--o-summary', output/'table.qzv',

207 '--o-feature-frequencies', output/'feature-frequencies.qza',

208 '--o-sample-frequencies', output/'sample-frequencies.qza')

209 invoke('feature-classifier', 'classify-sklearn', '--i-classifier', args.classifier.resolve(),

210 '--i-reads', output/'denoise/representative_sequences.qza', '--p-n-jobs', args.threads,

211 '--o-classification', output/'taxonomy.qza')

212 invoke('taxa', 'barplot', '--i-table', output/'denoise/table.qza', '--i-taxonomy', output/'taxonomy.qza',

213 '--m-metadata-file', args.metadata.resolve(), '--o-visualization', output/'taxa.qzv')

214 invoke('tools', 'export', '--input-path', output/'denoise/denoising_stats.qza',

215 '--output-path', output/'stats')

216 artifacts = sorted(output.rglob('*.qza'))

217 for artifact in artifacts:

218 invoke('tools', 'validate', artifact, '--level', 'max')

219 result = retention(output/'stats/stats.tsv', {s['sample_id']: s['raw_pairs'] for s in qc['samples']})

220 result.update({'classifier_sha256': digest(args.classifier),

221 'manifest_sha256': digest(args.manifest), 'metadata_sha256': digest(args.metadata),

222 'artifacts': {str(p.relative_to(output)): digest(p) for p in artifacts}})

223 (output/'retention-qc.json').write_text(json.dumps(result, indent=2)+'\n')

224 return result

225

226

227def digest(path: Path) -> str:

228 value = hashlib.sha256()

229 with path.open('rb') as handle:

230 for block in iter(lambda: handle.read(1024*1024), b''):

231 value.update(block)

232 return value.hexdigest()

233

234

235def main() -> None:

236 parser = argparse.ArgumentParser(description=__doc__)

237 sub = parser.add_subparsers(dest='action', required=True)

238 for action in ['validate', 'run']:

239 command = sub.add_parser(action)

240 command.add_argument('--manifest', type=Path, required=True)

241 command.add_argument('--metadata', type=Path, required=True)

242 command.add_argument('--primer-f', required=True)

243 command.add_argument('--primer-r', required=True)

244 command.add_argument('--trunc-f', type=int, required=True)

245 command.add_argument('--trunc-r', type=int, required=True)

246 command.add_argument('--amplicon-max', type=int, required=True, help='Maximum expected biological insert length AFTER removing primers')

247 command.add_argument('--min-overlap', type=int, default=12)

248 if action == 'run':

249 command.add_argument('--classifier', type=Path, required=True)

250 command.add_argument('--output', type=Path, required=True)

251 command.add_argument('--qiime', default='qiime')

252 command.add_argument('--threads', type=int, default=1)

253 stats = sub.add_parser('retention')

254 stats.add_argument('stats_tsv', type=Path)

255 args = parser.parse_args()

256 try:

257 if args.action == 'run':

258 if args.threads < 1:

259 raise ValueError('threads must be positive')

260 result = run(args)

261 elif args.action == 'validate':

262 result = validate(args.manifest, args.metadata, args.primer_f, args.primer_r,

263 args.trunc_f, args.trunc_r, args.amplicon_max, args.min_overlap)

264 else:

265 result = retention(args.stats_tsv)

266 except (ValueError, OSError, EOFError, KeyError, subprocess.CalledProcessError) as exc:

267 parser.exit(1, f'Error: {exc}\n')

268 print(json.dumps(result, indent=2))

269

270

271if __name__ == '__main__':

272 main()