samtools tview type consensus sequence to stdout

Question

[timmak]

Hi,

I found samtools tview quite amazing in that on the third line from the top, it produces a kind of consensus sequence using only what's available in the sorted .bam file. It does not require a reference sequence, and moreover, can identify regions with variation and mark these accordingly (e.g. with the appropriate nucleic acid notation such as Y, K, R, etc.). However, it's designed for interactive use and not so useful for batch queries.

Are there tools out there that can give me a consensus sequence for a particular region that is command-line friendly (i.e output to stdout)? I know samtools tview has a -d T option, but it only outputs 80 columns or so.

I have tried bcftools consensus and GATK FastaAlternateReferenceMaker but these tend to only give one consensus sequence masking the variation that may exist.

Many thanks for your help,

Tim

Answer 1

[Pierre Lindenbaum]

I know samtools tview has a -d T option, but it only outputs 80 columns or so.

you can override the ${COLUMNS} shell's variable.

$ COLUMNS=1000  samtools tview -d T in.bam

but beware for the memory required to print a large matrix.

[timmak]

samtools apparently resets the COLUMNS variable to 110 every time I run it. Which version are you using? I tried both 1.3.1 and 1.6, and it's the same.

$ COLUMNS=1000
[statgenpro ~/WORK/DCM/For Carol]
$ echo $COLUMNS
1000
[statgenpro ~/WORK/DCM/For Carol]
$ $samtoolsdir/samtools tview -p chr14:23898059 -d T $dir/$samp.sorted.dedup.realign.recal.bam | sed -n '3p'
CCACTGCCTTCCCATGTCTGGTCCACAGCTGGCTCTAAGCAAATAGCTGTTGAATGTGGGAGCGAGTGAGTGATTGTTCT
[statgenpro ~/WORK/DCM/For Carol]
$ echo $COLUMNS
110

[Pierre Lindenbaum]

this is not my cmd-line: I redefine/override COLUMNS each time I run samtools:

 $ COLUMNS=1000  samtools tview -d  ...

[timmak]

Thanks! That worked.