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

Question

[user3425664]

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!

[jv]

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

Answer 1

[user3425664]

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

Answer 2

[Pierre Lindenbaum]

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

[user3425664]

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)

[Pierre Lindenbaum]

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

this is not the command I wrote.

[user3425664]

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