Tenho usado o STAR para nossas amostras de RNA-Seq. O arquivo de registro final.out relata a porcentagem de leituras mapeadas exclusivamente junto com a porcentagem de leituras que mapeiam para loci múltiplos (menor ou igual a 10) e a porcentagem de leituras mapeadas para muitos loci (maior que 10). No entanto, quero dividir a parte vários locais em contagens individuais: Lê o mapeamento para 2 locais, 3 locais, 4 locais .. 10 locais.
O NH A tag parece ser usada por STAR . No entanto, uma abordagem ingênua de contagem de leituras resulta em relatar mais número de leituras do que leituras totais.
Por exemplo, meu final.out tem a seguinte aparência:
Velocidade de mapeamento, milhões de leituras por hora | 1403,36 Número de leituras de entrada | 53015978 Comprimento médio de leitura de entrada | 26 LEITURAS EXCLUSIVAS: Número de leituras mapeadas exclusivamente | 368916% de leituras mapeadas exclusivamente | 0,70% Comprimento médio mapeado | 26,45 Número de emendas: Total | 1057 Número de emendas: Anotado (sjdb) | 0 Número de emendas: GT / AG | 802 Número de emendas: GC / AG | 1 Número de emendas: AT / AC | 0 Número de emendas: Não canônico | 254 Taxa de incompatibilidade por base,% | 0,31% Taxa de exclusão por base | 0,00% Comprimento médio da eliminação | 1,45 Taxa de inserção por base | 0,00% Comprimento médio de inserção | 1.00 LEITURAS DE MULTI-MAPEAMENTO: Número de leituras mapeadas para vários locais | 45766732% das leituras mapeadas para vários locais | 86,33% Número de leituras mapeadas para muitos loci | 3757890
% de leituras mapeadas para muitos loci | 7,09% LEITURAS NÃO MAPEADAS:% de leituras não mapeadas: muitas incompatibilidades | 0,00%% de leituras não mapeadas: muito curtas | 5,89%% de leituras não mapeadas: outro | 0,00% Contando o histograma do número de posições que uma leitura mapeia usando pysam:
def get_reads_hist (bam): bam = pysam.AlignmentFile (bam, 'rb') counts = Counter () para consulta em bam.fetch (): nh_count = Counter (dict (query.get_tags ()) ['NH']) contagens + = nh_count retornos contagens resulta em
Counter ({1: 330606, 2: 86772164, 3: 329, 4: 38083, 5: 31, 6: 1094, 7 : 129, 8: 50, 10: 50}) As leituras de contagem 1 são boas, embora não correspondam às contagens em final. out file, pois estou contando uma certa categoria de leituras (digamos, aquelas mapeando apenas para tRNA ), mas as leituras mapeando para 2 locais são altamente superestimadas. Por que isso?