Escrevendo um script collate personalizado

Os arquivos resultantes de cada execução individual são armazenados por padrão em arquivos json (terminando em .results.json). Estes podem ser agrupados em um arquivo mesclado para exibir as mutações de resistência às drogas usando o tb-profiler collate. No entanto, você pode querer personalizar quais informações são extraídas dos arquivos resultantes e exibidas no output agrupado (collated). Esta página fornecerá um esboço de como escrever seu próprio analisador.

Esqueleto do script

Abra seu editor favorito, cole o código a seguir e salve o arquivo como tbprofiler_custom_script.py

#! /usr/bin/env python

import sys
import csv
import json
import argparse
import os
from collections import defaultdict
from tqdm import tqdm
import pathogenprofiler as pp
import tbprofiler as tbprofiler

def main(args):
    bed_file = "%s/share/tbprofiler/%s.bed" % (sys.base_prefix,args.db)
    locus_tag2drugs = tbprofiler.get_lt2drugs(bed_file)

    if args.samples:
        samples = [x.rstrip() for x in open(args.samples).readlines()]
    else:
        samples = [x.replace(args.suffix,"") for x in os.listdir(args.dir) if x[-len(args.suffix):]==args.suffix]

    for s in tqdm(samples):
        data = json.load(open(pp.filecheck("%s/%s%s" % (args.dir,s,args.suffix))))


parser = argparse.ArgumentParser(description='tbprofiler script',formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('--samples',type=str,help='File with samples')
parser.add_argument('--dir',default="results/",type=str,help='Directory containing results')
parser.add_argument('--db',default="tbdb",type=str,help='Database name')
parser.add_argument('--suffix',default=".results.json",type=str,help='File suffix')
parser.set_defaults(func=main)

args = parser.parse_args()
args.func(args)

Isso formará a estrutura básica do analisador. Vamos ver o que cada seção faz.

Linhas 3-11: Aqui importamos bibliotecas que são úteis no script.

Linhas 13: esta linha inicia a definição da função principal (main function). Recebe um argumento na forma de um objeto de argumentos que é construído usando o analisador de argumentos.

Linha 14: A localização do arquivo bed é armazenada na variável bed_file. Isso é construído usando o prefixo do banco de dados fornecido nos argumentos da linha de comando (tbdb por padrão).

Linha 15: Usamos a função tbprofiler.get_lt2drugs para extrair o locus_tag para o mapeamento de droga associado em um dicionário com o seguinte formato:

{
    "Rv0667":["rifampicin"],
    "Rv1484":["isoniazid","ethionamide"],
}

Linha 17-20: Se o argumento --samples foi fornecido, o arquivo é lido e armazenado como uma lista. Caso contrário, a lista é derivada da procura do diretório de resultados (especificado com --dir).

Linhas 22-23: Um loop é construído para percorrer as amostras. Para cada iteração, o arquivo resultante json é carregado em uma estrutura de dicionário. A estrutura de valor-chave é explicada abaixo.

Linhas 26-34: Essas linhas analisam os argumentos fornecidos ao script (por exemplo ---samples). Você pode adicionar novos argumentos e eles aparecerão como propriedades do objeto args na função principal (main function).

Estrutura de dados do resultado

Os arquivos de resultados individuais são lidos em um dicionário na linha 23. Aqui estão algumas (mas não todas) as chaves interessantes no dicionário:

{
    "main_lin": "lineage4",
    "sublin": "lineage4.9",
    "dr_variants": [
      {
        "sample": "por1A",
        "gene_name": "rpoB",
        "chr": "Chromosome",
        "genome_pos": 761155,
        "type": "missense",
        "change": "p.Ser450Leu",
        "freq": 1,
        "nucleotide_change": "761155C>T",
        "locus_tag": "Rv0667",
        "gene": "rpoB",
        "_internal_change": "450S>450L",
        "drug": "rifampicin",
        "literature": "10.1128/AAC.01093-18",
        "confidence": "high"
      }
    ],
    "other_variants": [
      {
        "sample": "DRR034528",
        "gene_name": "gyrA",
        "chr": "Chromosome",
        "genome_pos": 7362,
        "type": "missense",
        "change": "p.Glu21Gln",
        "freq": 1,
        "nucleotide_change": "7362G>C",
        "locus_tag": "Rv0006",
        "gene": "gyrA",
        "_internal_change": "21E>21Q"
      }
    ]
}

As mutações são divididas em duas listas “dr_variants” e “other_variants”. As “dr_variants” estão todas associadas à resistência às drogas. As outras variantes (“other_variants”) são mutações em genes resistentes às drogas, mas que não foram (ainda!) associadas à resistência. Algumas das outras mutações podem ser mutações filogenéticas que não têm nada a ver com resistência, enquanto outras podem ser, mas podem não ter evidências suficientes ainda.

No momento, se você executar o script e apontar para um diretório com alguns arquivos de resultado usando --dir ele lerá todos os dados, mas não fará nada com eles.

Roteiro de exemplo de trabalho

Como exemplo, estenderei sua funcionalidade para gerar um arquivo CSV simples com o nome da amostra e o tipo de resistência à droga.

Primeiro, criarei um objeto de gravação CSV. Adicione o seguinte código antes do for loop.

OUT = open(args.out,"w")
writer = csv.DictWriter(OUT, fieldnames = ["sample","dr-class"])
writer.writeheader()

Feche o identificador de arquivo após o término do loop usando

OUT.close()

Usamos a propriedade args.out ao criar o filehandle. Para que isso funcione, temos que adicionar um argumento extra ao analisador de argumentos. Adicione o seguinte para o código junto com os outros argumentos:

parser.add_argument('--out', type=str, help='Name of CSV output', required=True)

Agora vamos adicionar naquele código que determina a classe de resistência às drogas. Por padrão, o tb-profiler classifica apenas as amostras em: sensível, MDR, XDR e droga resistente. Vamos estender isso para encontrar pré-MDR e pré-XDR.

Primeiro, podemos criar um conjunto contendo uma lista dos medicamentos aos quais a amostra é resistente. Podemos fazer isso usando uma compreensão de lista que é convertida em um conjunto. Adicione o seguinte código dentro do for loop:

resistant_drugs = set([d["drug"] for d in data["dr_variants"]])

Definições das classes de drogas

As amostras podem ser classificadas em diferentes tipos usando as seguintes definições

Tipo	Droga resistência
Sensível	Sem droga resistance
Pré-MDR	Rifampicina ou isoniazida
MDR	Rifampicina e isoniazida
Pré-XDR	MDR e (qualquer fluoroquinolona ou qualquer segunda-linha injetável)
XDR	MDR e (qualquer fluoroquinolona e qualquer segunda-linha injetável)
Outras	Resistência para qualquer droga mas nenhuma das categorias acima

Codificando isso

Podemos ver imediatamente que as definições são apenas algumas operações booleanas. Podemos armazenar os valores booleanos procurando as drogas no conjunto resistant_drugs que criamos antes. Primeiro, criaremos uma lista de fluoroquinolonas (FLQ) e de segunda linha injetáveis (SLI). Adicione este código em qualquer lugar na função principal, mas antes do loop:

FLQ_set = set(["moxifloxacin","levofloxacin","ciprofloxacin","ofloxacin"])
SLI_set = set(["amikacin","capreomycin","kanamycin"])

Em seguida, crie os valores booleanos. Coloque-os dentro do loop depois de criar a variável de data:

resistant_drugs = set([d["drug"] for d in data["dr_variants"]])
rif = "rifampicin" in resistant_drugs
inh = "isoniazid" in resistant_drugs
flq = len(FLQ_set.intersection(resistant_drugs)) > 0
sli = len(SLI_set.intersection(resistant_drugs)) > 0

Para testar a rifampicina e a isoniazida, verificamos simplesmente se elas existem no resistant_drugs. Para verificar a resistência FLQ e SLI, encontramos as interseções entre FLQ_set e resistant_drugs e count. Se o valor for maior que 0, então temos resistência nessa categoria.

Agora estamos prontos para realizar a construção das classes usando nossas variáveis booleanas. Podemos percorrer a tabela acima usando uma série de instruções if-else.

if len(resistant_drugs)==0:
    drtype = "Sensitive"
elif (rif and not inh) or (inh and not rif):
    drtype = "Pre-MDR"
elif (rif and inh) and (not flq and not sli):
    drtype = "MDR"
elif (rif and inh) and ( (flq and not sli) or (sli and not flq) ):
    drtype = "Pre-XDR"
elif (rif and inh) and (flq and sli):
    drtype = "XDR"
else:
    drtype = "Other"

Armazenamos a classe de resistência aos medicamentos na variável drtype. Agora, tudo o que resta fazer é gravar isso no arquivo CSV. Podemos fazer isso usando:

writer.writerow({"sample":s, "dr-class":drtype})

Seu script agora deve ter a seguinte aparência:

 #! /usr/bin/env python

import sys
import csv
import json
import argparse
import os
from collections import defaultdict
from tqdm import tqdm
import pathogenprofiler as pp
import tbprofiler as tbprofiler

def main(args):
    bed_file = "%s/share/tbprofiler/%s.bed" % (sys.base_prefix,args.db)
    locus_tag2drugs = tbprofiler.get_lt2drugs(bed_file)

    if args.samples:
        samples = [x.rstrip() for x in open(args.samples).readlines()]
    else:
        samples = [x.replace(args.suffix,"") for x in os.listdir(args.dir) if x[-len(args.suffix):]==args.suffix]

    FLQ_set = set(["moxifloxacin","levofloxacin","ciprofloxacin","ofloxacin"])
    SLI_set = set(["amikacin","capreomycin","kanamycin"])

    OUT = open(args.out,"w")
    writer = csv.DictWriter(OUT, fieldnames = ["sample","dr-class"])
    writer.writeheader()

    for s in tqdm(samples):
        data = json.load(open(pp.filecheck("%s/%s%s" % (args.dir,s,args.suffix))))
        resistant_drugs = set([d["drug"] for d in data["dr_variants"]])

        resistant_drugs = set([d["drug"] for d in data["dr_variants"]])
        rif = "rifampicin" in resistant_drugs
        inh = "isoniazid" in resistant_drugs
        flq = len(FLQ_set.intersection(resistant_drugs)) > 0
        sli = len(SLI_set.intersection(resistant_drugs)) > 0

        if len(resistant_drugs)==0:
            drtype = "Sensitive"
        elif (rif and not inh) or (inh and not rif):
            drtype = "Pre-MDR"
        elif (rif and inh) and (not flq and not sli):
            drtype = "MDR"
        elif (rif and inh) and ( (flq and not sli) or (sli and not flq) ):
            drtype = "Pre-XDR"
        elif (rif and inh) and (flq and sli):
            drtype = "XDR"
        else:
            drtype = "Other"

        writer.writerow({"sample":s, "dr-class":drtype})

    OUT.close()

parser = argparse.ArgumentParser(description='tbprofiler script',formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('--out', type=str, help='Name of CSV output', required=True)
parser.add_argument('--samples', type=str, help='File with samples')
parser.add_argument('--dir', default="results/", type=str, help='Directory containing results')
parser.add_argument('--db', default="tbdb", type=str, help='Database name')
parser.add_argument('--suffix', default=".results.json", type=str, help='File suffix')
parser.set_defaults(func=main)

args = parser.parse_args()
args.func(args)

Tente executá-lo e verifique se o output parece como deveria:

python tbprofiler_custom_script.py --dir /path/to/dir --out test.csv

Se estiver funcionando conforme o esperado, tente adicionar outras informações ao output. Por exemplo, você pode extrair a sublinhagem acessando data["sublin"]. Dê uma olhada na estrutura do arquivo de resultado (você pode abri-lo com o Firefox para direcioná-lo a um formato mais legível). Esperançosamente, você pode usar este modelo como um padrão para seus próprios usos. Existem alguns exemplos disponíveis aqui se você estiver interessado. Se você tiver alguma dúvida ou se tiver um script legal que gostaria de disponibilizar para a comunidade, não hesite em me contatar por e-mail ou pelo GitHub. Obrigado por ler!