Em formação

É possível realinhar os arquivos SAM / BAM?


Estou realizando algumas comparações metagenômicas usando a abordagem de sequenciamento do genoma completo. Gostaria de filtrar as leituras de meus arquivos de maneira sequencial, portanto, gostaria de alinhar com um primeiro genoma de referência e, em seguida, usar o arquivo SAM / BAM gerado para extrair as leituras mapeadas com visualização de samtols -F4, regenerar fastq arquivos e alinhá-los a outro genoma de referência. Isso diminui o tamanho da entrada e acelera o procedimento. Usando BWA, porém, recebo um erro de sincronização de arquivo ao fazer o segundo alinhamento:

[mem_sam_pe] leituras emparelhadas têm nomes diferentes: "HISEQ: 90: C3UNJACXX: 1: 1107: 16388: 61141", "HISEQ: 90: C3UNJACXX: 3: 1303: 2008: 58095" [mem_sam_pe] [mem_sam_pe] leituras emparelhadas têm diferentes nomes: "HWI-ST1437: 64: C3UM1ACXX: 4: 1214: 3833: 62731", "HWI-ST1437: 64: C3UM1ACXX: 4: 1202: 18519: 45906"

Isso é estranho, já que os arquivos fastq foram gerados pelo arquivo BAM filtrado, os samtools deveriam saber como manter a sincronização dos arquivos. Minha pergunta é: é possível realinhar um arquivo BAM? talvez não com BWA, pelo menos com outra ferramenta?

Obrigada.

O procedimento que realizei é:

# align bwa mem -R     -o  # convert samtools view -Sb  >  # sort samtools sort  -o  # seleciona a visualização das samtools mapeadas -h -F 4  -o  # converter para fastq samtools fastq -1  -2   # realinhar bwa mem -R     | samtools sort -o 

Parece que resolvi usando PICARD: extraindo as leituras alinhadas, mas sem classificar, usei

java -jar picard.jar SamToFastq I = FASTQ = SECOND_END_FASTQ = bwa mem -R    

A partir desamtools flagstata saída parece boa.


É possível realinhar os arquivos SAM / BAM? - Biologia

Este site é usado apenas para teste. Visite: https://www.biostars.org para fazer uma pergunta.

Extraí a região específica de um arquivo BAM gerado pelo alinhamento fastq (Illumina / sequenciamento do genoma completo) em hg38 e converti para leituras fastq correspondentes dessa região por Picard (SamTofastq) e tentei realinhar para outra (minha) sequência de referência por bwa mem. Mas, parece que o bwa mem não consegue entender que o sequenciamento é como final pareado, uma vez que freqüentemente diz "pule a orientação FF, pois não há pares suficientes" durante o mapeamento, isso também aconteceu para outras orientações (FR, RF e RR), no entanto, há 1259352 leitura emparelhada. Posteriormente, a porcentagem emparelhada corretamente é muito baixa. Você poderia, por favor, me ajudar como posso resolver o problema?

Editar e atualizar: Eu também descobri que durante a conversão de bam para fastq por SamTofastq (picard), ele me retornou: Erro de validação de SAM: ERRO: Encontrado 13548 companheiros desemparelhados. Acho que pode não ser um problema real, na verdade, parece que Picard apenas considerou a leitura emparelhada para criar o arquivo Fastq. Então, o verdadeiro problema é com bwa mem, certo? Por favor, compartilhe sua ideia comigo?

Aqui estão os comandos usados:

Mostre-nos os cabeçalhos dos arquivos extraídos do alinhamento original. Eles provavelmente perderam as informações de origem de leitura.

Embora seja trivial, esse erro surgiu no passado quando as pessoas tentaram usar o mesmo arquivo (por exemplo, R1 / R1 em vez de R1 / R2) ao fazer o segundo alinhamento. Algo para verificar.

Desculpe, o arquivo está no cluster e não pode ser transferido. Você poderia me dizer onde estão as informações de orientação de leitura no cabeçalho BAM, o que devo procurar? Fiz alinhamento com dois arquivos fastq R1 / R2.

O que zcat filename_R1.fq.gz | head -4 e zcat filename_R2.fq.gz | head -4 show? Se seus arquivos estiverem descompactados, substitua zcat por cat no comando.

O nome de leitura em dois arquivos fastq gerados por SamTofastq é o mesmo, exceto para / 1 e / 2. Além disso, a contagem de leitura de dois arquivos fastq é a mesma. Então, o problema é outro lugar. Além disso, alguém diz "bwa espera que a primeira leitura do arquivo um seja o companheiro da primeira leitura do arquivo 2. Essa suposição obviamente não se mantém em todo o arquivo, porque classificar por coordenadas e jogar fora as leituras estragam isso" neste publicar. O que você acha?

Você deve postar o exato comandos que você usou para extrair as leituras do arquivo bam e, posteriormente, mapear os arquivos fastq novamente

Você também deve fazer o melhor para fornecer as informações solicitadas quando alguém solicitar informações adicionais.

Eu adicionei os comandos usados ​​à postagem original.

O alinhamento original ou a conversão de samtofastq eliminou os identificadores que mostram qual das leituras é de R1 e qual é de R2. Como você nos garante que as leituras de R1 / R2 estão em sincronia, você pode corrigir esse problema usando a seguinte ferramenta do pacote BBMap:

adicione uma das duas opções abaixo. Qualquer um deve funcionar.

é a saída de head -4 de dois arquivos fastq (desculpe se a resolução não for alta o suficiente), como você pode ver o nome lido em ambos os arquivos fastq é o mesmo e está junto com / 1 e / 2, respectivamente, então o nome lido e seu identificador (/ 1 e / 2) está correto. É possível que o nome lido ou os identificadores relacionados fiquem errados em todos os arquivos fastq? você ainda sugere o uso de repair.sh?

Me desculpe. Eu quis dizer reformat.sh em vez de repair.sh (corrigido agora).

BTW: não vemos seus exemplos fastq. Você pode copiar e colar as linhas (não poste como uma imagem) e, em seguida, formatar usando o botão de código.


Suas leituras têm / 1 e / 2. Portanto, é possível que simplesmente não haja leituras suficientes para bwa fazer essa estimativa.

Sim, as leituras têm / 1 e / 2 com base no cabeçalho de dois arquivos fastq. Existe 1259352 leitura emparelhada, não é suficiente para bwa mem na sua opinião? Eu também testei o bowtie2, mas a taxa geral de alinhamento foi muito ruim (cerca de 5,4%). Com bwa mem, a porcentagem de mapeamento é de cerca de 95%, mas o emparelhado corretamente é muito baixo (cerca de 10%). Por favor, compartilhe-me alguma sugestão e ideia?

Você está recebendo aquele erro / aviso com a) dados originais E b) leituras extraídas (

13K) para uma região específica. Achei que fosse apenas o último e, portanto, pode não haver dados suficientes.

Por que você está escolhendo dados como este BTW? A sua referência é apenas para esta região menor? Se for para todo o genoma, por que você não pode usar todos os dados para alinhar?

/ 1 e / 2 não são usados ​​há algum tempo. São dados antigos?

Obrigado por acompanhar o assunto. Em relação às suas perguntas, recebo o erro / aviso com as leituras extraídas, há mais de 1M de leituras PE (não 13K) no arquivo bam extraído, conforme mencionei em meu comentário anterior. Assumindo que o Bwa mem pode prever a distribuição errada do tamanho do inserto, usei CollectInsertSizeMetrics (de Picard) para prever a distribuição do tamanho do inserto no arquivo bam extraído, no entanto, não tenho certeza de como esses números devem ser alimentados para o bwa mem via I flag. Você poderia, por favor, me mostrar um exemplo?

2) sim, minha própria referência está relacionada apenas a esta pequena região. Na verdade, esta região é HLA e eu gostaria de visualizar os alelos relatados da tipagem HLA (de dados de sequenciamento do genoma inteiro) no IGV. Então, extraí essa região, converti em arquivos fastq e tentei realinhar para sequências genômicas hla como referência. Será ótimo se você me compartilhar qualquer ideia / sugestão sobre este assunto.

Quando você realinha estes, você não deveria realinhar para as sequências HLA FASTA de, por exemplo, IMGT?

Olá, Kevin, não realinhe a sequência fasta do HLA! por favor, diga-me qual é a sequência de referência para realinhar? Estou ficando mais confuso

Acho que mostrei a você (em algum outro tópico) um trabalho publicado por meus colegas da UCL (?) - nele, há alguns passos a seguir.

Sim, sua média é LOHHLA, eu acho. No entanto, no momento, gostaria de focar apenas na visualização de alelos HLA (a tipagem de HLA do sequenciamento do genoma inteiro foi feita por outra pessoa). Depois de resolver o problema atual, eu certamente testarei o LOHHLA, mas ele requer o Novalign que não é gratuito !. Portanto, agora apenas para a visualização dos alelos HLA digitados no IGV, realinhei as leituras fastq extraídas para a sequência fasta HLA como referência por bwa mem e enfrentei um mapeamento emparelhado adequado muito baixo. Considerando que as leituras em dois arquivos fastq gerados são pareadas com base na saída do cabeçote, que mostrei nos comentários anteriores deste post, acho que bwa mem pode estimar incorretamente a distribuição do tamanho da inserção. É por isso que eu calculei por Picard, mas não tenho certeza de como alimentá-lo para bwa mem através da opção I. Kevin, você poderia gentilmente me ajudar a resolver o problema atual?

Eu tentei replicar esse pipeline sozinho (LOHHLA), mas acabei ficando muito ocupado. Na verdade, tive uma reunião com o desenvolvedor (Dr. McGranahan), mas não consigo me lembrar dos detalhes mais sutis, neste estágio. Sim, NovoAlign foi uma escolha estranha, e continua sendo. No entanto, acredito que as versões mais antigas sejam gratuitas?

A tipagem HLA tem sido uma área obscura e estranha em bioinformática nos últimos anos. Minha primeira tentativa foi em 2015 com os dados do Roche 454 e fiz tudo manualmente - foi doloroso.

Você pode dar uma olhada na sugestão do Wouter, talvez? - Veja abaixo.

Obrigado. Como descobri que existem vários programas para digitação HLA, não há necessidade de trabalhar manualmente atualmente, mas prefiro trabalhar manualmente em uma amostra para um melhor aprendizado. Por favor, compartilhe qualquer guia / tutorial bem documentado, se você souber. A ferramenta sugerida por Wouter é outro programa de tipagem HLA como eu disse a vocês, no momento eu tenho que me concentrar apenas na visualização dos alelos HLA.

Faça login antes de adicionar sua resposta.

O uso deste site constitui aceitação de nosso Acordo de Usuário e Política de Privacidade.


Estimando comprimento de inserção de leitura emparelhada de arquivos SAM / BAM

Eu escrevi um único script Python para estimar o comprimento de inserção de leitura de extremidade pareada (ou comprimento de fragmento) das informações de mapeamento de leitura (ou seja, arquivos SAM / BAM). O algoritmo é simples: verifique o campo TLEN no formato SAM, descarte leituras de pares cujos pares estão muito distantes e use-as para estimar a média e a variância do comprimento da inserção.

Este script também é capaz de fornecer uma distribuição detalhada do comprimento de leitura e amplitude de leitura para sua conveniência. Consulte o uso detalhado abaixo.

Este script é distribuído no GitHub agora.

Observação : De acordo com a definição do SAM, o intervalo de leitura "é igual ao número de bases da base mapeada mais à esquerda até a base mapeada mais à direita". Este intervalo é a distância entre duas leituras em uma leitura emparelhada MAIS 2 vezes o comprimento de leitura. A extensão de leitura é diferente de "mate-inner-distance" no Tophat (opção -r), que mede apenas a distância entre duas leituras em uma leitura de extremidade pareada.

47 comentários:

Olá, wei, estou usando seu script para calcular o tamanho da inserção, mas está dizendo isto gaurav @ gaurav-OptiPlex-980:

/ Desktop $ samtools visualizar bowtie_out.nameSorted.PropMapPairsForRSEM.bam | head -n 1000000 | python getinsertsize.py
1M.
Comprimento de leitura: média 100,0, STD = 0,0
Comprimento de leitura possível e suas contagens:
<100: 1000000>
Nenhuma leitura de fim emparelhada qualificada encontrada. São leituras single-end?

Você pode me dizer como verificar se minhas leituras são de ponta única ou ponta emparelhada.

Provavelmente, essas leituras são leituras de extremidade única. Você pode verificar algumas linhas de seu arquivo BAM e examinar o 7º e o 8º campos de cada linha. Se forem leituras de extremidade única, você verá marcas como & quot0 * & quot, o que significa que não há informações emparelhadas correspondentes.

Eu verifiquei se minhas leituras estão emparelhadas. Na verdade, eu construí um transcriptoma denovo e mapeei essas leituras de volta para o transcriptoma e fiz o arquivo BAM. Portanto, este script funcionará neste caso também ou foi escrito especificamente para um assembly baseado em referência.

Também recebo uma mensagem dizendo:
Nenhuma leitura de fim emparelhada qualificada encontrada. São leituras single-end?

mas minhas leituras são emparelhadas. Eu posso ver o emparelhamento com igv. Você tem alguma sugestão?

é possível que seu formato bam seja uma versão mais antiga. você pode colar algumas linhas do seu arquivo bam? use o & quotsamtools ver test.bam | head -n 10 & quot imprime as primeiras 10 linhas de seu registro bam (test.bam).

Olá. Obrigado por compartilhar o script! A propósito, posso perguntar por que você usou o campo PNEXT (8ª coluna em um arquivo SAM) em vez do campo TLEN (9ª coluna em um arquivo SAM). Do meu entendimento, TLEN indica um tamanho de inserção inferido.

Na maioria dos casos, você pode usar qualquer uma das formas para calcular o tamanho da pastilha.

Obrigado pelo script Wei.

Tive que modificá-lo um pouco para fazê-lo funcionar, e acho que esse é o mesmo problema que os comentadores acima de mim provavelmente estavam tendo. O campo NH: i: é um campo opcional e muitos mapeadores de leitura curta realmente não o incluem, incluindo alguns mapeadores populares.

Como seu script exigia que o campo opcional estivesse presente e igual a 1 para contar essa leitura, todas as leituras estavam sendo ignoradas (dando essa mensagem: nenhuma leitura de fim emparelhada qualificada). Comentei essas linhas pertencentes a este campo NH: i: e, em seguida, o script foi executado com êxito e contou as leituras.

Você pode encaçapar o código aqui. Ainda estou recebendo o erro.

Apenas comente da linha 63 à linha 66 e deve estar tudo bem com os problemas & quotNH: i & quot. Eu também atualizei o script para que você possa experimentar o mais recente.

Estou tendo o erro a seguir:

$ samtools ver 10081.marked.realigned.recal.bam | head -n 1000000 | python getinsertsize.py
Arquivo & quotgetinsertsize.py & quot, linha 45
imprimir (str (nline / 1000000) + & # 39M. & # 39, arquivo = sys.stderr)
^
SyntaxError: sintaxe inválida

Alguém sabe qual é o problema?

qual é a sua versão de python? Funciona no 2.7, mas não tenho certeza se funciona no 2.6 ou inferior.

Ok, atualizei minha versão e obtive o seguinte resultado:

1.0M.
Comprimento de leitura: média 101,0, STD = 0,0
Comprimento de leitura possível e suas contagens:
<101: 951678>
Amplitude de leitura: média 0,0, STD = 0,0

Tenho certeza de que essas são leituras finais emparelhadas, por que obtenho intervalo 0?

Aqui está um pouco do meu arquivo bam também:

HWI-1KL153: 29: D1N72ACXX: 8: 1209: 12276:? 12090 1187 CHRM 9 60 1 01M = 200 292 GTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTA TTTTCGTCTGGGGGGTGTGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAG @ & GTB & gtEBBACB @ ?? @ @ @@ & gt BC??? B & GTB ADCABBCBCADCC & gt @ DC = & lt @ A @ C & lt @ DBCABBBC:??? = & Amp & lt @ & gtFFCB @ CF9BFBE @ G & lt = CA BB9 & gt & GTB X 0:? I: 1 X1: i: 0 BD: Z: KKOLLKNLJLILLJJIJKLMJLLLJMMIMMONLKLMNNNONNKBMNMMKKBBLNMOLON JJJJMKJKOOKNNNONKNPOOLNPNPONOOPQQQQQTTMNNN MD: Z: 101 RG: Z: Exome_10081 XG: i: 0 BI: Z: NNRPQNPPMOKPPNMNLNORMPLPMQOLQPRQLOLORQRQQQNFPPQQONGGPRQRPRQMMMMRNOO RRORQQOQPRSRSOQSRPRRSTSUUTUUWUQPRQ AM: i: 37 NM: i: 0 SM: i: 37 XM: i: 0 XO: i: 0 MQ: i: 60 XT: Um :VOCÊ

Você se importaria de me enviar as primeiras 10.000 linhas do seu arquivo bam (use head -n 10000)? Posso ver o que acontece. Você pode compactar o arquivo de texto e enviá-lo para li.david.wei AT gmail.com. Obrigado!

Olá, David, eu enviei para você. Informe se você recebeu.

Eu olhei seus dados. O problema é que há muitas leituras com tamanho de inserção 0, o que confunde o programa e filtra todas as leituras. Na verdade, este é um bug no meu programa, e eu atualizei o código-fonte e ele deve funcionar bem agora. Obrigado e tente!

PS: você precisa executar o programa como o & quotsamtools view [arquivo BAM] | head -n 1000000 | getinsertsize.py - & quot agora em vez de & quotsamtools view [arquivo BAM] | head -n 1000000 | getinsertsize.py & quot.

parece funcionar bem, obrigado por suas respostas muito rápidas!

Para o meu único exemplo, tenho o seguinte

Extensão de leitura: média 364,3316085392387, STD = 123,087158684682

O STD parece ser bastante alto, mas eu não tenho certeza. Alguém tem experiência com o que deveria ser a DST esperada?

Provavelmente, usar mais leituras terá uma estimativa melhor do intervalo. Além disso, parece que muitas de suas primeiras leituras de 10M são de chrM, o que pode não ser uma boa escolha. Você pode tentar usar leituras de outros cromossomos?

Tentei usar o bam inteiro e obtive praticamente o mesmo número (

100 milhões de leituras). Portanto, este parece ser um resultado real.

Eu também estava me perguntando se o script pode ser modificado para que o arquivo de distribuição de amplitude seja gerado de forma que as contagens de leitura sejam para caixas de 10 bp. Eu sou novo em python e a ajuda seria muito apreciada.

Você provavelmente pode carregar o arquivo em R e verificar o histograma. Isso deve lhe dar uma visão geral intuitiva das distribuições span.

Desculpe, a seta deve estar apontando para o sinal de igual de & quotfile = sys.stderr & quot

Oi! e muito obrigado pelo seu roteiro.

Eu sugeriria melhorar seu desempenho para eliminar & # 39.keys () & # 39 nas linhas 61 e 76.
if readlen em plrdlen.keys () se tornaria if readlen em plrdlen.

A propósito, você poderia definir exatamente o que entende por extensão de leitura?

Obrigado pela sua resposta. extensão de leitura é a distância entre duas leituras em uma leitura de extremidade emparelhada mais 2 * comprimento de leitura. É igual ao tamanho de inserção do fragmento.

Ótimo roteiro !! Muito obrigado da Malásia.

Quer saber se você saberia integrar o intervalo de leitura (assumindo que esse intervalo de leitura se refere ao comprimento de inserção) com o alinhamento tophat2?

Mais especificamente, um dos parâmetros tophat2 é & quot-r & quot (também conhecido como mate-inner-distance) http://tophat.cbcb.umd.edu/manual.html

Portanto, se, por exemplo, o valor do intervalo de leitura do seu script for 200, você saberia o que usar para o parâmetro tophat2 -r ??

Você precisaria considerar o comprimento de leitura da amostra (exemplo 100bp) para determinar um valor tophat2 -r correto.

o valor tophat2 -r é o intervalo de leitura menos 2 * comprimento de leitura, portanto, depende do comprimento de leitura. Como é a distância entre duas leituras em uma leitura emparelhada, é chamada de & quotmate-inner-distance & quot.

Hmmmmmmm. seguindo a partir da pergunta,

Incluir o parâmetro -r parece simples graças ao seu esclarecimento. Agora que incluí o valor (intervalo de leitura - 2 * comprimento de leitura) no Tophat. Para adicionar uma camada de complicação, se eu quiser também incluir o valor do desvio padrão nos parâmetros de mapeamento. Não posso simplesmente usar o valor do seu script, estou correto?

Por exemplo, a saída do seu script
amplitude de leitura: 242 SD: 90

Portanto, o valor de -r seria: 42 (242-2 * 100). No entanto, o mesmo não pode ser feito com SD, não podemos apenas menos um valor .. (se meu entendimento estiver correto)

A pergunta que estou fazendo. seria possível obter um valor SD ajustado com base no comprimento do intervalo de leitura modificado ?? Estou pensando em simplesmente subtrair 200 de cada comprimento de leitura para obter uma nova estimativa de SD.

Eu estava olhando para o script e vi onde seria um lugar para inserir o valor readlength * 2 e obter um novo SD correto. mas eu não consigo apontar onde.

Olá, o SD é o mesmo (90 em sua amostra). Isso ocorre porque para uma variável aleatória xe uma constante a, SD (x) = SD (x-a). O span de leitura vem diretamente da definição do campo [8], que já inclui o comprimento de leitura.

Obrigado pelo insight Wei!

Infelizmente, ainda não entendi bem. Acho que preciso me aprofundar nas estatísticas, bem como no campo Sam :)

Ainda não conseguirei ver por que, se ajustarmos o intervalo de leitura (subtraindo 200), o SD permanecerá o mesmo ??

Um fluxo de trabalho alternativo seria
1) Execute o seu script primeiro, obtenha uma estimativa para o intervalo de leitura
2) Subtraia o arquivado [8] com o valor de readspan-2 * read-length
3) Execute novamente um script ajustado e obtenha uma estimativa mais precisa para SD

Ou só preciso ler mais sobre estatísticas básicas e SD :)

Obrigado, acabei de usar isso, muito útil!

Obrigado por este post.
Eu tenho uma pergunta: por que com diferentes subamostras do mesmo bam posso ter tantas diferenças no intervalo de leitura?
Por exemplo :

1) visualização do samtools accept_hits_chr22.bam | head -n 100000 | python getinsertsize.py -
Comprimento de leitura: média 101,0, STD = 0,0
Extensão de leitura: média 236,150649351, STD = 88,6332394609
2) visualização do samtools accept_hits_chr22.bam | head -n 400000 | python getinsertsize.py -
Comprimento de leitura: média 101,0, STD = 0,0
Extensão de leitura: média 221,92447156, STD = 79,0142612548
3) visualização do samtools accept_hits_chr22.bam | head -n 800000 | python getinsertsize.py -
Comprimento de leitura: média 101,0, STD = 0,0
Extensão de leitura: média 3383,2321526, STD = 3460,01998406
4) visualização do samtools accept_hits_chr22.bam | head -n 1000000 | python getinsertsize.py -
Comprimento de leitura: média 101,0, STD = 0,0
Extensão de leitura: média 3111,18687862, STD = 3442,69889313

Portanto, isso significa que há algumas leituras emparelhadas realmente longas que atrapalham o cálculo do intervalo de leitura. Você pode simplesmente usar os resultados das 100 mil leituras principais.

Este comentário foi removido pelo autor.

Estou tentando seu script no formato Sam. & gt igave seguinte comando

python getinsertsize.py xyz.sam | -

ele mostrou um erro na linha 16 # nenhum módulo argparse

O comando acima precisa de qualquer modificação ou comando whats para executar seu código no formato sam. Eu tentei com a opção de visualização em samtools, mas é para o formato bam não para o sam. Por favor, você pode me orientar Como posso executar para o formato sam para calcular a média, ler o comprimento e SD do arquivo de alinhamento sam como entrada

Olá, isso significa que não há módulo argparse em seu sistema Python. Qual versão do python você usou? O Python 2.7 não deve ter problemas para executar este script.

Obrigado, irei atualizar o Python. Meu python versão 2.6.6 Obrigado pela resposta.

Eu tentei atualizar o Python 2.7 agora. Ele foi executado completamente, mas, finalmente, em vez de imprimir a saída. Ele mostrou o seguinte erro.

Falha no fechamento do destruidor do objeto de arquivo:
sys.excepthook está faltando
perdeu sys.stderr

Você poderia me sugerir por que esse erro apareceu

Muito obrigado Wei Li, agora está funcionando perfeitamente

Oi Wei,
Muito obrigado por isso. você saberia como adaptar o script para que se pudesse selecionar ou excluir leituras específicas com base no comprimento do intervalo? Ex: digamos que você queira apenas ler pares com span & lt500.
obrigado!

Eu estou mapeando algumas leituras de Illumina para uma espécie intimamente relacionada usando BWA e obtendo alguns resultados estranhos. Meu intervalo médio de leitura é diferente de zero, mas é menor do que o comprimento médio de leitura. Alguma ideia do que pode causar isso? Sou bastante novo nisso e não conheço python.

Comprimento de leitura: média 94,7053421448, STD = 15,8237819615
Extensão de leitura: média 74,962962963, STD = 48,8084808847

Provavelmente, seu arquivo bam inclui uma mistura de leituras com comprimentos diferentes. Isso pode complicar o cálculo do intervalo de leitura.

Eu tenho exatamente a mesma questão.

Legal..blog você também pode visitar nosso website que é um software de estimativa.http: //www.accurateestimator.com/

Alguma melhoria para leituras de comprimento de leitura desiguais?
Tenho problemas com diferentes bibliotecas (MP / PE), mas os tamanhos de inserção são todos

150bp! Não tenho certeza se isso & # 39s devido ao comprimento de leitura irregular.

Obrigado por boas ferramentas de trabalho
Eu tenho esta saída: Comprimento de leitura de leituras de extremidades emparelhadas: média 296,22, STD = 22,51
Extensão de leitura: média 484,54, STD = 108,36
Isso significa que o tamanho da minha inserção é 484-296 = 184?


Como recuperar uma ou mais sequências biológicas?

Uma das operações mais frequentes que faço durante o dia é a busca e recuperação de sequências biológicas de interesse. Portanto, acho que é necessário conversar com você sobre como isso pode ser feito, especialmente para novatos, e, portanto, ser capaz de & # 8220 jogar & # 8221 com alguns genes ou proteínas recuperados. Lembro-me que no início dos meus estudos de bioinformática gostava de recuperar sequências de genes e fiquei encantado ao observar aquela multidão de letras (relativas às bases nitrogenadas) que se seguiam dentro da sequência. Resumindo, obter uma sequência biológica de interesse é uma habilidade básica para um bioinformático e com este artigo explicarei como fazê-lo.

Vamos começar com uma suposição simples. As sequências biológicas de interesse podem ser recuperadas pesquisando em bancos de dados genéricos ou dedicados à espécie em questão ou usando comandos específicos lançados do terminal para obter a sequência desejada de arquivos baixados anteriormente. Mas vamos fazer alguns exemplos práticos imediatamente para entender melhor como fazê-lo.

RECUPERAR SEQUÊNCIAS DA BASE DE DADOS

Em relação à recuperação de uma ou mais sequências de um banco de dados, duas premissas importantes devem ser feitas:

  1. Existem diferentes tipos de bases de dados, mais ou menos especializadas na recuperação de determinadas sequências, sendo possível encontrar na web bases de dados específicas para algumas espécies ou famílias, como a base de dados de sequências da família Solanaceae, ou seja, Sol Genomics Network.
  2. Para agilizar e otimizar a busca no banco de dados é necessário o uso de alguns operadores booleanos como:
    • E, permite que você conecte duas ou mais cadeias de caracteres de pesquisa. Seu significado lógico pode ser entendido como o cruzamento entre duas consultas que nos interessa pesquisar de forma combinada.
    • OU, permite que você pesquise um elemento discriminando outros quando mais elementos são procurados por vez.
    • NOT, permite que você pesquise todos os elementos que não correspondam a uma determinada consulta. Portanto, ele tem uma função oposta ao operador AND.

Mas vamos fazer alguns exemplos:

Quando você deseja pesquisar algo em um banco de dados, precisa saber em qual banco de dados é mais fácil encontrar uma sequência específica. Deixe-me explicar melhor. Se você deseja pesquisar uma proteína específica, é útil consultar e depois pesquisar no banco de dados UniProt, que é gerenciado por um consórcio privado, ou no banco de dados Protein gerenciado pelo NCBI, certamente não faria sentido pesquisar para a proteína acima mencionada no banco de dados Gene. Certo? Ou ainda, se sua intenção é baixar o genoma do tomate, é mais fácil encontrá-lo no banco de dados do Genoma do que no banco de dados Assembly, onde encontramos principalmente informações e estatísticas relacionadas à montagem desse genoma. Em suma, a mensagem para levar para casa a respeito da pesquisa de sequências ou dados nos bancos de dados é:

Pesquise o que você precisa no banco de dados certo para ter mais sucesso em sua pesquisa.

Na verdade, ninguém jamais sairia e compraria um bife no verdureiro e uma maçã no açougueiro ... Você entende o que quero dizer?

RECUPERE SEQUÊNCIAS DO ARQUIVO USANDO COMANDOS EXECUTADOS NO TERMINAL

Em alguns casos, encontramos em nossas mãos um arquivo contendo diferentes informações e sequências, mas queremos derivar apenas algumas informações ou sequências a partir dele. Para isso contamos com a ajuda de comandos que, se executados no terminal, nos permitem chegar facilmente ao nosso objetivo.

Vamos apresentar algumas situações ideais para entender melhor como podemos fazer o acima.

  1. Baixamos o genoma de Arabidopsis thaliana do banco de dados Phytozome e recuperamos a partir dele a sequência & # 8220TGTAGGGATGAAGTCTTTCTTCGTTGT & # 8221. Como podemos fazer isso? grep pode ser usado, mas você deve levar em consideração que o grep permite que você extrapole apenas a sequência fornecida como uma consulta e a linha na qual ela é encontrada.

No vídeo abaixo, mostrei como usar essas três ferramentas:

Em relação ao grep, haveria muito mais a dizer, mas prefiro não fazer isso agora porque me afastaria do verdadeiro objetivo deste artigo, então prometerei apenas falar sobre isso em profundidade em um artigo futuro.

  1. Agora, vamos considerar outra situação possível. Digamos que sempre temos o arquivo do genoma de Arabidopsis thaliana e queremos extrapolar todo o cromossomo 3 dele. Para fazer isso, é útil usar o samtools faidx ferramenta. Assista ao vídeo abaixo para entender como:
  1. A situação apresentada no ponto 2 é muito frequente no trabalho de um bioinformático, mas como podemos derivar mais de uma sequência de um determinado genoma? Digamos que temos o conjunto de genes do tomate e queremos extrapolar dele todos os genes que estão envolvidos de alguma forma na síntese e na atividade da rubisco. No vídeo abaixo eu mostro como:
  1. Vejamos outro exemplo de situação que pode ser bastante comum. Dado o genoma de Solanum tuberosum, encontre os cabeçalhos de estrutura 1 a 3 e extrapole as sequências. Para fazer isso é conveniente usar expressões regulares, veja como no vídeo abaixo:

Também neste caso, o uso de expressões regulares deve ser aprofundado, então, como no caso do grep, tentarei falar melhor sobre isso em um artigo dedicado no futuro.

  1. Outra situação muito comum para um bioinformático. Temos para nossas & # 8220hands & # 8221 um arquivo vcf (formato de chamada de variante) que é um arquivo que mostra em uma matriz, consistindo de linhas e colunas, as várias informações (especialmente posicionais) relacionadas aos polimorfismos rastreados dentro de um determinado genoma. Às vezes pode ser útil extrapolar a partir deste arquivo apenas informações específicas e isso pode ser feito com o comando awk (também irei falar sobre isso em um artigo separado para não aborrecê-lo muito) ou através do uso de nossas queridas expressões regulares. No vídeo abaixo, mostrei um exemplo de awk.
  1. Finalmente, é útil dizer que com samtools faidx também é possível extrapolar uma sequência de um arquivo fasta graças às coordenadas deste último e, portanto, à sua posição no genoma considerado, visto que esta ferramenta é capaz de anotar o genoma. antes de procurar a sequência.

Aqui estamos nós no final deste artigo. Como sempre, espero ter conseguido explicar conceitos complexos da maneira mais simples possível. Deixe-me saber com um comentário o que você achou do artigo ou se há algo que você acha que deve ser adicionado ou corrigido.


Conteúdo

Um arquivo FASTQ normalmente usa quatro linhas por sequência.

  • A linha 1 começa com um caractere '@' e é seguida por um identificador de sequência e um opcional descrição (como uma linha de título FASTA).
  • A linha 2 é a sequência de letras brutas.
  • A linha 3 começa com um caractere '+' e é opcionalmente seguido pelo mesmo identificador de sequência (e qualquer descrição) novamente.
  • A linha 4 codifica os valores de qualidade para a sequência na linha 2 e deve conter o mesmo número de símbolos que letras na sequência.

Um arquivo FASTQ contendo uma única sequência pode ter a seguinte aparência:

O byte que representa a qualidade vai de 0x21 (qualidade mais baixa '!' Em ASCII) a 0x7e (qualidade mais alta '

'em ASCII). Aqui estão os caracteres de valor de qualidade em ordem crescente de qualidade da esquerda para a direita (ASCII):

Os arquivos Sanger FASTQ originais também permitiam que as strings de sequência e qualidade fossem quebradas (divididas em várias linhas), mas isso geralmente é desencorajado [ citação necessária ], pois pode tornar a análise complicada devido à escolha infeliz de "@" e "+" como marcadores (esses caracteres também podem ocorrer na string de qualidade).

Identificadores de sequência Illumina Editar

As sequências do software Illumina usam um identificador sistemático:

HWUSI-EAS100R o nome único do instrumento
6 pista da célula de fluxo
73 número do ladrilho dentro da pista da célula de fluxo
941 coordenada 'x' do cluster dentro do bloco
1973 coordenada 'y' do cluster dentro do bloco
#0 número de índice para uma amostra multiplexada (0 para nenhuma indexação)
/1 o membro de um par, / 1 ou / 2 (apenas leituras emparelhadas ou pares)

Versões do pipeline Illumina desde 1.4 parecem usar #NNNNNN ao invés de #0 para o ID multiplex, onde NNNNNN é a sequência da etiqueta multiplex.

Com Casava 1.8, o formato da linha '@' mudou:

EAS139 o nome único do instrumento
136 o ID de execução
FC706VJ a identificação da célula de fluxo
2 pista da célula de fluxo
2104 número do ladrilho dentro da pista da célula de fluxo
15343 coordenada 'x' do cluster dentro do bloco
197393 coordenada 'y' do cluster dentro do bloco
1 o membro de um par, 1 ou 2 (apenas leituras emparelhadas ou pares)
Y Y se a leitura for filtrada (não passou), N caso contrário
18 0 quando nenhum dos bits de controle está ligado, caso contrário, é um número par
ATCACG sequência de índice

Observe que as versões mais recentes do software Illumina geram um número de amostra (conforme retirado da folha de amostra) no lugar de uma sequência de índice. Por exemplo, o seguinte cabeçalho pode aparecer na primeira amostra de um lote:

Edição de arquivo de leitura de sequência NCBI

Os arquivos FASTQ do arquivo de leitura de sequência INSDC geralmente incluem uma descrição, por exemplo,

Neste exemplo, há um identificador atribuído ao NCBI e a descrição contém o identificador original da Solexa / Illumina (conforme descrito acima) mais o comprimento de leitura. O sequenciamento foi realizado em modo pareado (

Tamanho do inserto de 500 bp), consulte SRR001666. O formato de saída padrão do fastq-dump produz pontos inteiros, contendo quaisquer leituras técnicas e, normalmente, leituras biológicas de extremidade única ou pareada.

Modern usage of FASTQ almost always involves splitting the spot into its biological reads, as described in submitter-provided metadata:

When present in the archive, fastq-dump can attempt to restore read names to original format. NCBI does not store original read names by default:

In the example above, the original read names were used rather than the accessioned read name. NCBI accessions runs and the reads they contain. Original read names, assigned by sequencers, are able to function as locally unique identifiers of a read, and convey exactly as much information as a serial number. The ids above were algorithmically assigned based upon run information and geometric coordinates. Early SRA loaders parsed these ids and stored their decomposed components internally. NCBI stopped recording read names because they are frequently modified from the vendors' original format in order to associate some additional information meaningful to a particular processing pipeline, and this caused name format violations that resulted in a high number of rejected submissions. Without a clear schema for read names, their function remains that of a unique read id, conveying the same amount of information as a read serial number. See various SRA Toolkit issues for details and discussions.

Also note that fastq-dump converts this FASTQ data from the original Solexa/Illumina encoding to the Sanger standard (see encodings below). This is because the SRA serves as a repository for NGS information, rather than format. The various *-dump tools are capable of producing data in several formats from the same source. The requirements for doing so have been dictated by users over several years, with the majority of early demand coming from the 1000 Genomes Project.

Quality Edit

A quality value Q is an integer mapping of p (i.e., the probability that the corresponding base call is incorrect). Two different equations have been in use. The first is the standard Sanger variant to assess reliability of a base call, otherwise known as Phred quality score:

The Solexa pipeline (i.e., the software delivered with the Illumina Genome Analyzer) earlier used a different mapping, encoding the odds p/(1-p) instead of the probability p:

Although both mappings are asymptotically identical at higher quality values, they differ at lower quality levels (i.e., approximately p > 0.05, or equivalently, Q < 13).

At times there has been disagreement about which mapping Illumina actually uses. The user guide (Appendix B, page 122) for version 1.4 of the Illumina pipeline states that: "The scores are defined as Q=10*log10(p/(1-p)) [sic], where p is the probability of a base call corresponding to the base in question". [2] In retrospect, this entry in the manual appears to have been an error. The user guide (What's New, page 5) for version 1.5 of the Illumina pipeline lists this description instead: "Important Changes in Pipeline v1.3 [sic] The quality scoring scheme has changed to the Phred [i.e., Sanger] scoring scheme, encoded as an ASCII character by adding 64 to the Phred value. A Phred score of a base is: Q phred = − 10 log 10 ⁡ e >=-10log _< ext<10>>e> , where e is the estimated probability of a base being wrong. [3]

Edição de codificação

  • Sanger format can encode a Phred quality score from 0 to 93 using ASCII 33 to 126 (although in raw read data the Phred quality score rarely exceeds 60, higher scores are possible in assemblies or read maps). Also used in SAM format. [4] Coming to the end of February 2011, Illumina's newest version (1.8) of their pipeline CASAVA will directly produce fastq in Sanger format, according to the announcement on seqanswers.com forum. [5]
  • PacBio HiFi reads, which are typically stored in SAM/BAM format, use the Sanger convention: Phred quality scores from 0 to 93 are encoded using ASCII 33 to 126. Raw PacBio subreads use the same convention but typically assign a placeholder base quality (Q0) to all bases in the read. [6]
  • Solexa/Illumina 1.0 format can encode a Solexa/Illumina quality score from -5 to 62 using ASCII 59 to 126 (although in raw read data Solexa scores from -5 to 40 only are expected)
  • Starting with Illumina 1.3 and before Illumina 1.8, the format encoded a Phred quality score from 0 to 62 using ASCII 64 to 126 (although in raw read data Phred scores from 0 to 40 only are expected).
  • Starting in Illumina 1.5 and before Illumina 1.8, the Phred scores 0 to 2 have a slightly different meaning. The values 0 and 1 are no longer used and the value 2, encoded by ASCII 66 "B", is used also at the end of reads as a Read Segment Quality Control Indicator. [7] The Illumina manual [8] (page 30) states the following: If a read ends with a segment of mostly low quality (Q15 or below), then all of the quality values in the segment are replaced with a value of 2 (encoded as the letter B in Illumina's text-based encoding of quality scores). This Q2 indicator does not predict a specific error rate, but rather indicates that a specific final portion of the read should not be used in further analyses. Also, the quality score encoded as "B" letter may occur internally within reads at least as late as pipeline version 1.6, as shown in the following example:

An alternative interpretation of this ASCII encoding has been proposed. [9] Also, in Illumina runs using PhiX controls, the character 'B' was observed to represent an "unknown quality score". The error rate of 'B' reads was roughly 3 phred scores lower the mean observed score of a given run.

  • Starting in Illumina 1.8, the quality scores have basically returned to the use of the Sanger format (Phred+33).

For raw reads, the range of scores will depend on the technology and the base caller used, but will typically be up to 41 for recent Illumina chemistry. Since the maximum observed quality score was previously only 40, various scripts and tools break when they encounter data with quality values larger than 40. For processed reads, scores may be even higher. For example, quality values of 45 are observed in reads from Illumina's Long Read Sequencing Service (previously Moleculo).

Color space Edit

For SOLiD data, the sequence is in color space, except the first position. The quality values are those of the Sanger format. Alignment tools differ in their preferred version of the quality values: some include a quality score (set to 0, i.e. '!') for the leading nucleotide, others do not. The sequence read archive includes this quality score.

Simulation Edit

FASTQ read simulation has been approached by several tools. [10] [11] A comparison of those tools can be seen here. [12]

Compression Edit

General compressors Edit

General-purpose tools such as Gzip and bzip2 regard FASTQ as a plain text file and result in suboptimal compression ratios. NCBI's Sequence Read Archive encodes metadata using the LZ-77 scheme. General FASTQ compressors typically compress distinct fields (read names, sequences, comments, and quality scores) in a FASTQ file separately these include Genozip, [13] DSRC and DSRC2, FQC, LFQC, Fqzcomp, and Slimfastq.

Reads Edit

Having a reference genome around is convenient because then instead of storing the nucleotide sequences themselves, one can just align the reads to the reference genome and store the positions (pointers) and mismatches the pointers can then be sorted according to their order in the reference sequence and encoded, e.g., with run-length encoding. When the coverage or the repeat content of the sequenced genome is high, this leads to a high compression ratio. Unlike the SAM/BAM formats, FASTQ files do not specify a reference genome. Alignment-based FASTQ compressors supports the use of either user-provided or de novo assembled reference: LW-FQZip uses a provided reference genome and Quip, Leon, k-Path and KIC perform de novo assembly using a de Bruijn graph-based approach. Genozip [13] can optionally use a reference if the user provides one, which may be a single- or multi-species reference file.

Explicit read mapping and de novo assembly are typically slow. Reordering-based FASTQ compressors first cluster reads that share long substrings and then independently compress reads in each cluster after reordering them or assembling them into longer contigs, achieving perhaps the best trade-off between the running time and compression rate. SCALCE is the first such tool, followed by Orcom and Mince. BEETL uses a generalized Burrows–Wheeler transform for reordering reads, and HARC achieves better performance with hash-based reordering. AssemblTrie instead assembles reads into reference trees with as few total number of symbols as possible in the reference. [14] [15]

Benchmarks for these tools are available in. [16]

Quality values Edit

Quality values account for about half of the required disk space in the FASTQ format (before compression), and therefore the compression of the quality values can significantly reduce storage requirements and speed up analysis and transmission of sequencing data. Both lossless and lossy compression are recently being considered in the literature. For example, the algorithm QualComp [17] performs lossy compression with a rate (number of bits per quality value) specified by the user. Based on rate-distortion theory results, it allocates the number of bits so as to minimize the MSE (mean squared error) between the original (uncompressed) and the reconstructed (after compression) quality values. Other algorithms for compression of quality values include SCALCE [18] and Fastqz. [19] Both are lossless compression algorithms that provide an optional controlled lossy transformation approach. For example, SCALCE reduces the alphabet size based on the observation that “neighboring” quality values are similar in general. For a benchmark, see. [20]

As of the HiSeq 2500 Illumina gives the option to output qualities that have been coarse grained into quality bins. The binned scores are computed directly from the empirical quality score table, which is itself tied to the hardware, software and chemistry that were used during the sequencing experiment. [21]

Genozip [13] uses its DomQual algorithm to compress binned quality scores, such as those generated by Illumina or by Genozip's own --optimize option which generates bins similar to Illumina.

Encryption Edit

Genozip [13] encrypts FASTQ files (as well as other genomic formats), by applying the standard AES encryption at its most secure level of 256 bits (--password option).

Cryfa [22] uses AES encryption and enables to compact data besides encryption. It can also address FASTA files.

There is no standard file extension for a FASTQ file, but .fq and .fastq are commonly used.


What is Numerical Taxonomy? How is it useful?

Classification of biological species is one of the important concern while studying taxonomic and or evolutionary relationships among various species. Classification is either based on only one / a few characters known as “Monothetic”, or based on multiple characters known as “Polythetic”.

It is obviously much more difficult to classify organisms on the basis of multiple characters rather than a few characters. The traditional approaches of taxonomists are tedious. The arrival of computer techniques in the field of biology has made the task easier for the taxonomists.

Numerical taxonomy is a system of grouping of species by numerical methods based on their character states. It was first initiated by Peter H.A.Sneath et al.

Before going further I would like to clear the difference between two common terms, namely, “Classification” & “Identification”. When the organisms are classified on the basis of like properties, then it is called Classification, and after the classification, when the additional unidentified objects are allocated, then it is known as Identification. The purpose of taxonomy is to group the objects to be classified into natural taxa. Conventional taxonomists equate the taxonomic relationships with evolutionary relationships, but the numerical taxonomists defined them as three kinds:

  • Phenetic: based on overall similarity.
  • Cladistic: based on a common line of descents.
  • Chronistic: temporal relation among various evolutionary branches.

How is the classification do by Numerical Taxonomy?

The objects to be classified are known as Operational Taxonomic Units (OTUs). They may be species, genera, family, higher ranking taxonomic groups,etc., The characters are numerically recorded either in the form of appropriate numbers or may be programmed in such a way that the differences among them are proportional to their dissimilarity. Lets say, a character called ‘hairness of leaf’, it may be recorded as:

  • hairless = 0
  • sparsely haired = 1
  • regularly haired = 2
  • densely haired = 3

Fig.1 OTUs (black dots) represented in a multidimensional space.

Such a numerical system imply that the dissimilarity between densely haired and hairless is 3 times than that of sparsely haired and hairless.

The other method of implementing numerical taxonomy is that the characters are always represented by only two states,i.e., 0 for the absence and 1 for the presence of a particular character. This method is usually implemented in the field of microbiology. After that, all the characters and the taxonomic units are arranged in the form of the data matrix and the similarity among all possible pairs of OTUs is computed based on the considered characters. The similarity ( more specifically, dissimilarity) is the distance between OTUs is represented in a multidimensional space, where the characters can be visualized as the coordinates. The objects that are very similar are plotted close to each other and those which are dissimilar are plotted farther apart. Then these straight lines are computed. The similarity among the OTUs is calculated by ‘simialrity matrix’ having few color schemes, where the dark-shaded areas are highly similar. This matrix is then rearranged to get the clusters of similar OTUs. The results of numerical taxonomy are generally represented in the form of phenograms.

An exhaustive list of references for this article is available with the author and is available on personal request, for more details write to [email protected]


Our 2020/2021 Highlights document

Over the past year we have provided new insights into human, parasite and microbe evolution, cellular development, and the mutational processes that lead to cancer. We supported the UK's national COVID-19 genomic surveillance efforts and delivered reference genomes for a wide range of UK species as part of the Darwin Tree of Life project . mais

Download our 2020/21 Highlights document

Watch Again: Dr Roser Vento-Tormo's online seminar and Q&A

Mapping tissues na Vivo e em vitro one cell at a time:

Tissue microenvironment shapes cellular identity and function, understanding how this happens na Vivo can help us engineer novel em vitro models.

Dr Roser Vento-Tormo’s research is focused on the adaptation of immune cells in tissues and their function in steady state and inflammation. Her team uses genomics, spatial transcriptomics and bioinformatics tools to reconstruct the microenvironment that will shape immune cellular identity and function. Watch her seminar and Q&A at the link below:


Is it possible to re-align of SAM/BAM files? - Biologia

Run stringtie from the command line like this:
The main input of the program is a BAM file with RNA-Seq read mappings which must be sorted by their genomic location (for example the accepted_hits.bam file produced by TopHat or the output of HISAT2 after sorting and converting it using samtools as explained below).

Input files

The file resulted from the above command (alns.sorted.bam) can be used as input for StringTie.

Every spliced read alignment (i.e. an alignment across at least one junction) in the input SAM file must contain the tag XS to indicate the genomic strand that produced the RNA from which the read was sequenced. Alignments produced by TopHat and HISAT2 (when run with --dta option) already include this tag, but if you use a different read mapper you should check that this XS tag is included for spliced alignments.

NOTE: be sure to run HISAT2 with the --dta option for alignment, or your results will suffer.

As an option, a reference annotation file in GTF/GFF3 format can be provided to StringTie. In this case, StringTie will prefer to use these "known" genes from the annotation file, and for the ones that are expressed it will compute coverage, TPM and FPKM values. It will also produce additional transcripts to account for RNA-seq data that aren't covered by (or explained by) the annotation. Note that if option -e is not used the reference transcripts need to be fully covered by reads in order to be included in StringTie's output. In that case, other transcripts assembled from the data by StringTie and not present in the reference file will be printed as well.

NOTE: we highly recommend that you provide annotation if you are analyzing a genome that is well-annotated, such as human, mouse, or other model organisms.

Output files

  1. Stringtie's main output is a GTF file containing the assembled transcripts
  2. Gene abundances in tab-delimited format
  3. Fully covered transcripts that match the reference annotation, in GTF format
  4. Files (tables) required as input to Ballgown, which uses them to estimate differential expression
  5. In merge mode, a merged GTF file from a set of GTF files

Description of each column's values:

  • seqname: Denotes the chromosome, contig, or scaffold for this transcript. Here the assembled transcript is on chromosome X.
  • source: The source of the GTF file. Since this example was produced by StringTie, this column simply shows 'StringTie'.
  • feature: Feature type e.g., exon, transcript, mRNA, 5'UTR).
  • start: Start position of the feature (exon, transcript, etc), using a 1-based index.
  • fim: End position of the feature, using a 1-based index.
  • pontuação: A confidence score for the assembled transcript. Currently this field is not used, and StringTie reports a constant value of 1000 if the transcript has a connection to a read alignment bundle.
  • strand: If the transcript resides on the forward strand, '+'. If the transcript resides on the reverse strand, '-'.
  • frame: Frame or phase of CDS features. StringTie does not use this field and simply records a ".".
  • attributes: A semicolon-separated list of tag-value pairs, providing additional information about each feature. Depending on whether an instance is a transcript or an exon and on whether the transcript matches the reference annotation file provided by the user, the content of the attributes field will differ. The following list describes the possible attributes shown in this column:
    • gene_id: A unique identifier for a single gene and its child transcript and exons based on the alignments' file name.
    • transcript_id: A unique identifier for a single transcript and its child exons based on the alignments' file name.
    • exon_number: A unique identifier for a single exon, starting from 1, within a given transcript.
    • reference_id: The transcript_id in the reference annotation (optional) that the instance matched.
    • ref_gene_id: The gene_id in the reference annotation (optional) that the instance matched.
    • ref_gene_name: The gene_name in the reference annotation (optional) that the instance matched.
    • cov: The average per-base coverage for the transcript or exon.
    • FPKM: Fragments per kilobase of transcript per million read pairs. This is the number of pairs of reads aligning to this feature, normalized by the total number of fragments sequenced (in millions) and the length of the transcript (in kilobases).
    • TPM: Transcripts per million. This is the number of transcripts from this particular gene normalized first by gene length, and then by sequencing depth (in millions) in the sample. A detailed explanation and a comparison of TPM and FPKM can be found here, and TPM was defined by B. Li and C. Dewey here.
    • Column 1 / Gene ID: The gene identifier comes from the reference annotation provided with the -G option. If no reference is provided this field is replaced with the name prefix for output transcripts (-l).
    • Column 2 / Gene Name: This field contains the gene name in the reference annotation provided with the -G option. If no reference is provided this field is populated with '-'.
    • Column 3 / Referência: Name of the reference sequence that was used in the alignment of the reads. Equivalent to the 3rd column in the .SAM alignment.
    • Column 4 / Strand: '+' denotes that the gene is on the forward strand, '-' for the reverse strand.
    • Column 5 / Start: Start position of the gene (1-based index).
    • Column 6 / End: End position of the gene (1-based index).
    • Column 7 / Coverage: Per-base coverage of the gene.
    • Column 8 / FPKM: normalized expression level in FPKM units (see previous section).
    • Column 9 / TPM: normalized expression level in RPM units (see previous section).

    3. Fully covered transcripts matching the reference annotation transcripts (in GTF format)

    If StringTie is run with the -C <cov_refs.gtf> option (requires -G <reference_annotation> ), it returns a file with all the transcripts in the reference annotation that are fully covered, end to end, by reads. The output format is a GTF file as described above. Each line of the GTF is corresponds to a gene or transcript in the reference annotation.

    4. Ballgown Input Table Files

    If StringTie is run with the -B option, it returns a Ballgown input table file, which contains coverage data for all transcripts. The output table files are placed in the same directory as the main GTF output. These tables have these specific names: (1) e2t.ctab, (2) e_data.ctab, (3) i2t.ctab, (4) i_data.ctab, and (5) t_data.ctab. A detailed description of each of these five required inputs to Ballgown can be found on the Ballgown site, at this link.

    5. Merge mode: Merged GTF

    If StringTie is run with the --merge option, it takes as input a list of GTF/GFF files and merges/assembles these transcripts into a non-redundant set of transcripts. This step creates a uniform set of transcripts for all samples to facilitate the downstream calculation of differentially expressed levels for all transcripts among the different experimental conditions. Output is a merged GTF file with all merged gene models, but without any numeric results on coverage, FPKM, and TPM. Then, with this merged GTF, StringTie can re-estimate abundances by running it again with the -e option on the original set of alignment files, as illustrated in the figure below.

    Evaluating transcript assemblies

    A simple way of getting more information about the transcripts assembled by StringTie (summary of gene and transcript counts, novel vs. known etc.), or even performing basic tracking of assembled isoforms across multiple RNA-Seq experiments, is to use the gffcompare program. Basic usage information and download options for this program can be found on the GFF utilities page.

    Differential expression analysis

    1. for each RNA-Seq sample, map the reads to the genome with HISAT2 using the --dta option. It is highly recommended to use the reference annotation information when mapping the reads, which can be either embedded in the genome index (built with the --ss e --exon options, see HISAT2 manual), or provided separately at run time (using the --known-splicesite-infile option of HISAT2). The SAM output of each HISAT2 run must be sorted and converted to BAM using samtools as explained above.
    2. for each RNA-Seq sample, run StringTie to assemble the read alignments obtained in the previous step it is recommended to run StringTie with the -G option if the reference annotation is available.
    3. run StringTie with --merge in order to generate a non-redundant set of transcripts observed in any of the RNA-Seq samples assembled previously. o stringtie --merge mode takes as input a list of all the assembled transcripts files (in GTF format) previously obtained for each sample, as well as a reference annotation file (-G option) if available.
    4. for each RNA-Seq sample, run StringTie using the -B/-b e -e options in order to estimate transcript abundances and generate read coverage tables for Ballgown. The -e option is not required but recommended for this run in order to produce more accurate abundance estimations of the input transcripts. Each StringTie run in this step will take as input the sorted read alignments (BAM file) obtained in step 1 for the corresponding sample and the -G option with the merged transcripts (GTF file) generated by stringtie --merge in step 3. Please note that this is the only case where the -G option is not used with a reference annotation, but with the global, merged set of transcripts as observed across all samples. (This step is the equivalent of the Tablemaker step described in the original Ballgown pipeline.)
    5. Ballgown can now be used to load the coverage tables generated in the previous step and perform various statistical analyses for differential expression, generate plots etc.

    An alternate, faster differential expression analysis workflow can be pursued if there is no interest in novel isoforms (i.e. assembled transcripts present in the samples but missing from the reference annotation), or if only a well known set of transcripts of interest are targeted by the analysis. This simplified protocol has only 3 steps (depicted below) as it bypasses the individual assembly of each RNA-Seq sample and the "transcript merge" step. This simplified workflow attempts to directly estimate and analyze the expression of a known set of transcripts as given in the reference annotation Arquivo.

    O pacote R IsoformSwitchAnalyzeR can be used to assign gene names to transcripts assembled by StringTie, which can be particularly helpful in cases where StringTie could not perform this assignment unambigiously.
    The importIsoformExpression() + importRdata() function of the package can be used to import the expression and annotation data into R. During this import the package will attempt to clean up and recover isoform annotations where possible. From the resulting switchAnalyzeRlist object, IsoformSwitchAnalyzeR can detect isoform switches along with predicted functional consequences. The extractGeneExpression() function can be used to get a gene expression (read count) matrix for analysis with other tools.
    More information and code examples can be found here.

    Using StringTie with DESeq2 and edgeR

    DESeq2 and edgeR are two popular Bioconductor packages for analyzing differential expression, which take as input a matrix of read counts mapped to particular genomic features (e.g., genes). We provide a Python script (prepDE.py, or the Python 3 version: prepDE.py3 ) that can be used to extract this read count information directly from the files generated by StringTie (run with the -e parameter).

    prepDE.py derives hypothetical read counts for each transcript from the coverage values estimated by StringTie for each transcript, by using this simple formula: reads_per_transcript = coverage * transcript_len / read_len

    • one option is to provide a path to a directory containing all sample sub-directories, with the same structure as the ballgown directory in the StringTie protocol paper in preparation for Ballgown. By default (no -i option), the script is going to look in the current directory for all sub-directories having .gtf files in them, as in this example:
    • Alternatively, one can provide a text file listing sample IDs and their respective paths (sample_lst.txt).

    Usage: prepDE.py [options]
    generates two CSV files containing the count matrices for genes and transcripts, using the coverage values found in the output of stringtie -e

    Options:
    -h, --help show this help message and exit
    -i INPUT, --input=INPUT, --in=INPUTa folder containing all sample sub-directories, or a text file with sample ID and path to its GTF file on each line [default: . ]
    -g Gwhere to output the gene count matrix [default: gene_count_matrix.csv]
    -t Twhere to output the transcript count matrix [default: transcript_count_matrix.csv]
    -l LENGTH, --length=LENGTHthe average read length [default: 75]
    -p PATTERN, --pattern=PATTERNa regular expression that selects the sample subdirectories
    -c, --clusterwhether to cluster genes that overlap with different gene IDs, ignoring ones with geneID pattern (see below)
    -s STRING, --string=STRINGif a different prefix is used for geneIDs assigned by StringTie [default: MSTRG]
    -k KEY, --key=KEYif clustering, what prefix to use for geneIDs assigned by this script [default: prepG]
    --legend=LEGENDif clustering, where to output the legend file mapping transcripts to assigned geneIDs [default: legend.csv]

    These count matrices (CSV files) can then be imported into R for use by DESeq2 and edgeR (using the DESeqDataSetFromMatrix e DGEList functions, respectively).

    Protocol: Using StringTie with DESeq2

    Given a list of GTFs, which were re-estimated upon merging, users can follow the below protocol to use DESeq2 for differential expression analysis. To generate GTFs from raw reads follow the StringTie protocol paper (up to the Ballgown step).


    Variant calling with GATK

    Different variant callers disagree a great deal, for single nucleotide polymorphisms (SNPs) and particularly for insertions and deletions (indels). Of the various methods available (samtools, varscan, freebayes, ReadXplorer etc) GATK, by the Broad Institute is the best. The HaplotypeCaller module which performs local de novo assemblies around indels has recently been updated to include non-diploid organisms, so it’s used here but there is little difference for bacterial genomes at least, between HaplotypeCaller and UnifiedGenotyper.

    Following conversations at a few conferences and meetings now, most recently the brilliant 2nd ISCB Bioinformatics Symposium at TGAC, I’ve realised that others have had the same issues I had with implementing GATK. This mostly arises in preparation of your data beforehand, so here is a brief run-through of exactly how to call variants from microbial genome short read data using GATK which I hope is useful!:

    You’ll need to install samtools, picardtools, GATK, bwa and optionally, vcffilter for this workflow. Picardtools and GATK are simply .jar files so that’s no problem while you probably already have bwa installed, otherwise installation is well documented!

    This workflow begins with short read (fastq) files and a fasta reference. First the reference is prepared, a sequence dictionary is created, short reads are aligned to the reference and read group information provided, resulting sequence alignment map (sam) file sorted and converted to binary alignment map (bam) format, duplicates marked, bam file sorted, indel targets identified, indels realigned and variants called. Simple!

    For simplicity an example set of commands are provided here, where the reference is reference.fasta and the short reads are R1.fastq.gz and R2.fastq.gz. You will need to enter the paths and versions of the software being used at each step and your actual file names. Ploidy is set to 1.

    See Updated GATK pipeline to HaplotypeCaller + gVCF for implementation of the more recent gVCF workflow. Otherwise HaplotypeCaller alone (below) will work just fine


    #Index reference
    bwa index reference.fasta

    #Sort reference
    samtools faidx reference.fasta

    #Create sequence dictionary
    java -jar

    /bin/picard-tools-1.8.5/CreateSequenceDictionary.jar REFERENCE=reference.fasta OUTPUT=reference.dict

    #Align reads and assign read group
    bwa mem -R &#[email protected] ID:FLOWCELL1.LANE1 PL:ILLUMINA LB:test SM:PA01” reference.fasta R1.fastq.gz R2.fastq.gz > aln.sam

    #Sort sam file
    java -jar

    /bin/picard-tools-1.8.5/SortSam.jar I=aln.sam O=sorted.bam SORT_ORDER=coordinate

    #Mark duplicates
    java -jar

    /bin/picard-tools-version/MarkDuplicates.jar I=sorted.bam O=dedup.bam METRICS_FILE=metrics.txt

    #Sort bam file
    java -jar

    #Create realignment targets
    java -jar

    /bin/GATK3.3/GenomeAnalysisTK.jar -T RealignerTargetCreator -R reference.fasta -I dedup.bam -o targetintervals.list

    #Indel realignment
    java -jar

    /bin/GATK3.3/GenomeAnalysisTK.jar -T IndelRealigner -R reference.fasta -I dedup.bam -targetIntervals targetintervals.list -o realigned.bam

    #Call variants (HaplotypeCaller)
    java -jar

    /bin/GATK3.3/GenomeAnalysisTK.jar -T HaplotypeCaller -R reference.fasta -I realigned.bam -ploidy 1 -stand_call_conf 30 -stand_emit_conf 10 -o raw.vcf

    The resulting vcf file will contain your variant calls!

    Then you can optionally filter the variants:

    #Filter variants

    /bin/vcflib/bin/vcffilter -f ‘DP > 9’ -f ‘QUAL > 10’ raw.vcf > filtered.vcf

    Or first split the raw.vcf file into SNPs and indels:

    #Extract SNPs
    java -jar

    /bin/GATK3.3/GenomeAnalysisTK.jar -T SelectVariants -R reference.fasta -V raw.vcf -selectType SNP -o snps.vcf

    #Extract Indels
    java -jar

    /bin/GATK/GenomeAnalysisTK.jar -T SelectVariants -R reference.fasta -V raw.vcf -selectType INDEL -o indels.vcf

    I also have a neat perl wrapper to automate this workflow over many short read files and would be happy to make this available if people are interested. Please do comment with any questions or issues and I’ll do my best to resolve them!