use picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads

READ_LIST_FILE (Arquivo) Lista de leitura Arquivo contendo leituras que serão incluídas ou excluídas do arquivo OUTPUT SAM ou BAM. Valor padrão: nulo.

(1) Não será muito rápido, mas você pode fornecer ao grep um arquivo de QNAMES.

  samtools view file.bam | grep -f 'qnames.txt > subset.sam  

onde qnames.txt tem

  EXEMPLO: QNAME1EXAMPLE: QNAME2EXAMPLE: QNAME3EXAMPLE: QNAME4EXAMPLE: QNAME5  

(2) Isso seria um pouco mais complicado, mas você pode dar um exemplo onde o grep pode estar com o QNAME correto?

(3) Para manter o cabeçalho BAM Eu usaria samtools -H file.bam > header.txt para obter o cabeçalho e, em seguida, cat o arquivo de cabeçalho com o arquivo grep'ed

Abaixo está um pequeno código python mostrando uma maneira ingênua, mas gerenciável de fazer isso (NB, espero que grep seja mais rápido, embora seja irritante fazer com que ele exiba o cabeçalho):

  import pysamqnames = set (...) # ler nomes go herebam = pysam.AlignmentFile ("algum arquivo de entrada.bam") obam = pysam. AlignmentFile ("algum arquivo de saída.bam", "w", template = bam) para b em bam.fetch (until_eof = True): se b.query_name em qnames: obam.write (b) bam.close () obam. close ()  

Agora isso vai funcionar, mas se você tiver muitos nomes para olhar, vai acabar ficando bem lento (assim como grep -f. .. ). Se você precisar selecionar uma grande lista de leituras por nome, uma estratégia mais eficiente é:

1. Classifique o arquivo BAM por nome ( samtools sort -n ou picard) , que pode ser feito com vários tópicos.
2. Classifique a lista de nomes lidos para corresponder (isso será mais complicado do que se poderia ingenuamente pensar, uma vez que samtools e picard nomearão de forma diferente, então certifique-se de colocar alguma reflexão sobre isso).
3. Execute o mesmo tipo de iteração acima, mas apenas compare com o nome mais lido em sua lista da etapa 2 (removendo esse elemento da pilha após uma correspondência).

Executar um grep em $ n $ alinhamentos com $ m $ qnames daria a você $ O (mn) $ operações. Se você decidir classificar seu BAM e qnames primeiro, pode reduzi-lo para $ O (\ min \ {m, n \}) $ retirando leituras e qnames suas pilhas. Você também não lê todo o BAM de uma vez dessa forma, limitando o consumo de memória. Tenho quase certeza de que Picard também faz algo assim.

Aqui está um exemplo em Python que escrevi exatamente para esse propósito:

  import pysamdef filter_qname ( bamfile, idfile, outfile): '' 'Filtrar leituras em um arquivo SAM / BAM por seus nomes de consulta Parâmetros ---------- bamfile: str Caminho do arquivo BAM classificado idfile: str Caminho do arquivo de texto contendo qnames para manter outfile: str Arquivo de saída para gravar em '' '# lido no arquivo BAM classificado por nome bam = pysam.AlignmentFile (bamfile, "rb") if outfile.endswith ('. sam '): output = pysam.AlignmentFile (outputfile, "w", template = bam) elif outfile.endswith ('. bam'): output = pysam.AlignmentFile (outputfile, "wb", template = bam) else: raise ValueError ("Formato de arquivo de saída desconhecido para ʻoutfile`: {} ". format (outfile)) # lido em IDs a serem removidos (list (set (...)) para manter apenas IDs únicos) ids = list (set ([l.rstrip () para l em aberto (idfile , 'r'). readlines ()])) # classificar IDs para corresponder BAM para processamento eficiente (função destrutiva) ids.sort (key = str.lower, reverse = True) n_ids = len (ids) # variável para garantir que o BAM seja classificado pelo nome da consulta last_q = Nenhum lê = bam.fetch (until_eof = True ) ler = próximo (leituras) enquanto Verdadeiro: # se o nome de leitura for maior que o topo da pilha atual se read.query_name > ids [-1]: # se este for o último ID se n_ids == 1: # escrever leituras restantes novo arquivo e sair do loop while output.write (ler)
para r em bam: output.write (read) break # caso contrário, pop esse ID, tente novamente else: ids.pop () n_ids - = 1 # ignorar leituras que correspondem ao topo da pilha elif read.query_name == ids [- 1]: read = next (lê) # se o nome de leitura for menor que o topo da pilha, escreva-o e siga em frente, senão: output.write (read) read = next (lê) output.close ()  

Eu implementei em meu próprio pacote de ferramentas de bioinformática diversas aqui se você quiser vê-lo em ação

Você pode usar o SAMsift:

  samsift -i file.bam -0 'q = open ("qnames .txt "). read (). splitlines () '-f' QNAME in q ' 

-i file.bam especifica o arquivo BAM de entrada. -0 'q = open ("qnames.txt"). read (). splitlines ()' carrega sua lista de QNAMEs. -f 'QNAME in q' especifica o filtro para alinhamentos.

Com o SAMsift, pode-se usar qualquer expressão Python para filtragem. A desvantagem é que ele é mais lento do que as outras ferramentas.

A maioria das soluções postadas até agora são lentas ou podem produzir resultados falsos em certos casos. (Não testei Picard, que presumo que funcione conforme o planejado. Mas, pessoalmente, tenho tendência a recuar quando tenho que configurar uma JVM!)

Vamos abordar o problema de um ângulo diferente:

Boas funções hash têm complexidade ~ O (1) para uma pesquisa. Se pudermos criar uma tabela hash de nossos QNAMEs, a pesquisa de um arquivo BAM deve ser O (n) no número de registros. Os desenvolvedores Python têm muito orgulho da função hash em cpython, então vamos usá-la:

  #! / Usr / bin / env python3import syswith open (sys.argv [1], 'r' ) como indexfile: ids = set (l.rstrip ('\ r \ n') para l em indexfile) para linha em sys.stdin: qname, _ = line.split ('\ t', 1) se qname em ids : sys.stdout.write (line)  

Com um conjunto de teste consistindo em um arquivo BAM de aproximadamente 3 GB e um qnames.txt com aproximadamente 65 mil entradas, executando samtools view input.bam | ./idfilter.py qnames.txt leva cerca de 6,5 segundos no meu servidor de teste. Em comparação, canalizar o SAM para fgrep -f qnames.txt (usando GNU grep) leva cerca de 8,5 segundos (mas pode produzir resultados espúrios).

(Por que grep produz resultados espúrios? Se você tiver um QNAME foo em seu qnames.txt , mas seu arquivo BAM também contém fooX e fooY , eles serão correspondidos por fgrep . A solução é ancorar cada um de seus padrões de pesquisa entre o início da linha e uma guia, por exemplo, ^ foo \ t , mas então você tem que usar o grep padrão que terá que construir um NFA para cada padrão, e não o fgrep que implementa Aho-Corasick, e você será ordens de magnitude mais lento.)

Reescrevi meu pequeno script Python em Rust usando std :: Collections :: HashMap para a função hash e usando a caixa rust_htslib para ler e escrever diretamente do BAM, e isso é mais que duas vezes mais rápido que o Python, mesmo ao ler e escrever BAM. No entanto, meu Rust não é realmente adequado para consumo público ...