Questão:
Gerando o alinhamento reconstruído do BAM
ShanZhengYang
2018-03-01 10:07:12 UTC
view on stackexchange narkive permalink

Eu tenho um (pequeno) arquivo BAM com os campos CIGAR e MD.

Pergunta 1: Quais ferramentas existem em Python e / ou R para reconstruir o alinhamento entre a referência e a leitura em um BAM? Dado que esta é uma análise muito padrão, acho que deveria haver pacotes R ou bibliotecas Python com funções que realizam essa tarefa com uma entrada BAM ...

Meu objetivo era criar um arquivo * tsv com duas colunas, uma com o alinhamento reconstruído e a outra coluna com a sequência correspondente. Isso seria útil em minha análise. (Eu usaria R data.table ou pandas em Python.)

Aqui está a leitura e a referência:

  CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC |||||||||| |||||| |||||||||||||||||||| CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG  

Aqui está o dataframe de saída desejado (com apenas uma linha de exemplo). Poderíamos usar um R data.table / data.frame ou um DataFrame do pandas. A coluna 1 são os alinhamentos (ou seja, a reconstrução do alinhamento entre a referência e a leitura com base nas cadeias CIGAR e MD), a coluna 2 são as sequências correspondentes. :

  Referência ReadAligned CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC ... ....  

Pergunta 2: Usando as tags CIGAR / MD, pode-se então criar colunas com o número de inserções e o número de exclusões. Esse problema parece muito mais simples, por exemplo, a biblioteca python cigar ou as funções de análise de charutos de pysam

Você poderia esclarecer sua pergunta: Você deseja armazenar os alinhamentos em uma matriz? Ou o formato CIGAR em uma matriz?
@Llopis Um arquivo tsv / csv com duas colunas: a primeira coluna tem os alinhamentos, a segunda coluna tem a seqüência. Isso faz sentido?
Você deseja que o sam2pairwise gere os alinhamentos para você ou deseja analisar as cordas do charuto por conta própria no R / python? Escrever um analisador de charutos é tedioso, mas capaz de fazer. Pode ser mais fácil se você apenas usar sam2pairwise para gerar os alinhamentos para cada leitura e, em seguida, analisar isso em uma tabela.
Agora sua pergunta se reduz a "como contar o número de inserções e eliminações" no charuto, o que é muito mais fácil. Com perl, é apenas um eco de uma linha `78M2I10M5D30M | perl -ne 's / (\ d +) ([ID]) / $ h {$ 2} + = $ 1 / eg; imprimir "$ h {I} \ t $ h {D} \ n" '`. Alguém provavelmente postará respostas python.
@user172818 Tentei editar isso para ser mais específico. Isso é mais organizado? O segundo problema é mais fácil; Estou motivado principalmente pelo primeiro problema, que está relacionado.
Trzy respostas:
heathobrien
2018-03-02 00:42:45 UTC
view on stackexchange narkive permalink

Se tudo o que você deseja é o número de inserções / exclusões, você pode usar isto para extrair essa informação:

  from cigar import Cigarc = c = Cigar ( '20M5I10M3I5M') num_insertions = len (list (filter (lambda x: x [1] == 'I', c.items ()))) inserção_length = sum ([x [0] para x na lista (filtro (lambda x: x [1] == 'I', c.items ()))])  

As exclusões seriam a mesma coisa com 'D' em lugar de 'I'

Existe um método padrão para pegar a string do BAM dentro do python? Parece que isso é analisado automaticamente com pysam
Eu acho que pysam seria a forma padrão
user172818
2018-03-02 12:29:51 UTC
view on stackexchange narkive permalink

Estou respondendo à pergunta 1. No entanto, como não estou familiarizado com o python, estou fornecendo uma solução javascript apenas para mostrar a lógica básica. Para tentar o código no final, você precisa instalar node.js e executá-lo com node this-script.js .

Dado que este é um padrão muito análise

Tendo a acreditar que esta é uma análise muito rara. Eu só vi isso ser feito duas vezes. Converter o CIGAR + MD em alinhamento preenchido é complicado na codificação, pois você deve manter a sequência, o charuto e o MD sincronizados o tempo todo. O código abaixo gera duas strings de uma vez. Pode ser logicamente mais fácil gerar a referência e as strings de consulta separadamente.

Em geral, MD é mal definido. É muito difícil trabalhar com ele. Se possível, evite essa tag. Você pode fazer muitas coisas com a combinação de CIGAR e NM.

  function md2pad (CIGAR, MD, SEQ) {var m, re_cigar = / (\ d +) ([MIDSHNX =]) / g, re_MD = / (\ d +) | (\ ^ [A-Za-z] +) | ([A-Za-z]) / g; var charuto = []; while ((m = re_cigar.exec (CIGAR))! = null) // analisa o CIGAR em uma matriz if (m [2]! = 'H') // não faz nada para o recorte rígido, pois não tem efeito charuto. empurre ([parseInt (m [1]), m [2]]); var k = 0, sx = '', sy = ''; // k: k-ésimo operador de charuto; x: ref; y: consulta // cx / cy: coordenada inicial de charuto [k] em ref / consulta; mx / my: início da op atual MD var cx = 0, cy = 0, mx = 0, my = 0; while ((m = re_MD.exec (MD))! = null) {if (m [2]! = null) {// exclusão da referência var len = m [2] .length - 1; // comprimento da deleção sx + = m [2] .substr (1), sy + = Array (len + 1) .join ("-"); mx + = len, cx + = len, ++ k; // consome uma exclusão D de ref} else {// cópia ou incompatibilidade var ml = m [1]! = null? parseInt (m [1]): 1; // duração desta operação MD while (k < cigar.length && charuto [k] [1]! = 'D') {var cl = charuto [k] [0], op = charuto [k] [1];
if (op == 'M' || op == 'X' || op == '=') {if (my + ml < cy + cl) {// um MD termina no meio do charuto M if ( ml > 0) {if (m [3]! = null) sx + = m [3], sy + = SEQ [my]; else sx + = SEQ.substr (my, ml), sy + = SEQ.substr (my, ml); } mx + = ml, my + = ml, ml = 0; quebrar; } else {var dl = cy + cl - my; sx + = SEQ.substr (my, dl), sy + = SEQ.substr (my, dl); cx + = cl, cy + = cl, ++ k; // consome um operador M mx + = dl, my + = dl, ml - = dl; }} else if (op == 'I') {sy + = SEQ.substr (cy, cl), sx + = Array (cl + 1) .join ("-"); cy + = cl, my + = cl, ++ k; // consumir uma inserção I para ref} else if (op == 'S') {cy + = cl, my + = cl, ++ k; // consome um soft-clipping S} senão lança Error ("MD inconsistente"); // não funciona com N} if (ml! = 0) throw Error ("MD inconsistente"); }} if (cx! = mx || cy! = my) throw Error ("MD inconsistente"); retorno [sx, sy];} var charuto = '16025H3M2I55M1D108M306H'; var MD = '7T15C24T9 ^ C0G2T6T16T3A1A7G21T17G6C1T6T5A4'; SEQ var = 'GATCACAGGCCTATCACCCTATTAATCACTCACGGGAGCTCTCCATGCATCTGGTATTTTTTCGGGGGGGGATGCACGCGATAGCATCGCGGGCCGCTGGAACCGGAGCACCCTATGTCGCAGGATCTGTCTTTGATTCCTACCTCATGCCATTATTAATCGCGCCTA'; var s = md2pad (charuto, MD, seq); console.log (s [0 ]); console.log (s [1]);  
Colin D
2018-11-26 07:59:37 UTC
view on stackexchange narkive permalink

Este programa em C ++ é capaz de realizar a tarefa https://github.com/mlafave/sam2pairwise



Estas perguntas e respostas foram traduzidas automaticamente do idioma inglês.O conteúdo original está disponível em stackexchange, que agradecemos pela licença cc by-sa 3.0 sob a qual é distribuído.
Loading...