How Do I Generate Nucleotides Present Either Per Site or One I Choose?

I'm quite new in the world of bioinformatics and am self learning (no science degree). As the title of the post asks, how do I generate a list of every nucleotide that is present at any one given site? (a range would be great too) I'm looking for results similar to how IGV represents each site when you click on coverage.

My current commands are

bcftools mpileup -Ou -f SARS-CoV-2_reference.fasta {accession}.sorted.bam | bcftools call -mv -Ob -o {accession}.bcf
bcftools query -f  '{accession} %REF%POS%ALT\n'  {accession}.bcf

Thank you for your help!

bcftools mpileup

• 1,876 views • link

updated 3.1 years ago by Ram

• written 4.3 years ago by user3425664 •

I'm not quite sure it is that you are looking for, but I have you considered a bedgraph or wig file?

I prefer bedtools coverage for parsing out coverage info from a bam file, for info at https://bedtools.readthedocs.io/en/latest/content/tools/coverage.html

• 0 views • link

written 4.3 years ago by jv

2 answers

In case someone is here in the future and wants the same output I have found a python script that does exactly what I needed.

allelecount

• 0 views • link

written 4.3 years ago by user3425664 •

have a look at sam2tsv : http://lindenb.github.io/jvarkit/Sam2Tsv.html

samtools view -O BAM indexed.bam "chr1:234-235" | java -jar dist/sam2tsv -R ref.fasta

• 0 views • link

written 4.3 years ago by Pierre Lindenbaum

Thank you Pierre! I just tried it and i get this

samtools view SRR11409417.sorted.bam "NC_045512.2:29090-29096" | java -jar JavaKit/jvarkit/dist/sam2tsv.jar -R SARS-CoV-2_reference.fasta

[SEVERE][MultiBamLauncher]Dictionary missing : A Sequence dictionary is missing for <bam>. A Bam should have a header with a set of lines starting with '@SQ' see https://samtools.github.io/hts-specs/SAMv1.pdf com.github.lindenb.jvarkit.lang.JvarkitException$BamDictionaryMissing: Dictionary missing : A Sequence dictionary is missing for <bam>. A Bam should have a header with a set of lines starting with '@SQ' see https://samtools.github.io/hts-specs/SAMv1.pdf at com.github.lindenb.jvarkit.util.bio.SequenceDictionaryUtils.extractRequired(SequenceDictionaryUtils.java:166) at com.github.lindenb.jvarkit.jcommander.MultiBamLauncher.doWork(MultiBamLauncher.java:242) at com.github.lindenb.jvarkit.util.jcommander.Launcher.instanceMain(Launcher.java:796) at com.github.lindenb.jvarkit.util.jcommander.Launcher.instanceMainWithExit(Launcher.java:959) at com.github.lindenb.jvarkit.tools.sam2tsv.Sam2Tsv.main(Sam2Tsv.java:441) [INFO][Launcher]sam2tsv Exited with failure (-1)

• 0 views • link

written 4.3 years ago by user3425664 •

samtools view SRR11409417.sorted.bam "NC_045512.2:29090-29096"

this is not the command I wrote.

• 0 views • link

written 4.3 years ago by Pierre Lindenbaum

sorry I thought the O was a 0 so the command wasn't working :/ apologies

• 0 views • link

written 4.3 years ago by user3425664 •