get consensus sequence from bam file with insertions
0
1
Entering edit mode
3.4 years ago
Vvedenskiia ▴ 30

Hello everyone! I want to get consensus sequence from my NGS data. I created sorted bam, and visualized it in IGV. I know that alligment has SNPs, insertions and deletions. After that I tried to get consensus sequence from my NGS data. I used two ways:

  1. From samtools manual: "Generate the consensus sequence for one diploid individual: samtools mpileup -uf ref.fa aln.bam | bcftools call -c | vcfutils.pl vcf2fq > cns.fq" - this coomand gave me only SNPs in consensus sequence.
  2. Copy consensus sequence from IGV. It gave me consensus with SNPs and deletions.

How can I get insertions in consensus sequence?

Thanks!

consensus • 2.3k views
ADD COMMENT
0
Entering edit mode

Please use the formatting bar (especially the code option) to present your post better. You can use backticks for inline code (`text` becomes text), or select a chunk of text and use the highlighted button to format it as a code block. I've done it for you this time.
code_formatting

ADD REPLY
0
Entering edit mode
ADD REPLY
1
Entering edit mode

Yes, but in this try result sequence has same length as a reference. But region with reads is much less.

ADD REPLY
0
Entering edit mode

Could you please show a screenshot of a region where you suppose an indel?

ADD REPLY
1
Entering edit mode

More information. I used command:

bcftools consensus -f ref.fasta -s sorted_03.bam vars3_1.vcf.gz >out.fa

But insertions and deletions were not find in consensus.

ADD REPLY
0
Entering edit mode

This is part with deletion: https://ibb.co/27jVv4S

ADD REPLY
0
Entering edit mode

I saw this post: https://www.biostars.org/p/179357/ With sam=1.3, but than I tried

samtools mpileup -uf SA05sequence.fasta sorted.bam | bcftools call -c > vars.vcf

with bam file from bbmap it was not work. The report was: The sequence "chr3B chr3B:12400161..12411680" not found (this is my reference name). Maybe problems with indels can be associated with cigar sting? Different formats or something like that?

ADD REPLY
0
Entering edit mode

Did you resolve this problem?

ADD REPLY
0
Entering edit mode

Please be specific about what problem.

ADD REPLY

Login before adding your answer.

Traffic: 1340 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6