pig-gut-profiler.ngl

#!/usr/bin/env ngless
ngless "1.0"
import "mocat" version "1.0"
import "specI" version "1.0"
import "PigGutCatalog" version "1.0"

input = load_mocat_sample(ARGV[1])
RESULTS = ARGV[2]

input = preprocess(input, keep_singles=False) using |read|:
    read = substrim(read, min_quality=25)
    if len(read) < 45:
        discard

mapped = map(input, reference='susScr11')

mapped = select(mapped) using |mr|:
    mr = mr.filter(min_match_size=45, min_identity_pc=90, action={unmatch})
    if mr.flag({mapped}):
        discard

input = as_reads(mapped)

mapped = map(input, reference='pig-gut', mode_all=True)
mapped = select(mapped) using |mr|:
    mr = mr.filter(min_match_size=45, min_identity_pc=95, action={drop})
    if not mr.flag({mapped}):
        discard

counted = count(mapped,
                    features=['OGs'],
                    multiple={dist1},
                    normalization={scaled})
write(counted,
        ofile=RESULTS </> 'eggNOG.traditional.counts.txt',
        auto_comments=[{hash}, {script}])

mapped_refmg = map(input, reference='refmg')
mapped_refmg = select(mapped_refmg) using |mr|:
    mr = mr.filter(min_match_size=45, min_identity_pc=97, action={drop})
    if not mr.flag({mapped}):
        discard

write(count(mapped_refmg,
                features=['specI_cluster']),
    ofile=RESULTS </> 'specI.raw.counts.txt')

write(count(mapped_refmg,
                    features=['specI_cluster'],
                    normalization={scaled}),
    ofile=RESULTS </> 'specI.scaled.counts.txt')